diff --git a/Root/HiggsWeightTool.cxx b/Root/HiggsWeightTool.cxx
index e166a1a9b0030b849c1ebd5da8fc82a8eba83f3a..2f3e462f5389ee22fe7b792d12bf29b19635f38b 100644
--- a/Root/HiggsWeightTool.cxx
+++ b/Root/HiggsWeightTool.cxx
@@ -147,7 +147,11 @@ namespace xAOD {
return evtInfo->mcEventWeights();
}
- HiggsWeights HiggsWeightTool::getHiggsWeights(double pTH, int Njets30, int STXS_Stage1) {
+ /// Access MC weight for uncertainty propagation
+ /// Note: input kinematics should be HTXS_Higgs_pt, HTXS_Njets_pTjet30, and HTXS_Stage1_Category_pTjet30
+ HiggsWeights HiggsWeightTool::getHiggsWeights(int STXS_Njets30, double STXS_pTH, int STXS_Stage1) {
+ // convert to GeV
+ double pTH = STXS_pTH/1000; int Njets=STXS_Njets30;
const std::vector<float> &weights = getMCweights();
const xAOD::EventInfo *evtInfo = nullptr;
@@ -162,28 +166,64 @@ namespace xAOD {
HiggsWeights hw;
// set kinematics
- hw.pTH=pTH; hw.Njets30=Njets30;
+ hw.pTH=pTH; hw.Njets30=Njets;
// 1. Nominal weight
hw.nominal = weights[m_nom];
// 2. PDF weights
+ hw.alphaS_up=hw.alphaS_dn=0;
if (m_pdfUnc.size()==31) {
double pdf0=weights[m_pdfUnc[0]];
// PDF uncertainty. Scale to nominal weight (since PDF0 is relative to Powheg NLO)
// w_PDF[i] = PDF[i] / PDF[0] * w_nominal
for (size_t i=1;i<=30;++i)
hw.pdf4lhc.push_back(weights[m_pdfUnc[i]]/pdf0*hw.nominal);
- }
- // 2b. alphaS weights
- if (m_aS_up&&m_aS_dn) {
- hw.alphaS_up = weights[m_aS_up];
- hw.alphaS_dn = weights[m_aS_dn];
+ // 2b. alphaS weights
+ if (m_aS_up&&m_aS_dn) {
+ hw.alphaS_up = weights[m_aS_up]/pdf0*hw.nominal;
+ hw.alphaS_dn = weights[m_aS_dn]/pdf0*hw.nominal;
+ }
}
- // 3. Quark mass variations
- m_tinf=m_bminlo=m_nnlopsNom=0;
+ // nominal
+ hw.nnpdf30_nnlo = m_nnlopsNom ? weights[m_nnlopsNom] : 0;
+ double nnlo = hw.nnpdf30_nnlo;
+ // 3. Quark mass variations
+ //m_tinf=m_bminlo=m_nnlopsNom=0;
+ hw.mt_inf = m_tinf ? weights[m_tinf]*hw.nominal/nnlo : 0;
+ hw.mb_minlo = m_bminlo ? weights[m_bminlo]*hw.nominal/nnlo : 0;
+
+ // WG1 QCD uncertainty
+ // Cross sections in the =0, =1, and >=2 jets of Powheg ggH after reweighing scaled to sigma(N3LO)
+ static std::vector<double> sig({30.26,13.12,5.14});
+
+ // BLPTW absolute uncertainties in pb
+ static std::vector<double> yieldUnc({ 1.12, 0.66, 0.42});
+ static std::vector<double> resUnc ({ 0.03, 0.57, 0.42});
+ static std::vector<double> cut01Unc({-1.22, 1.00, 0.21});
+ static std::vector<double> cut12Unc({ 0,-0.86, 0.86});
+
+ // account for missing EW+quark mass effects by scaling BLPTW total cross section to sigma(N3LO)
+ double sf = 48.52/47.4;
+
+ int jetBin = (Njets > 1 ? 2 : Njets);
+ hw.qcd_wg1_mu = (1.0 + yieldUnc[jetBin]/sig[jetBin]*sf)*hw.nominal;
+ hw.qcd_wg1_res = (1.0 + resUnc[jetBin]/sig[jetBin]*sf)*hw.nominal;
+ hw.qcd_wg1_mig01 = (1.0 + cut01Unc[jetBin]/sig[jetBin]*sf)*hw.nominal;
+ hw.qcd_wg1_mig12 = (1.0 + cut12Unc[jetBin]/sig[jetBin]*sf)*hw.nominal;
+
+ // High pT uncertainty
+ static double y1_1 = 0.88, y2_1 = 1.16, x2_1 = 150;
+ static double y1_ge2 = 0.88, y2_ge2 = 1.16, x2_ge2 = 225;
+ double pTH_unc = 1.0;
+ if (Njets>=2) pTH_unc = pTH>x2_ge2?y2_ge2:y1_ge2+(y2_ge2-y1_ge2)*pTH/x2_ge2;
+ else if (Njets==1) pTH_unc = pTH>x2_1?y2_1:y1_1+(y2_1-y1_1)*pTH/x2_1;
+ hw.qcd_wg1_pTH = pTH_unc*hw.nominal;
+
+ // Quark-mass uncdertainty - TODO
+ hw.qcd_wg1_qm = hw.nominal;
return hw;
}
diff --git a/TruthWeightTools/HiggsWeightTool.h b/TruthWeightTools/HiggsWeightTool.h
index 15e2a63c488961602f8e4652969686fbd7476446..ffa720b852037396217e36378858766f99d79f32 100644
--- a/TruthWeightTools/HiggsWeightTool.h
+++ b/TruthWeightTools/HiggsWeightTool.h
@@ -46,10 +46,10 @@ namespace xAOD {
int Njets30, STXScat;
/// methods to print weights to the screen
- char *uncStr(double var, double nom) { return Form("%s%.3f",var>nom?"+":"",(var-nom)/nom); }
+ char *uncStr(double var, double nom) { return Form("%s%.1f%%",var>nom?"+":"",(var-nom)/nom*100); }
void print() {
double n=nominal;
- printf("\n Higgs MC weights of current event, pTH %.1f GeV, Njets=%i\n",
+ printf("\n------\n Higgs MC weights of current event, pTH = %.1f GeV, Njets = %i\n",
pTH,Njets30);
printf(" Nominal weight: %.3f\n",nominal);
if (pdf4lhc.size()==30) {
@@ -61,6 +61,13 @@ namespace xAOD {
for (size_t i=20;i<30;++i) printf(" %s",uncStr(pdf4lhc[i],n));
printf("\n alphaS up: %s, down: %s\n",
uncStr(alphaS_up,n),uncStr(alphaS_dn,n));
+ printf("\n Quark mass varations (m_top=inf): %s (m_b minlo): %s\n",
+ uncStr(mt_inf,n),uncStr(mb_minlo,n));
+ printf("\n WG1 proposed QCD uncertainty scheme\n");
+ printf(" mu: %s, res: %s, mig01: %s, mig12: %s\n",
+ uncStr(qcd_wg1_mu,n),uncStr(qcd_wg1_res,n),uncStr(qcd_wg1_mig01,n),uncStr(qcd_wg1_mig12,n));
+ printf(" pTH: %s, quark-mass: %s\n",
+ uncStr(qcd_wg1_pTH,n),uncStr(qcd_wg1_qm,n));
}
}
@@ -97,7 +104,7 @@ namespace xAOD {
/// Index of MC event weight
size_t getWeightIndex(std::string weightName);
- HiggsWeights getHiggsWeights(double HTXS_pTH=-1, int HTXS_Njets30=-1,int HTXS_cat=-1);
+ HiggsWeights getHiggsWeights(int HTXS_Njets30=-1, double HTXS_pTH=-99.0, int HTXS_cat=-1);
// returns hardcoded list of weight names matching expecation
// of ATLAS-default ggF NNLOPS
diff --git a/util/HiggsWeightTest.cxx b/util/HiggsWeightTest.cxx
index 2d9cb71f343c2e8bef43a1665f7096afded19ecf..ac2a71f4aacdcc7ef4f41c5e746c749b3f0d90d4 100644
--- a/util/HiggsWeightTest.cxx
+++ b/util/HiggsWeightTest.cxx
@@ -75,7 +75,7 @@ int main( int argc, char* argv[] ) {
// in reality we should pass the Higgs pT and Njets30 from HTXS
// but standard input files don't have these, so let's randomly sample them..
int HTXS_Njets30 = gRandom->Poisson(0.9);
- double HTXS_pTH = std::abs(gRandom->Gaus(0.0,40.0));
+ double HTXS_pTH = std::abs(gRandom->Gaus(0.0,50.0))*1000; // convert to MeV
if ( entry == 0 ) {
::Info(APP_NAME,"There are %lu weights in EventInfo and %lu in TruthEvents",