diff --git a/Root/HiggsWeightTool.cxx b/Root/HiggsWeightTool.cxx
index 2541a6e3661aba2cc5f575e646ecd95e055fb185..7e0f24e2ad9b896a9d2677ced20809a14e33bbc4 100644
--- a/Root/HiggsWeightTool.cxx
+++ b/Root/HiggsWeightTool.cxx
@@ -202,8 +202,10 @@ namespace xAOD {
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_b,hw.qcd_wg1_qm_t);
+ updateWeights(hw.nominal,hw.qcd_wg1_vbf2j,hw.qcd_wg1_vbf3j);
updateWeights(hw.nominal,hw.qcd_nnlops_nnlo,hw.qcd_nnlops_pow);
updateWeights(hw.nominal,hw.qcd_stxs);
+ updateWeights(hw.nominal,hw.qcd_jve);
}
void HiggsWeightTool::updateWeight(const double &w_nom, double &w) {
@@ -342,7 +344,8 @@ namespace xAOD {
// 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});
- static std::vector<double> sig({30.38,12.69,5.45}); // NNLOPS test sample
+ //static std::vector<double> sig({30.117,12.928,5.475}); // NNLOPS 2M
+ static std::vector<double> sig({30.117,12.928,4.845}); // NNLOPS subtracting VBF
// BLPTW absolute uncertainties in pb
static std::vector<double> yieldUnc({ 1.12, 0.66, 0.42});
@@ -376,10 +379,47 @@ namespace xAOD {
hw.qcd_wg1_qm_b = w_nom * qmSF;
hw.qcd_wg1_qm_t = w_nom * linInter(pTH,160,1.0,500,1.37);
+ hw.qcd_wg1_vbf2j = w_nom; hw.qcd_wg1_vbf3j = w_nom;
+ if (STXS_Stage1==101) { // GG2H_VBFTOPO_JET3VETO, tot unc 38%
+ hw.qcd_wg1_mu = hw.qcd_wg1_res = hw.qcd_wg1_mig01 = hw.qcd_wg1_mig12 = w_nom;
+ hw.qcd_wg1_vbf2j = w_nom * 1.20;
+ hw.qcd_wg1_vbf3j = w_nom * (1.0-0.32);
+ } if (STXS_Stage1==102) { // GG2H_VBFTOPO_JET3, tot unc 30.4%
+ hw.qcd_wg1_mu = hw.qcd_wg1_res = hw.qcd_wg1_mig01 = hw.qcd_wg1_mig12 = w_nom;
+ hw.qcd_wg1_vbf2j = w_nom * 1.20;
+ hw.qcd_wg1_vbf3j = w_nom * 1.235;
+ }
// Quark-mass uncdertainty - TODO
//if (pTH>500) qm=1.5; else if (pTH>160) qm=1.0+0.5*(pTH-160)/340;
//hw.qcd_wg1_qm = w_nom;
+ /********
+ * JVE as a cross check
+ */
+ // Taking eps0 and eps1 from Powheg NNLOPS
+ // and setting inclusive uncertainty to 3.9% (YR4 for N3LO)
+ // setting uncertatinies such that similar to the previous numbers
+ static double sig_ge1 = sig[1]+sig[2], sig_tot=sig[0]+sig_ge1,
+ Dsig_tot=sig_tot*0.039, D01=sig_tot*0.0242, D12=sig[1]*0.105;
+ //Dsig_tot=2.48, D01=1.25, D12=0.88;
+ hw.qcd_jve.push_back(w_nom*(1.0+Dsig_tot/sig_tot)); // incl
+ // eps0
+ if (Njets==0) hw.qcd_jve.push_back(w_nom*(1.0-D01/sig[0]));
+ else hw.qcd_jve.push_back(w_nom*(1.0+D01/sig_ge1));
+ // eps1
+ if (Njets==0) hw.qcd_jve.push_back(w_nom);
+ else if (Njets==1) hw.qcd_jve.push_back(w_nom*(1.0-D12/sig[1]));
+ else hw.qcd_jve.push_back(w_nom*(1.0+D12/sig[2]));
+
+ if (STXS_Stage1==101||STXS_Stage1==102)
+ hw.qcd_jve[0]=hw.qcd_jve[1]=hw.qcd_jve[2]=w_nom;
+ hw.qcd_jve.push_back(hw.qcd_wg1_vbf2j);
+ hw.qcd_jve.push_back(hw.qcd_wg1_vbf3j);
+ hw.qcd_jve.push_back(hw.qcd_wg1_pTH);
+ hw.qcd_jve.push_back(hw.qcd_wg1_qm_t);
+
+
+
/*********
* Tackmann STXS uncertainty scheme
*/
@@ -389,6 +429,10 @@ namespace xAOD {
hw.qcd_stxs.push_back(hw.qcd_wg1_res);
hw.qcd_stxs.push_back(hw.qcd_wg1_mig01);
hw.qcd_stxs.push_back(hw.qcd_wg1_mig12);
+
+ // Add in the VBF uncertainties here
+ hw.qcd_stxs.push_back(hw.qcd_wg1_vbf2j);
+ hw.qcd_stxs.push_back(hw.qcd_wg1_vbf3j);
// This is followed by Dsig60, Dsig120 and Dsig200
// These are extracted from Powheg NNLOPS scale variations (using this tool!)
diff --git a/TruthWeightTools/HiggsWeightTool.h b/TruthWeightTools/HiggsWeightTool.h
index 0638f205f7aa48fc654685feeadd7ea3bac9facd..2d7acfb3c1194b6139dd49a562a4f45973813870 100644
--- a/TruthWeightTools/HiggsWeightTool.h
+++ b/TruthWeightTools/HiggsWeightTool.h
@@ -35,11 +35,14 @@ 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_b, qcd_wg1_qm_t;
+ double qcd_wg1_pTH, qcd_wg1_qm_b, qcd_wg1_qm_t, qcd_wg1_vbf2j, qcd_wg1_vbf3j;
/// Tackmann proposed QCD uncertainty scheme, TODO
std::vector<double> qcd_stxs;
+ /// Tackmann proposed QCD uncertainty scheme, TODO
+ std::vector<double> qcd_jve;
+
/// Powheg NNLOPS possible scheme TODO
double qcd_nnlops_nnlo, qcd_nnlops_pow;
@@ -51,7 +54,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}; }
+ std::vector<double> qcd_wg1() { return {qcd_wg1_mu,qcd_wg1_res,qcd_wg1_mig01,qcd_wg1_mig12,
+ qcd_wg1_vbf2j,qcd_wg1_vbf3j,qcd_wg1_pTH,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); }
diff --git a/util/HiggsWeightTest.cxx b/util/HiggsWeightTest.cxx
index 6523ef683614e65fa822a492e74a2e7f848d5b10..6ccc757233fef6b8c0718c2dd41f43163c9d4ce0 100644
--- a/util/HiggsWeightTest.cxx
+++ b/util/HiggsWeightTest.cxx
@@ -53,7 +53,7 @@ int main( int argc, char* argv[] ) {
else if (arg=="--output") ofn=argv[++i];
else if (arg=="--weightCutOff") weightMax=atof(argv[++i]);
else if (arg.Contains(".root")) files.push_back(arg);
- else std::runtime_error(TString("Cannot intepret argument: "+arg).Data());
+ else throw std::runtime_error(TString("Cannot intepret argument: "+arg).Data());
}
if (files.size()==0) {
@@ -86,6 +86,14 @@ int main( int argc, char* argv[] ) {
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(7,"pTH_wg1qcd",Nbins,min,max,ptTit); // WG1 propsed scheme for ggF
+ HistV h_pTH_qcd_stxs = makeHistos(7,"pTH_qcd_stxs",Nbins,min,max,ptTit); // WG1 propsed scheme for ggF
+
+ TH1F *h_pTHGE1J = new TH1F("pTH_ge1J",ptTit,Nbins,min,max);
+ HistV h_pTH_qcdGE1J = makeHistos(8,"pTH_ge1J_qcd",Nbins,min,max,ptTit); // Default Powheg QCD variations (NLO)
+ HistV h_pTH_nnlops_qcdGE1J = makeHistos(26,"pTH_ge1J_nnlops_qcd",Nbins,min,max,ptTit); // NNLOPS internal QCD vars
+ HistV h_pTH_nnlops_qcd2GE1J = makeHistos(2,"pTH_ge1J_nnlo_qcd",Nbins,min,max,ptTit); // NNLO and Powheg vars for NNLOPS
+ HistV h_pTH_wg1qcdGE1J = makeHistos(7,"pTH_ge1J_wg1qcd",Nbins,min,max,ptTit); // WG1 propsed scheme for ggF
+ HistV h_pTH_qcd_stxsGE1J = makeHistos(7,"pTH_ge1J_qcd_stxs",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);
@@ -157,10 +165,13 @@ int main( int argc, char* argv[] ) {
event.getEntry( entry );
const xAOD::EventInfo *evtInfo;
RETURN_CHECK( APP_NAME, event.retrieve( evtInfo, "EventInfo" ) );
- std::vector<float> weights = evtInfo->mcEventWeights();
-
- const xAOD::TruthEventContainer *tevts;
- RETURN_CHECK( APP_NAME, event.retrieve( tevts, "TruthEvents" ) );
+ const xAOD::EventInfo *hgamtEvtInfo;
+ RETURN_CHECK( APP_NAME, event.retrieve( hgamtEvtInfo, "HGamTruthEventInfo" ) );
+ const xAOD::EventInfo *hgamEvtInfo;
+ RETURN_CHECK( APP_NAME, event.retrieve( hgamEvtInfo, "HGamEventInfo" ) );
+
+ bool dal = hgamEvtInfo->auxdata<char>("isDalitz");
+ if (dal) continue;
TLorentzVector h;
h.SetPtEtaPhiM(gRandom->Gaus(0,40),0,0,125);
@@ -171,12 +182,23 @@ int main( int argc, char* argv[] ) {
HTXS_Njets30 = evtInfo->auxdata<int>("HTXS_Njets_pTjet30");
HTXS_Stage1 = evtInfo->auxdata<int>("HTXS_Stage1_Category_pTjet30");
HTXS_index = evtInfo->auxdata<int>("HTXS_Stage1_FineIndex_pTjet30");
- HTXS_pTH = evtInfo->auxdata<float>("HTXS_Higgs_pt");
- double HTXS_etaH = evtInfo->auxdata<float>("HTXS_Higgs_eta");
- double HTXS_phiH = evtInfo->auxdata<float>("HTXS_Higgs_phi");
- double HTXS_mH = evtInfo->auxdata<float>("HTXS_Higgs_m");
+ //HTXS_pTH = evtInfo->auxdata<float>("HTXS_Higgs_pt");
+ HTXS_pTH = hgamtEvtInfo->auxdata<float>("pT_yy");
+ double HTXS_yH = hgamtEvtInfo->auxdata<float>("yAbs_yy");
+ double m = hgamtEvtInfo->auxdata<float>("m_yy");
+ double mT = sqrt(m*m+HTXS_pTH*HTXS_pTH);
+ double E = mT*cosh(HTXS_yH);
+ double pz = mT*sinh(HTXS_yH);
+ if (E<m) fatal("BAD");
+ if (HTXS_pTH<0) fatal("BAD, pT=0");
+ double p = sqrt(E*E-m*m);
+ double eta = 0.5*log((p+pz)/(p-pz));
+ //printf("pT = %.2f, y = %.3f, m=%.3f, eta=%.3f\n",HTXS_pTH,HTXS_yH,m,eta);
+ //double HTXS_etaH = evtInfo->auxdata<float>("HTXS_Higgs_eta");
+ //double HTXS_phiH = evtInfo->auxdata<float>("HTXS_Higgs_phi");
+ //double HTXS_mH = evtInfo->auxdata<float>("HTXS_Higgs_m");
- h.SetPtEtaPhiM(HTXS_pTH,HTXS_etaH,HTXS_phiH,HTXS_mH);
+ h.SetPtEtaPhiM(HTXS_pTH,eta,0,m);
h*=1e-3; // convert to GeV
} else HTXS_pTH=h.Pt()*1000;
@@ -196,6 +218,7 @@ int main( int argc, char* argv[] ) {
// 1. Nominal histograms
h_pTH -> Fill(h.Pt(),hw.nominal);
+ if (HTXS_Njets30>=1) h_pTHGE1J -> Fill(h.Pt(),hw.nominal);
h_Njets-> Fill(HTXS_Njets30,hw.nominal);
h_yH -> Fill(h.Rapidity(),hw.nominal);
h_STXS -> Fill(HTXS_index,hw.nominal);
@@ -234,6 +257,7 @@ int main( int argc, char* argv[] ) {
// without considering the NNLO correction. I.e. don't use for NNLOPS.
// Don't treat all these as NPs!! Take envelope. Or perhaps carefully chose 1 or 2 of them as NPs.
fillHistos(h_pTH_qcd,h.Pt(),hw.qcd);
+ if (HTXS_Njets30>=1) fillHistos(h_pTH_qcdGE1J,h.Pt(),hw.qcd);
fillHistos(h_yH_qcd,h.Rapidity(),hw.qcd);
fillHistos(h_Njets_qcd,HTXS_Njets30,hw.qcd);
fillHistos(h_STXS_qcd,HTXS_index,hw.qcd);
@@ -244,6 +268,7 @@ int main( int argc, char* argv[] ) {
// This should give an NNLO accurate normalization uncertainty (8-11%)
// Note: The NNLOPS paper takes the envelope of all these variations
fillHistos(h_pTH_nnlops_qcd,h.Pt(),hw.qcd_nnlops);
+ if (HTXS_Njets30>=1) fillHistos(h_pTH_nnlops_qcdGE1J,h.Pt(),hw.qcd_nnlops);
fillHistos(h_yH_nnlops_qcd,h.Rapidity(),hw.qcd_nnlops);
fillHistos(h_Njets_nnlops_qcd,HTXS_Njets30,hw.qcd_nnlops);
fillHistos(h_STXS_nnlops_qcd,HTXS_index,hw.qcd_nnlops);
@@ -253,6 +278,7 @@ int main( int argc, char* argv[] ) {
// (nnloNom-PowDnDn) as an uncorrelation source (similar to >=1 in ST: affects high pT)
NumV nnlops_2np_qcd={hw.qcd_nnlops_nnlo,hw.qcd_nnlops_pow};
fillHistos(h_pTH_nnlops_qcd2,h.Pt(),nnlops_2np_qcd);
+ if (HTXS_Njets30>=1) fillHistos(h_pTH_nnlops_qcd2GE1J,h.Pt(),nnlops_2np_qcd);
fillHistos(h_yH_nnlops_qcd2,h.Rapidity(),nnlops_2np_qcd);
fillHistos(h_Njets_nnlops_qcd2,HTXS_Njets30,nnlops_2np_qcd);
fillHistos(h_STXS_nnlops_qcd2,HTXS_index,nnlops_2np_qcd);
@@ -261,10 +287,14 @@ int main( int argc, char* argv[] ) {
//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);
+ if (HTXS_Njets30>=1) fillHistos(h_pTH_wg1qcdGE1J,h.Pt(),wg1_qcd);
fillHistos(h_yH_wg1qcd,h.Rapidity(),wg1_qcd);
fillHistos(h_Njets_wg1qcd,HTXS_Njets30,wg1_qcd);
fillHistos(h_STXS_wg1qcd,HTXS_index,wg1_qcd);
+ fillHistos(h_pTH_qcd_stxs,h.Pt(),hw.qcd_stxs);
+ if (HTXS_Njets30>=1) fillHistos(h_pTH_qcd_stxsGE1J,h.Pt(),hw.qcd_stxs);
+
// 4. Other weights
fillHistos(h_yH_other,h.Rapidity(),
{hw.weight0,hw.nnpdf30_nnlo,hw.pdf4lhc_nnlo,hw.nnpdf30_nlo,hw.pdf4lhc_nlo,
@@ -275,12 +305,6 @@ int main( int argc, char* argv[] ) {
// Print stuff to the screen for the first event in each file
if ( entry == 0 ) {
- ::Info(APP_NAME,"There are %lu weights in EventInfo and %lu in TruthEvents",
- weights.size(),tevts->at(0)->weights().size());
- std::vector<float> ws = tevts->at(0)->weights();
- for (size_t i=0;i<10;++i)
- ::Info(APP_NAME,"Weight %lu %.3f and %.3f. %lu weights and %lu names",
- i,weights[i],ws[i],ws.size(),higgsMCtool->getWeightNames().size());
hw.print();
}
diff --git a/util/STXSacc.cxx b/util/STXSacc.cxx
index e84a384b4c821288eca6f5bd46fc4046759b1008..33703a70c9672b26637b9698d7d0e2524cd83815 100644
--- a/util/STXSacc.cxx
+++ b/util/STXSacc.cxx
@@ -58,7 +58,7 @@ int main( int argc, char* argv[] ) {
xAOD::HiggsWeightTool *higgsMCtool = new xAOD::HiggsWeightTool( "HiggsWeightTool" );
higgsMCtool->setProperty( "OutputLevel", MSG::DEBUG ).ignore();
higgsMCtool->setProperty("ForceNNLOPS",true).ignore();
- higgsMCtool->setProperty("WeightCutOff",100.0).ignore();
+ //higgsMCtool->setProperty("WeightCutOff",200.0).ignore();
higgsMCtool->initialize().ignore();
// Open the input file: