5 #ifndef TRUTHRIVETTOOLS_HIGGSTEMPLATECROSSSECTIONS_H
6 #define TRUTHRIVETTOOLS_HIGGSTEMPLATECROSSSECTIONS_H
9 #include "Rivet/Analysis.hh"
10 #include "Rivet/Particle.hh"
11 #include "Rivet/Projections/FastJets.hh"
18 #include "TLorentzVector.h"
37 :
Analysis(
"HiggsTemplateCrossSections"),
49 if ( ptcl.genParticle()->end_vertex() ) {
50 if ( !
hasChild(ptcl.genParticle(),ptcl.pid()) )
return ptcl;
58 auto prodVtx =
p.genParticle()->production_vertex();
59 if (prodVtx ==
nullptr)
return false;
62 for (
const auto &
part:ptcls )
63 if ( ancestor==
part.genParticle() )
return true;
77 if (child.pid()==pdgID)
return true;
84 if (
parent->pdg_id()==pdgID)
return true;
99 #
if RIVET_VERSION_CODE >= 40000
100 PID::isChargedLepton(child.pid())
111 std::string
msg=
"",
int NmaxWarnings=20)
const {
117 static std::atomic<int> Nwarnings = 0;
118 if ( !
msg.empty() && ++Nwarnings < NmaxWarnings )
129 HiggsClassification cat;
130 cat.prodMode = prodMode;
142 "Unkown Higgs production mechanism. Cannot classify event."
143 " Classification for all events will most likely fail.");
158 if ( ptcl->end_vertex() && !
hasChild(ptcl,PID::HIGGS) ) {
159 cat.higgs =
Particle(ptcl); ++Nhiggs;
163 if ( HSvtx==
nullptr && ptcl->production_vertex() && !
hasParent(ptcl,PID::HIGGS) )
164 HSvtx = ptcl->production_vertex();
170 "Current event has "+
std::to_string(Nhiggs)+
" Higgs bosons. There must be only one.");
171 if (cat.higgs.children().size()<2)
173 "Could not identify Higgs boson decay products.");
175 if (HSvtx ==
nullptr)
184 bool is_uncatdV =
false;
186 FourMomentum uncatV_p4(0,0,0,0);
187 FourVector uncatV_v4(0,0,0,0);
189 if (
isVH(prodMode) ) {
199 uncatV_p4 +=
Particle(ptcl).momentum();
200 uncatV_v4 +=
Particle(ptcl).origin();
203 is_uncatdV =
true; cat.V =
Particle(24,uncatV_p4,uncatV_v4);
209 if (
isVH(prodMode) && !cat.V.genParticle()->end_vertex() )
212 if (
isVH(prodMode) && cat.V.children().size()<2 )
215 if ( ( prodMode==
HTXS::WH && (nZs>0||nWs!=1) ) ||
218 std::to_string(nZs)+
" Z-bosons. Inconsitent with VH expectation.");
228 if (
top.genParticle()->end_vertex() )
229 for (
const auto &child:
top.children())
235 if ( (prodMode==
HTXS::TTH && Ws.size()<2) || (prodMode==
HTXS::TH && Ws.size()<1 ) )
250 }
else leptonicVs = std::move(uncatV_decays);
251 for (
const auto &
W:Ws )
if (
W.genParticle()->end_vertex() && !
quarkDecay(
W) ) leptonicVs +=
W;
256 FourMomentum
sum(0,0,0,0), vSum(0,0,0,0), hSum(0,0,0,0);
270 cat.p4decay_higgs = hSum;
271 cat.p4decay_V = vSum;
274 FastJets
jets(fps_temp,
275 #
if RIVET_VERSION_CODE >= 40000
302 cat.isZ2vvDecay =
false;
305 cat.stage1_cat_pTjet25GeV =
getStage1Category(prodMode,cat.higgs,cat.jets25,cat.V);
306 cat.stage1_cat_pTjet30GeV =
getStage1Category(prodMode,cat.higgs,cat.jets30,cat.V);
325 throw std::invalid_argument(
"Input value is out of bin range or bins vector is empty.");
328 for (
size_t i = 1;
i <
bins.size(); ++
i) {
330 return static_cast<int>(
i - 1);
334 return static_cast<int>(
bins.size() - 1);
341 if (
jets.size()<2)
return 0;
342 const FourMomentum &
j1=
jets[0].momentum(), &
j2=
jets[1].momentum();
343 bool VBFtopo = (
j1+
j2).
mass() > 400.0 && std::abs(
j1.rapidity()-
j2.rapidity()) > 2.8;
344 return VBFtopo ? (
j1+
j2+higgs.momentum()).
pt()<25 ? 2 : 1 : 0;
351 if (
jets.size()<2)
return 0;
352 const FourMomentum &
j1=
jets[0].momentum(), &
j2=
jets[1].momentum();
354 if(mjj>350 && mjj<=700)
return (
j1+
j2+higgs.momentum()).
pt()<25 ? 1 : 2;
355 else if(mjj>700)
return (
j1+
j2+higgs.momentum()).
pt()<25 ? 3 : 4;
365 if (
jets.size()<2)
return 0;
366 const FourMomentum &
j1=
jets[0].momentum(), &
j2=
jets[1].momentum();
368 if(mjj>350 && mjj<=700)
return (
j1+
j2+higgs.momentum()).
pt()<25 ? 1 : 2;
369 else if(mjj>700 && mjj<=1000)
return (
j1+
j2+higgs.momentum()).
pt()<25 ? 3 : 4;
370 else if(mjj>1000 && mjj<=1500)
return (
j1+
j2+higgs.momentum()).
pt()<25 ? 5 : 6;
371 else if(mjj>1500)
return (
j1+
j2+higgs.momentum()).
pt()<25 ? 7 : 8;
383 int ctrlHiggs = std::abs(higgs.rapidity())<2.5;
391 return Category(prodMode*10 + ctrlHiggs);
400 int Njets=
jets.size(), ctrlHiggs = std::abs(higgs.rapidity())<2.5, fwdHiggs = !ctrlHiggs;
401 double pTj1 =
jets.size() ?
jets[0].momentum().pt() : 0;
422 if (std::abs(higgs.rapidity())>2.5)
return QQ2HQQ_FWDH;
466 int Njets=
jets.size(), ctrlHiggs = std::abs(higgs.rapidity())<2.5, fwdHiggs = !ctrlHiggs;
486 if (std::abs(higgs.rapidity())>2.5)
return QQ2HQQ_FWDH;
487 int Njets=
jets.size();
495 else if ( mjj > 350 ) {
543 int Njets=
jets.size(), ctrlHiggs = std::abs(higgs.rapidity())<2.5, fwdHiggs = !ctrlHiggs;
550 if ( higgs.pt()>200 ){
552 double pTHj = (
jets[0].momentum()+higgs.momentum()).
pt();
562 double pTHjj = (
jets[0].momentum()+
jets[1].momentum()+higgs.momentum()).
pt();
573 if (std::abs(higgs.rapidity())>2.5)
return QQ2HQQ_FWDH;
574 int Njets=
jets.size();
579 double pTHjj = (
jets[0].momentum()+
jets[1].momentum()+higgs.momentum()).
pt();
595 int Njets=
jets.size();
603 int Njets=
jets.size();
611 int Njets=
jets.size();
639 printf(
"==============================================================\n");
640 printf(
"======== HiggsTemplateCrossSections Initialization =========\n");
641 printf(
"==============================================================\n");
645 char *pm_env =
getenv(
"HIGGSPRODMODE");
646 string pm(pm_env==
nullptr?
"":pm_env);
657 MSG_WARNING(
"No HIGGSPRODMODE shell variable found. Needed when running Rivet stand-alone.");
668 printf(
"==============================================================\n");
670 printf(
"======== Sucessful Initialization =========\n");
671 printf(
"==============================================================\n");
684 int F=cat.stage0_cat%10,
P=cat.stage1_cat_pTjet30GeV/100;
691 static const vector<int> offset1_2({0,1,18,29,35,41,47,53,55,57});
692 int off1_2 = offset1_2[
P];
694 static const vector<int> offset1_2_Fine({0,1,29,54,70,86,102,109,111,113});
695 int off1_2_Fine = offset1_2_Fine[
P];
716 if (cat.jets30.size()>=2) {
717 const FourMomentum &
j1 = cat.jets30[0].momentum(), &
j2 = cat.jets30[1].momentum();
725 MSG_INFO (
" ====================================================== ");
726 MSG_INFO (
" Higgs Template X-Sec Categorization Tool ");
728 MSG_INFO (
" ====================================================== ");
729 bool allSuccess = (numEvents()==m_errorCount[
HTXS::SUCCESS]);
730 if ( allSuccess )
MSG_INFO (
" >>>> All "<< m_errorCount[
HTXS::SUCCESS] <<
" events successfully categorized!");
733 MSG_INFO (
" >>>> --> the following errors occured:");
741 MSG_INFO (
" ====================================================== ");
742 MSG_INFO (
" ====================================================== ");
801 #ifdef RIVET_ANALYSIS_PATH