diff --git a/Root/HiggsWeightTool.cxx b/Root/HiggsWeightTool.cxx
index ffedbdeae3de3d84eb368b49b9f10b084d007668..2541a6e3661aba2cc5f575e646ecd95e055fb185 100644
--- a/Root/HiggsWeightTool.cxx
+++ b/Root/HiggsWeightTool.cxx
@@ -201,7 +201,7 @@ namespace xAOD {
updateWeights(hw.nominal,hw.ct14nlo,hw.ct14nlo_0118);
updateWeights(hw.nominal,hw.qcd_wg1_mu,hw.qcd_wg1_res);
updateWeights(hw.nominal,hw.qcd_wg1_mig01,hw.qcd_wg1_mig12);
- updateWeights(hw.nominal,hw.qcd_wg1_pTH,hw.qcd_wg1_qm);
+ updateWeights(hw.nominal,hw.qcd_wg1_pTH,hw.qcd_wg1_qm_b,hw.qcd_wg1_qm_t);
updateWeights(hw.nominal,hw.qcd_nnlops_nnlo,hw.qcd_nnlops_pow);
updateWeights(hw.nominal,hw.qcd_stxs);
}
@@ -362,13 +362,23 @@ namespace xAOD {
// 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*w_nom;
+ double pTH_uncSF = 1.0;
+ if (Njets==1) pTH_uncSF = linInter(pTH,0,y1_1,x2_1,y2_1);
+ else if (Njets>=2) pTH_uncSF = linInter(pTH,0,y1_ge2,x2_ge2,y2_ge2);
+ //pTH>x2_ge2?y2_ge2:y1_ge2+(y2_ge2-y1_ge2)*pTH/x2_ge2;
+ //else if (Njets==1) pTH_uncSF = pTH>x2_1?y2_1:y1_1+(y2_1-y1_1)*pTH/x2_1;
+ hw.qcd_wg1_pTH = pTH_uncSF*w_nom;
+
+ double qmSF=1.0;
+ if (Njets==0) qmSF = linInter(pTH,6,0.92,24,1.05);
+ else if (Njets==1) qmSF = linInter(pTH,35,0.92,62.5,1.05);
+ else if (Njets>=2) qmSF = linInter(pTH,40,0.92,120,1.05);
+ hw.qcd_wg1_qm_b = w_nom * qmSF;
+ hw.qcd_wg1_qm_t = w_nom * linInter(pTH,160,1.0,500,1.37);
// Quark-mass uncdertainty - TODO
- hw.qcd_wg1_qm = w_nom;
+ //if (pTH>500) qm=1.5; else if (pTH>160) qm=1.0+0.5*(pTH-160)/340;
+ //hw.qcd_wg1_qm = w_nom;
/*********
* Tackmann STXS uncertainty scheme
@@ -383,12 +393,13 @@ namespace xAOD {
// This is followed by Dsig60, Dsig120 and Dsig200
// These are extracted from Powheg NNLOPS scale variations (using this tool!)
// As the envelope of hw.qcd_nnlops, which gives (all x-sec below have Njets>=1):
- // sig60 = 9.687 +/- 1.566 pb
- // sig120 = 2.534 +/- 0.526 pb
- // sig200 = 0.562 +/- 0.115 pb
- static double sig0_60=8.458, sig60_200=9.125, sig120_200=1.972,
- sig0_120=sig0_60+sig60_200-sig120_200, sig200_plus=0.562;
- static double Dsig60_200=1.457, Dsig120_200=0.411, Dsig200_plus=0.115;
+ // sig(60,200) = 9.095 +/- 1.445 pb, BLPTW 10.9%
+ // sig(120,200) = 1.961 +/- 0.401 pb, BLPTW 13.1%
+ // sig(200,inf) = 0.582 +/- 0.121 pb, BLPTW 15.1%
+ static double sig0_60=8.719, sig60_200=9.095, sig120_200=1.961,
+ sig0_120=sig0_60+sig60_200-sig120_200, sig200_plus=0.582; // 0.121 (-) 0.151*0.582
+ // static double Dsig60_200=1.457, Dsig120_200=0.411, Dsig200_plus=0.115; // with 40k events, and no BLPTW subraction
+ static double Dsig60_200=1.055, Dsig120_200=0.206, Dsig200_plus=0.0832; // with 2M evts, and subtraction
double dsig60=0, dsig120=0, dsig200=0;
if (Njets>=1) {
if (pTH<60) dsig60=-Dsig60_200/sig0_60; // -17.2%
@@ -397,7 +408,7 @@ namespace xAOD {
if (pTH<120) dsig120 = -Dsig120_200/sig0_120; // -2.6%
else if (pTH<200) dsig120 = Dsig120_200/sig120_200; // +20.8%
- if (pTH>200) dsig200=Dsig200_plus/sig200_plus; // +20.5%
+ if (pTH>200) dsig200=Dsig200_plus/sig200_plus; // +14.3%
}
hw.qcd_stxs.push_back((1.0+dsig60)*w_nom);
hw.qcd_stxs.push_back((1.0+dsig120)*w_nom);
diff --git a/TruthWeightTools/HiggsWeightTool.h b/TruthWeightTools/HiggsWeightTool.h
index 0cc9cf53d5b072b0f7fe27a239feec8764031e4e..0638f205f7aa48fc654685feeadd7ea3bac9facd 100644
--- a/TruthWeightTools/HiggsWeightTool.h
+++ b/TruthWeightTools/HiggsWeightTool.h
@@ -35,7 +35,7 @@ namespace xAOD {
/// WG1 proposed QCD uncertainty scheme
double qcd_wg1_mu, qcd_wg1_res, qcd_wg1_mig01, qcd_wg1_mig12;
- double qcd_wg1_pTH, qcd_wg1_qm;
+ double qcd_wg1_pTH, qcd_wg1_qm_b, qcd_wg1_qm_t;
/// Tackmann proposed QCD uncertainty scheme, TODO
std::vector<double> qcd_stxs;
@@ -51,6 +51,8 @@ namespace xAOD {
double pTH;
int Njets30, STXS;
+ std::vector<double> qcd_wg1() { return {qcd_wg1_mu,qcd_wg1_res,qcd_wg1_mig01,qcd_wg1_mig12,qcd_wg1_pTH,qcd_wg1_qm_b,qcd_wg1_qm_t}; }
+
/// methods to print weights to the screen
char *uncStr(double var, double nom) { return var==0?Form(" N/A"):Form("%s%.1f%%",var>=nom?"+":"",(var-nom)/nom*100); }
void print() {
@@ -84,8 +86,8 @@ namespace xAOD {
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));
+ printf(" pTH: %s, quark-mass, b: %s, t: %s\n",
+ uncStr(qcd_wg1_pTH,n),uncStr(qcd_wg1_qm_b,n),uncStr(qcd_wg1_qm_b,n));
}
if (qcd.size()) {
@@ -151,6 +153,13 @@ namespace xAOD {
private:
+ /// linear interpolation
+ double linInter(double x, double x1, double y1, double x2, double y2) {
+ if (x<x1) return y1;
+ if (x>x2) return y2;
+ return y1+(y2-y1)*(x-x1)/(x2-x1);
+ }
+
/// Access the HiggsWeights
HiggsWeights getHiggsWeightsInternal(int HTXS_Njets30=-1, double HTXS_pTH=-99.0, int HTXS_cat=-1);
diff --git a/util/HiggsWeightTest.cxx b/util/HiggsWeightTest.cxx
index 090f1125971948334a09521717425689d22c6300..6523ef683614e65fa822a492e74a2e7f848d5b10 100644
--- a/util/HiggsWeightTest.cxx
+++ b/util/HiggsWeightTest.cxx
@@ -85,7 +85,7 @@ int main( int argc, char* argv[] ) {
HistV h_pTH_qcd = makeHistos(8,"pTH_qcd",Nbins,min,max,ptTit); // Default Powheg QCD variations (NLO)
HistV h_pTH_nnlops_qcd = makeHistos(26,"pTH_nnlops_qcd",Nbins,min,max,ptTit); // NNLOPS internal QCD vars
HistV h_pTH_nnlops_qcd2 = makeHistos(2,"pTH_nnlo_qcd",Nbins,min,max,ptTit); // NNLO and Powheg vars for NNLOPS
- HistV h_pTH_wg1qcd = makeHistos(6,"pTH_wg1qcd",Nbins,min,max,ptTit); // WG1 propsed scheme for ggF
+ HistV h_pTH_wg1qcd = makeHistos(7,"pTH_wg1qcd",Nbins,min,max,ptTit); // WG1 propsed scheme for ggF
Nbins=10; min=-0.5; max=9.5; Str tit=";#it{N}_{jets}";
TH1F *h_Njets = new TH1F("Njets30",tit,Nbins,min,max);
@@ -96,7 +96,7 @@ int main( int argc, char* argv[] ) {
HistV h_Njets_qcd = makeHistos(8,"Njets30_qcd",Nbins,min,max,tit); // Default Powheg QCD variations (NLO)
HistV h_Njets_nnlops_qcd = makeHistos(26,"Njets30_nnlops_qcd",Nbins,min,max,tit); // NNLOPS internal QCD vars
HistV h_Njets_nnlops_qcd2 = makeHistos(2,"Njets30_nnlo_qcd",Nbins,min,max,tit); // NNLO and Powheg vars for NNLOPS
- HistV h_Njets_wg1qcd = makeHistos(6,"Njets30_wg1qcd",Nbins,min,max,tit); // WG1 propsed scheme for ggF
+ HistV h_Njets_wg1qcd = makeHistos(7,"Njets30_wg1qcd",Nbins,min,max,tit); // WG1 propsed scheme for ggF
Nbins=52; min=1; max=53; tit=";STXS fine index";
TH1F *h_STXS = new TH1F("STXS",tit,Nbins,min,max);
@@ -107,7 +107,7 @@ int main( int argc, char* argv[] ) {
HistV h_STXS_qcd = makeHistos(8,"STXS_qcd",Nbins,min,max,tit); // Default Powheg QCD variations (NLO)
HistV h_STXS_nnlops_qcd = makeHistos(26,"STXS_nnlops_qcd",Nbins,min,max,tit); // NNLOPS internal QCD vars
HistV h_STXS_nnlops_qcd2 = makeHistos(2,"STXS_nnlo_qcd",Nbins,min,max,tit); // NNLO and Powheg vars for NNLOPS
- HistV h_STXS_wg1qcd = makeHistos(6,"STXS_wg1qcd",Nbins,min,max,tit); // WG1 propsed scheme for ggF
+ HistV h_STXS_wg1qcd = makeHistos(7,"STXS_wg1qcd",Nbins,min,max,tit); // WG1 propsed scheme for ggF
Nbins=80; min=-4; max=4; tit=";#it{y_{H}}";
TH1F *h_yH = new TH1F("yH",tit,Nbins,min,max);
@@ -118,7 +118,7 @@ int main( int argc, char* argv[] ) {
HistV h_yH_qcd = makeHistos(8,"yH_qcd",Nbins,min,max,tit); // Default Powheg QCD variations (NLO)
HistV h_yH_nnlops_qcd = makeHistos(26,"yH_nnlops_qcd",Nbins,min,max,tit); // NNLOPS internal QCD vars
HistV h_yH_nnlops_qcd2 = makeHistos(2,"yH_nnlo_qcd",Nbins,min,max,tit); // NNLO and Powheg vars for NNLOPS
- HistV h_yH_wg1qcd = makeHistos(6,"yH_wg1qcd",Nbins,min,max,tit); // WG1 propsed scheme for ggF
+ HistV h_yH_wg1qcd = makeHistos(7,"yH_wg1qcd",Nbins,min,max,tit); // WG1 propsed scheme for ggF
HistV h_yH_other = makeHistos(8,"yH_other",Nbins,min,max,tit); // WG1 propsed scheme for ggF
// Initialise the application:
@@ -258,7 +258,8 @@ int main( int argc, char* argv[] ) {
fillHistos(h_STXS_nnlops_qcd2,HTXS_index,nnlops_2np_qcd);
// 3.d WG1 propsed uncertainty scheme. Recommended.
- NumV wg1_qcd={hw.qcd_wg1_mu,hw.qcd_wg1_res,hw.qcd_wg1_mig01,hw.qcd_wg1_mig12,hw.qcd_wg1_pTH,hw.qcd_wg1_qm};
+ //NumV wg1_qcd={hw.qcd_wg1_mu,hw.qcd_wg1_res,hw.qcd_wg1_mig01,hw.qcd_wg1_mig12,hw.qcd_wg1_pTH,hw.qcd_wg1_qm};
+ NumV wg1_qcd=hw.qcd_wg1();
fillHistos(h_pTH_wg1qcd,h.Pt(),wg1_qcd);
fillHistos(h_yH_wg1qcd,h.Rapidity(),wg1_qcd);
fillHistos(h_Njets_wg1qcd,HTXS_Njets30,wg1_qcd);