Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-07-03 04:17:53 

0001 #ifndef TRUTHRIVETTOOLS_HIGGSTEMPLATECROSSSECTIONS_CC
0002 #define TRUTHRIVETTOOLS_HIGGSTEMPLATECROSSSECTIONS_CC
0003 
0004 // -*- C++ -*-
0005 #include "Rivet/Analysis.hh"
0006 #include "Rivet/Particle.hh"
0007 #include "Rivet/Projections/FastJets.hh"
0008 
0009 // Definition of the StatusCode and Category enums
0010 //#include "HiggsTemplateCrossSections.h"
0011 #include "SimDataFormats/HTXS/interface/HiggsTemplateCrossSections.h"  //
0012 
0013 #include <atomic>
0014 
0015 namespace Rivet {
0016 
0017   /// @class HiggsTemplateCrossSections
0018   /// @brief  Rivet routine for classifying MC events according to the Higgs template cross section categories
0019   /// @author Jim Lacey (DESY) <james.lacey@cern.ch,jlacey@desy.de>
0020   /// @author Dag Gillberg (Carleton University) <dag.gillberg@cern.ch>
0021   class HiggsTemplateCrossSections : public Analysis {
0022   public:
0023     // Constructor
0024     HiggsTemplateCrossSections() : Analysis("HiggsTemplateCrossSections"), m_HiggsProdMode(HTXS::UNKNOWN) {}
0025 
0026   public:
0027     /// @name Utility methods
0028     /// Methods to identify the Higgs boson and
0029     /// associated vector boson and to build jets
0030     /// @{
0031 
0032     /// follow a "propagating" particle and return its last instance
0033     Particle getLastInstance(Particle ptcl) {
0034       if (ptcl.genParticle()->end_vertex()) {
0035         if (!hasChild(ptcl.genParticle(), ptcl.pid()))
0036           return ptcl;
0037         else
0038           return getLastInstance(ptcl.children()[0]);
0039       }
0040       return ptcl;
0041     }
0042 
0043     /// @brief Whether particle p originate from any of the ptcls
0044     bool originateFrom(const Particle &p, const Particles &ptcls) {
0045       const ConstGenVertexPtr prodVtx = p.genParticle()->production_vertex();
0046       if (prodVtx == nullptr)
0047         return false;
0048       // for each ancestor, check if it matches any of the input particles
0049       for (const auto &ancestor : HepMCUtils::particles(prodVtx, Relatives::ANCESTORS)) {
0050         for (const auto &part : ptcls)
0051           if (ancestor == part.genParticle())
0052             return true;
0053       }
0054       // if we get here, no ancetor matched any input particle
0055       return false;
0056     }
0057 
0058     /// @brief Whether particle p originates from p2
0059     bool originateFrom(const Particle &p, const Particle &p2) {
0060       Particles ptcls = {p2};
0061       return originateFrom(p, ptcls);
0062     }
0063 
0064     /// @brief Checks whether the input particle has a child with a given PDGID
0065     bool hasChild(const ConstGenParticlePtr &ptcl, int pdgID) {
0066       for (const auto &child : Particle(*ptcl).children()) {
0067         if (child.pid() == pdgID) {
0068           return true;
0069         }
0070       }
0071       return false;
0072     }
0073 
0074     /// @brief Checks whether the input particle has a parent with a given PDGID
0075     bool hasParent(const ConstGenParticlePtr &ptcl, int pdgID) {
0076       for (auto parent : HepMCUtils::particles(ptcl->production_vertex(), Relatives::PARENTS))
0077         if (parent->pdg_id() == pdgID)
0078           return true;
0079       return false;
0080     }
0081 
0082     /// @brief Return true is particle decays to quarks
0083     bool quarkDecay(const Particle &p) {
0084       for (const auto &child : p.children()) {
0085         if (PID::isQuark(child.pid())) {
0086           return true;
0087         }
0088       }
0089       return false;
0090     }
0091 
0092     /// @brief Return true if particle decays to charged leptons.
0093     bool ChLeptonDecay(const Particle &p) {
0094       for (const auto &child : p.children()) {
0095         if (PID::isChargedLepton(child.pid())) {
0096           return true;
0097         }
0098       }
0099       return false;
0100     }
0101 
0102     /// @brief Returns the classification object with the error code set.
0103     ///        Prints an warning message, and keeps track of number of errors
0104     HiggsClassification error(HiggsClassification &cat,
0105                               HTXS::ErrorCode err,
0106                               std::string msg = "",
0107                               int NmaxWarnings = 20) {
0108       // Set the error, and keep statistics
0109       cat.errorCode = err;
0110       ++m_errorCount[err];
0111 
0112       // Print warning message to the screen/log
0113       static std::atomic<int> Nwarnings{0};
0114       if (!msg.empty() && ++Nwarnings < NmaxWarnings)
0115         MSG_WARNING(msg);
0116 
0117       return cat;
0118     }
0119     /// @}
0120 
0121     /// @brief Main classificaion method.
0122     HiggsClassification classifyEvent(const Event &event, const HTXS::HiggsProdMode prodMode) {
0123       if (m_HiggsProdMode == HTXS::UNKNOWN)
0124         m_HiggsProdMode = prodMode;
0125 
0126       // the classification object
0127       HiggsClassification cat;
0128       cat.prodMode = prodMode;
0129       cat.errorCode = HTXS::UNDEFINED;
0130       cat.stage0_cat = HTXS::Stage0::UNKNOWN;
0131       cat.stage1_cat_pTjet25GeV = HTXS::Stage1::UNKNOWN;
0132       cat.stage1_cat_pTjet30GeV = HTXS::Stage1::UNKNOWN;
0133       cat.stage1_1_cat_pTjet25GeV = HTXS::Stage1_1::UNKNOWN;
0134       cat.stage1_1_cat_pTjet30GeV = HTXS::Stage1_1::UNKNOWN;
0135       cat.stage1_1_fine_cat_pTjet25GeV = HTXS::Stage1_1_Fine::UNKNOWN;
0136       cat.stage1_1_fine_cat_pTjet30GeV = HTXS::Stage1_1_Fine::UNKNOWN;
0137       cat.stage1_2_cat_pTjet25GeV = HTXS::Stage1_2::UNKNOWN;
0138       cat.stage1_2_cat_pTjet30GeV = HTXS::Stage1_2::UNKNOWN;
0139       cat.stage1_2_fine_cat_pTjet25GeV = HTXS::Stage1_2_Fine::UNKNOWN;
0140       cat.stage1_2_fine_cat_pTjet30GeV = HTXS::Stage1_2_Fine::UNKNOWN;
0141 
0142       if (prodMode == HTXS::UNKNOWN)
0143         return error(cat,
0144                      HTXS::PRODMODE_DEFINED,
0145                      "Unkown Higgs production mechanism. Cannot classify event."
0146                      " Classification for all events will most likely fail.");
0147 
0148       /*****
0149        * Step 1. 
0150        *  Idenfify the Higgs boson and the hard scatter vertex
0151        *  There should be only one of each.
0152        */
0153 
0154       ConstGenVertexPtr HSvtx = nullptr;
0155       int Nhiggs = 0;
0156       for (const ConstGenParticlePtr &ptcl : HepMCUtils::particles(event.genEvent())) {
0157         // a) Reject all non-Higgs particles
0158         if (!PID::isHiggs(ptcl->pdg_id()))
0159           continue;
0160         // b) select only the final Higgs boson copy, prior to decay
0161         if (ptcl->end_vertex() && !hasChild(ptcl, PID::HIGGS)) {
0162           cat.higgs = Particle(ptcl);
0163           ++Nhiggs;
0164         }
0165         // c) HepMC3 does not provide signal_proces_vertex anymore
0166         //    set hard-scatter vertex based on first Higgs boson
0167         if (HSvtx == nullptr && ptcl->production_vertex() && !hasParent(ptcl, PID::HIGGS))
0168           HSvtx = ptcl->production_vertex();
0169       }
0170 
0171       // Make sure things are in order so far
0172       if (Nhiggs != 1)
0173         return error(cat,
0174                      HTXS::HIGGS_IDENTIFICATION,
0175                      "Current event has " + std::to_string(Nhiggs) + " Higgs bosons. There must be only one.");
0176       if (cat.higgs.children().size() < 2)
0177         return error(cat, HTXS::HIGGS_DECAY_IDENTIFICATION, "Could not identify Higgs boson decay products.");
0178 
0179       if (HSvtx == nullptr)
0180         return error(cat, HTXS::HS_VTX_IDENTIFICATION, "Cannot find hard-scatter vertex of current event.");
0181 
0182       /*****
0183        * Step 2. 
0184        *   Identify associated vector bosons
0185        */
0186 
0187       // Find associated vector bosons
0188       bool is_uncatdV = false;
0189       Particles uncatV_decays;
0190       FourMomentum uncatV_p4(0, 0, 0, 0);
0191       FourVector uncatV_v4(0, 0, 0, 0);
0192       int nWs = 0, nZs = 0;
0193       if (isVH(prodMode)) {
0194         for (auto ptcl : HepMCUtils::particles(HSvtx, Relatives::CHILDREN)) {
0195           if (PID::isW(ptcl->pdg_id())) {
0196             ++nWs;
0197             cat.V = Particle(ptcl);
0198           }
0199           if (PID::isZ(ptcl->pdg_id())) {
0200             ++nZs;
0201             cat.V = Particle(ptcl);
0202           }
0203         }
0204         if (nWs + nZs > 0)
0205           cat.V = getLastInstance(cat.V);
0206         else {
0207           for (auto ptcl : HepMCUtils::particles(HSvtx, Relatives::CHILDREN)) {
0208             if (!PID::isHiggs(ptcl->pdg_id())) {
0209               uncatV_decays += Particle(ptcl);
0210               uncatV_p4 += Particle(ptcl).momentum();
0211               // uncatV_v4 += Particle(ptcl).origin();
0212             }
0213           }
0214           // is_uncatdV = true; cat.V = Particle(24,uncatV_p4,uncatV_v4);
0215           is_uncatdV = true;
0216           cat.V = Particle(24, uncatV_p4);
0217         }
0218       }
0219 
0220       if (!is_uncatdV) {
0221         if (isVH(prodMode) && !cat.V.genParticle()->end_vertex())
0222           return error(cat, HTXS::VH_DECAY_IDENTIFICATION, "Vector boson does not decay!");
0223 
0224         if (isVH(prodMode) && cat.V.children().size() < 2)
0225           return error(cat, HTXS::VH_DECAY_IDENTIFICATION, "Vector boson does not decay!");
0226 
0227         if ((prodMode == HTXS::WH && (nZs > 0 || nWs != 1)) ||
0228             ((prodMode == HTXS::QQ2ZH || prodMode == HTXS::GG2ZH) && (nZs != 1 || nWs > 0)))
0229           return error(cat,
0230                        HTXS::VH_IDENTIFICATION,
0231                        "Found " + std::to_string(nWs) + " W-bosons and " + std::to_string(nZs) +
0232                            " Z-bosons. Inconsitent with VH expectation.");
0233       }
0234 
0235       // Find and store the W-bosons from ttH->WbWbH
0236       Particles Ws;
0237       if (prodMode == HTXS::TTH || prodMode == HTXS::TH) {
0238         // loop over particles produced in hard-scatter vertex
0239         for (auto ptcl : HepMCUtils::particles(HSvtx, Relatives::CHILDREN)) {
0240           if (!PID::isTop(ptcl->pdg_id()))
0241             continue;
0242           Particle top = getLastInstance(Particle(ptcl));
0243           if (top.genParticle()->end_vertex())
0244             for (const auto &child : top.children())
0245               if (PID::isW(child.pid()))
0246                 Ws += getLastInstance(child);
0247         }
0248       }
0249 
0250       // Make sure result make sense
0251       if ((prodMode == HTXS::TTH && Ws.size() < 2) || (prodMode == HTXS::TH && Ws.empty()))
0252         return error(cat, HTXS::TOP_W_IDENTIFICATION, "Failed to identify W-boson(s) from t-decay!");
0253 
0254       /*****
0255        * Step 3.
0256        *   Build jets
0257        *   Make sure all stable particles are present
0258        */
0259 
0260       // Create a list of the vector bosons that decay leptonically
0261       // Either the vector boson produced in association with the Higgs boson,
0262       // or the ones produced from decays of top quarks produced with the Higgs
0263       Particles leptonicVs;
0264       if (!is_uncatdV) {
0265         if (isVH(prodMode) && !quarkDecay(cat.V))
0266           leptonicVs += cat.V;
0267       } else
0268         leptonicVs = uncatV_decays;
0269       for (const auto &W : Ws)
0270         if (W.genParticle()->end_vertex() && !quarkDecay(W))
0271           leptonicVs += W;
0272 
0273       // Obtain all stable, final-state particles
0274       const Particles FS = apply<FinalState>(event, "FS").particles();
0275       Particles hadrons;
0276       FourMomentum sum(0, 0, 0, 0), vSum(0, 0, 0, 0), hSum(0, 0, 0, 0);
0277       for (const Particle &p : FS) {
0278         // Add up the four momenta of all stable particles as a cross check
0279         sum += p.momentum();
0280         // ignore particles from the Higgs boson
0281         if (originateFrom(p, cat.higgs)) {
0282           hSum += p.momentum();
0283           continue;
0284         }
0285         // Cross-check the V decay products for VH
0286         if (isVH(prodMode) && !is_uncatdV && originateFrom(p, Ws))
0287           vSum += p.momentum();
0288         // ignore final state particles from leptonic V decays
0289         if (!leptonicVs.empty() && originateFrom(p, leptonicVs))
0290           continue;
0291         // All particles reaching here are considered hadrons and will be used to build jets
0292         hadrons += p;
0293       }
0294 
0295       cat.p4decay_higgs = hSum;
0296       cat.p4decay_V = vSum;
0297 
0298       FinalState fps_temp;
0299       FastJets jets(fps_temp, JetAlg::ANTIKT, 0.4);
0300       jets.calc(hadrons);
0301 
0302       cat.jets25 = jets.jetsByPt(Cuts::pT > 25.0);
0303       cat.jets30 = jets.jetsByPt(Cuts::pT > 30.0);
0304 
0305       // check that four mometum sum of all stable particles satisfies momentum consevation
0306       /*
0307       if ( sum.pt()>0.1 )
0308     return error(cat,HTXS::MOMENTUM_CONSERVATION,"Four vector sum does not amount to pT=0, m=E=sqrt(s), but pT="+
0309              std::to_string(sum.pt())+" GeV and m = "+std::to_string(sum.mass())+" GeV");
0310 */
0311       // check if V-boson was not included in the event record but decay particles were
0312       // EFT contact interaction: return UNKNOWN for category but set all event/particle kinematics
0313       if (is_uncatdV)
0314         return error(cat, HTXS::VH_IDENTIFICATION, "Failed to identify associated V-boson!");
0315 
0316       /*****
0317        * Step 4.
0318        *   Classify and save output
0319        */
0320 
0321       // Apply the categorization categorization
0322       cat.isZ2vvDecay = false;
0323       if ((prodMode == HTXS::GG2ZH || prodMode == HTXS::QQ2ZH) && !quarkDecay(cat.V) && !ChLeptonDecay(cat.V))
0324         cat.isZ2vvDecay = true;
0325       cat.stage0_cat = getStage0Category(prodMode, cat.higgs, cat.V);
0326       cat.stage1_cat_pTjet25GeV = getStage1Category(prodMode, cat.higgs, cat.jets25, cat.V);
0327       cat.stage1_cat_pTjet30GeV = getStage1Category(prodMode, cat.higgs, cat.jets30, cat.V);
0328       cat.stage1_1_cat_pTjet25GeV = getStage1_1_Category(prodMode, cat.higgs, cat.jets25, cat.V);
0329       cat.stage1_1_cat_pTjet30GeV = getStage1_1_Category(prodMode, cat.higgs, cat.jets30, cat.V);
0330       cat.stage1_1_fine_cat_pTjet25GeV = getStage1_1_Fine_Category(prodMode, cat.higgs, cat.jets25, cat.V);
0331       cat.stage1_1_fine_cat_pTjet30GeV = getStage1_1_Fine_Category(prodMode, cat.higgs, cat.jets30, cat.V);
0332       cat.stage1_2_cat_pTjet25GeV = getStage1_2_Category(prodMode, cat.higgs, cat.jets25, cat.V);
0333       cat.stage1_2_cat_pTjet30GeV = getStage1_2_Category(prodMode, cat.higgs, cat.jets30, cat.V);
0334       cat.stage1_2_fine_cat_pTjet25GeV = getStage1_2_Fine_Category(prodMode, cat.higgs, cat.jets25, cat.V);
0335       cat.stage1_2_fine_cat_pTjet30GeV = getStage1_2_Fine_Category(prodMode, cat.higgs, cat.jets30, cat.V);
0336 
0337       cat.errorCode = HTXS::SUCCESS;
0338       ++m_errorCount[HTXS::SUCCESS];
0339       ++m_sumevents;
0340 
0341       return cat;
0342     }
0343 
0344     /// @name Categorization methods
0345     /// Methods to assign the truth category based
0346     /// on the identified Higgs boson and associated
0347     /// vector bosons and/or reconstructed jets
0348     /// @{
0349 
0350     /// @brief Return bin index of x given the provided bin edges. 0=first bin, -1=underflow bin.
0351     int getBin(double x, std::vector<double> bins) {
0352       for (size_t i = 1; i < bins.size(); ++i)
0353         if (x < bins[i])
0354           return i - 1;
0355       return bins.size() - 1;
0356     }
0357 
0358     /// @brief VBF topolog selection
0359     /// 0 = fail loose selction: m_jj > 400 GeV and Dy_jj > 2.8
0360     /// 1 pass loose, but fail additional cut pT(Hjj)<25. 2 pass tight selection
0361     int vbfTopology(const Jets &jets, const Particle &higgs) {
0362       if (jets.size() < 2)
0363         return 0;
0364       const FourMomentum &j1 = jets[0].momentum(), &j2 = jets[1].momentum();
0365       bool VBFtopo = (j1 + j2).mass() > 400.0 && std::abs(j1.rapidity() - j2.rapidity()) > 2.8;
0366       return VBFtopo ? (j1 + j2 + higgs.momentum()).pt() < 25 ? 2 : 1 : 0;
0367     }
0368 
0369     /// @brief VBF topology selection Stage1.1 and Stage1.2
0370     /// 0 = fail loose selection: m_jj > 350 GeV
0371     /// 1 pass loose, but fail additional cut pT(Hjj)<25. 2 pass pT(Hjj)>25 selection
0372     /// 3 pass tight (m_jj>700 GeV), but fail additional cut pT(Hjj)<25. 4 pass pT(Hjj)>25 selection
0373     int vbfTopology_Stage1_X(const Jets &jets, const Particle &higgs) {
0374       if (jets.size() < 2)
0375         return 0;
0376       const FourMomentum &j1 = jets[0].momentum(), &j2 = jets[1].momentum();
0377       double mjj = (j1 + j2).mass();
0378       if (mjj > 350 && mjj <= 700)
0379         return (j1 + j2 + higgs.momentum()).pt() < 25 ? 1 : 2;
0380       else if (mjj > 700)
0381         return (j1 + j2 + higgs.momentum()).pt() < 25 ? 3 : 4;
0382       else
0383         return 0;
0384     }
0385 
0386     /// @brief VBF topology selection for Stage1.1 and Stage 1.2 Fine
0387     /// 0 = fail loose selection: m_jj > 350 GeV
0388     /// 1 pass loose, but fail additional cut pT(Hjj)<25. 2 pass pT(Hjj)>25 selection
0389     /// 3 pass 700<m_jj<1000 GeV, but fail additional cut pT(Hjj)<25. 4 pass pT(Hjj)>25 selection
0390     /// 5 pass 1000<m_jj<1500 GeV, but fail additional cut pT(Hjj)<25. 6 pass pT(Hjj)>25 selection
0391     /// 7 pass m_jj>1500 GeV, but fail additional cut pT(Hjj)<25. 8 pass pT(Hjj)>25 selection
0392     int vbfTopology_Stage1_X_Fine(const Jets &jets, const Particle &higgs) {
0393       if (jets.size() < 2)
0394         return 0;
0395       const FourMomentum &j1 = jets[0].momentum(), &j2 = jets[1].momentum();
0396       double mjj = (j1 + j2).mass();
0397       if (mjj > 350 && mjj <= 700)
0398         return (j1 + j2 + higgs.momentum()).pt() < 25 ? 1 : 2;
0399       else if (mjj > 700 && mjj <= 1000)
0400         return (j1 + j2 + higgs.momentum()).pt() < 25 ? 3 : 4;
0401       else if (mjj > 1000 && mjj <= 1500)
0402         return (j1 + j2 + higgs.momentum()).pt() < 25 ? 5 : 6;
0403       else if (mjj > 1500)
0404         return (j1 + j2 + higgs.momentum()).pt() < 25 ? 7 : 8;
0405       else
0406         return 0;
0407     }
0408 
0409     /// @brief Whether the Higgs is produced in association with a vector boson (VH)
0410     bool isVH(HTXS::HiggsProdMode p) { return p == HTXS::WH || p == HTXS::QQ2ZH || p == HTXS::GG2ZH; }
0411 
0412     /// @brief Stage-0 HTXS categorization
0413     HTXS::Stage0::Category getStage0Category(const HTXS::HiggsProdMode prodMode,
0414                                              const Particle &higgs,
0415                                              const Particle &V) {
0416       using namespace HTXS::Stage0;
0417       int ctrlHiggs = std::abs(higgs.rapidity()) < 2.5;
0418       // Special cases first, qq→Hqq
0419       if ((prodMode == HTXS::WH || prodMode == HTXS::QQ2ZH) && quarkDecay(V)) {
0420         return ctrlHiggs ? VH2HQQ : VH2HQQ_FWDH;
0421       } else if (prodMode == HTXS::GG2ZH && quarkDecay(V)) {
0422         return Category(HTXS::GGF * 10 + ctrlHiggs);
0423       }
0424       // General case after
0425       return Category(prodMode * 10 + ctrlHiggs);
0426     }
0427 
0428     /// @brief Stage-1 categorization
0429     HTXS::Stage1::Category getStage1Category(const HTXS::HiggsProdMode prodMode,
0430                                              const Particle &higgs,
0431                                              const Jets &jets,
0432                                              const Particle &V) {
0433       using namespace HTXS::Stage1;
0434       int Njets = jets.size(), ctrlHiggs = std::abs(higgs.rapidity()) < 2.5, fwdHiggs = !ctrlHiggs;
0435       double pTj1 = !jets.empty() ? jets[0].momentum().pt() : 0;
0436       int vbfTopo = vbfTopology(jets, higgs);
0437 
0438       // 1. GGF Stage 1 categories
0439       //    Following YR4 write-up: XXXXX
0440       if (prodMode == HTXS::GGF || (prodMode == HTXS::GG2ZH && quarkDecay(V))) {
0441         if (fwdHiggs)
0442           return GG2H_FWDH;
0443         if (Njets == 0)
0444           return GG2H_0J;
0445         else if (Njets == 1)
0446           return Category(GG2H_1J_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200}));
0447         else if (Njets >= 2) {
0448           // events with pT_H>200 get priority over VBF cuts
0449           if (higgs.pt() <= 200) {
0450             if (vbfTopo == 2)
0451               return GG2H_VBFTOPO_JET3VETO;
0452             else if (vbfTopo == 1)
0453               return GG2H_VBFTOPO_JET3;
0454           }
0455           // Njets >= 2jets without VBF topology
0456           return Category(GG2H_GE2J_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200}));
0457         }
0458       }
0459       // 2. Electroweak qq->Hqq Stage 1 categories
0460       else if (prodMode == HTXS::VBF || (isVH(prodMode) && quarkDecay(V))) {
0461         if (std::abs(higgs.rapidity()) > 2.5)
0462           return QQ2HQQ_FWDH;
0463         if (pTj1 > 200)
0464           return QQ2HQQ_PTJET1_GT200;
0465         if (vbfTopo == 2)
0466           return QQ2HQQ_VBFTOPO_JET3VETO;
0467         if (vbfTopo == 1)
0468           return QQ2HQQ_VBFTOPO_JET3;
0469         double mjj = jets.size() > 1 ? (jets[0].mom() + jets[1].mom()).mass() : 0;
0470         if (60 < mjj && mjj < 120)
0471           return QQ2HQQ_VH2JET;
0472         return QQ2HQQ_REST;
0473       }
0474       // 3. WH->Hlv categories
0475       else if (prodMode == HTXS::WH) {
0476         if (fwdHiggs)
0477           return QQ2HLNU_FWDH;
0478         else if (V.pt() < 150)
0479           return QQ2HLNU_PTV_0_150;
0480         else if (V.pt() > 250)
0481           return QQ2HLNU_PTV_GT250;
0482         // 150 < pTV/GeV < 250
0483         return jets.empty() ? QQ2HLNU_PTV_150_250_0J : QQ2HLNU_PTV_150_250_GE1J;
0484       }
0485       // 4. qq->ZH->llH categories
0486       else if (prodMode == HTXS::QQ2ZH) {
0487         if (fwdHiggs)
0488           return QQ2HLL_FWDH;
0489         else if (V.pt() < 150)
0490           return QQ2HLL_PTV_0_150;
0491         else if (V.pt() > 250)
0492           return QQ2HLL_PTV_GT250;
0493         // 150 < pTV/GeV < 250
0494         return jets.empty() ? QQ2HLL_PTV_150_250_0J : QQ2HLL_PTV_150_250_GE1J;
0495       }
0496       // 5. gg->ZH->llH categories
0497       else if (prodMode == HTXS::GG2ZH) {
0498         if (fwdHiggs)
0499           return GG2HLL_FWDH;
0500         if (V.pt() < 150)
0501           return GG2HLL_PTV_0_150;
0502         else if (jets.empty())
0503           return GG2HLL_PTV_GT150_0J;
0504         return GG2HLL_PTV_GT150_GE1J;
0505       }
0506       // 6.ttH,bbH,tH categories
0507       else if (prodMode == HTXS::TTH)
0508         return Category(TTH_FWDH + ctrlHiggs);
0509       else if (prodMode == HTXS::BBH)
0510         return Category(BBH_FWDH + ctrlHiggs);
0511       else if (prodMode == HTXS::TH)
0512         return Category(TH_FWDH + ctrlHiggs);
0513       return UNKNOWN;
0514     }
0515 
0516     /// @brief Stage-1.1 categorization
0517     HTXS::Stage1_1::Category getStage1_1_Category(const HTXS::HiggsProdMode prodMode,
0518                                                   const Particle &higgs,
0519                                                   const Jets &jets,
0520                                                   const Particle &V) {
0521       using namespace HTXS::Stage1_1;
0522       int Njets = jets.size(), ctrlHiggs = std::abs(higgs.rapidity()) < 2.5, fwdHiggs = !ctrlHiggs;
0523       int vbfTopo = vbfTopology_Stage1_X(jets, higgs);
0524 
0525       // 1. GGF Stage 1 categories
0526       //    Following YR4 write-up: XXXXX
0527       if (prodMode == HTXS::GGF || (prodMode == HTXS::GG2ZH && quarkDecay(V))) {
0528         if (fwdHiggs)
0529           return GG2H_FWDH;
0530         if (higgs.pt() > 200)
0531           return GG2H_PTH_GT200;
0532         if (Njets == 0)
0533           return higgs.pt() < 10 ? GG2H_0J_PTH_0_10 : GG2H_0J_PTH_GT10;
0534         if (Njets == 1)
0535           return Category(GG2H_1J_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200}));
0536         if (Njets > 1) {
0537           //VBF topology
0538           if (vbfTopo)
0539             return Category(GG2H_MJJ_350_700_PTHJJ_0_25 + vbfTopo - 1);
0540           //Njets >= 2jets without VBF topology (mjj<350)
0541           return Category(GG2H_GE2J_MJJ_0_350_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200}));
0542         }
0543       }
0544 
0545       // 2. Electroweak qq->Hqq Stage 1.1 categories
0546       else if (prodMode == HTXS::VBF || (isVH(prodMode) && quarkDecay(V))) {
0547         if (std::abs(higgs.rapidity()) > 2.5)
0548           return QQ2HQQ_FWDH;
0549         int Njets = jets.size();
0550         if (Njets == 0)
0551           return QQ2HQQ_0J;
0552         else if (Njets == 1)
0553           return QQ2HQQ_1J;
0554         else if (Njets >= 2) {
0555           double mjj = (jets[0].mom() + jets[1].mom()).mass();
0556           if (mjj < 60)
0557             return QQ2HQQ_MJJ_0_60;
0558           else if (60 < mjj && mjj < 120)
0559             return QQ2HQQ_MJJ_60_120;
0560           else if (120 < mjj && mjj < 350)
0561             return QQ2HQQ_MJJ_120_350;
0562           else if (mjj > 350) {
0563             if (higgs.pt() > 200)
0564               return QQ2HQQ_MJJ_GT350_PTH_GT200;
0565             if (vbfTopo)
0566               return Category(QQ2HQQ_MJJ_GT350_PTH_GT200 + vbfTopo);
0567           }
0568         }
0569       }
0570       // 3. WH->Hlv categories
0571       else if (prodMode == HTXS::WH) {
0572         if (fwdHiggs)
0573           return QQ2HLNU_FWDH;
0574         else if (V.pt() < 75)
0575           return QQ2HLNU_PTV_0_75;
0576         else if (V.pt() < 150)
0577           return QQ2HLNU_PTV_75_150;
0578         else if (V.pt() > 250)
0579           return QQ2HLNU_PTV_GT250;
0580         // 150 < pTV/GeV < 250
0581         return jets.empty() ? QQ2HLNU_PTV_150_250_0J : QQ2HLNU_PTV_150_250_GE1J;
0582       }
0583       // 4. qq->ZH->llH categories
0584       else if (prodMode == HTXS::QQ2ZH) {
0585         if (fwdHiggs)
0586           return QQ2HLL_FWDH;
0587         else if (V.pt() < 75)
0588           return QQ2HLL_PTV_0_75;
0589         else if (V.pt() < 150)
0590           return QQ2HLL_PTV_75_150;
0591         else if (V.pt() > 250)
0592           return QQ2HLL_PTV_GT250;
0593         // 150 < pTV/GeV < 250
0594         return jets.empty() ? QQ2HLL_PTV_150_250_0J : QQ2HLL_PTV_150_250_GE1J;
0595       }
0596       // 5. gg->ZH->llH categories
0597       else if (prodMode == HTXS::GG2ZH) {
0598         if (fwdHiggs)
0599           return GG2HLL_FWDH;
0600         else if (V.pt() < 75)
0601           return GG2HLL_PTV_0_75;
0602         else if (V.pt() < 150)
0603           return GG2HLL_PTV_75_150;
0604         else if (V.pt() > 250)
0605           return GG2HLL_PTV_GT250;
0606         return jets.empty() ? GG2HLL_PTV_150_250_0J : GG2HLL_PTV_150_250_GE1J;
0607       }
0608       // 6.ttH,bbH,tH categories
0609       else if (prodMode == HTXS::TTH)
0610         return Category(TTH_FWDH + ctrlHiggs);
0611       else if (prodMode == HTXS::BBH)
0612         return Category(BBH_FWDH + ctrlHiggs);
0613       else if (prodMode == HTXS::TH)
0614         return Category(TH_FWDH + ctrlHiggs);
0615       return UNKNOWN;
0616     }
0617 
0618     /// @brief Stage-1_1 categorization
0619     HTXS::Stage1_1_Fine::Category getStage1_1_Fine_Category(const HTXS::HiggsProdMode prodMode,
0620                                                             const Particle &higgs,
0621                                                             const Jets &jets,
0622                                                             const Particle &V) {
0623       using namespace HTXS::Stage1_1_Fine;
0624       int Njets = jets.size(), ctrlHiggs = std::abs(higgs.rapidity()) < 2.5, fwdHiggs = !ctrlHiggs;
0625       int vbfTopo = vbfTopology_Stage1_X_Fine(jets, higgs);
0626 
0627       // 1. GGF Stage 1.1 categories
0628       //    Following YR4 write-up: XXXXX
0629       if (prodMode == HTXS::GGF || (prodMode == HTXS::GG2ZH && quarkDecay(V))) {
0630         if (fwdHiggs)
0631           return GG2H_FWDH;
0632         if (higgs.pt() > 200)
0633           return GG2H_PTH_GT200;
0634         if (Njets == 0)
0635           return higgs.pt() < 10 ? GG2H_0J_PTH_0_10 : GG2H_0J_PTH_GT10;
0636         if (Njets == 1)
0637           return Category(GG2H_1J_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200}));
0638         if (Njets > 1) {
0639           //double mjj = (jets[0].mom()+jets[1].mom()).mass();
0640           double pTHjj = (jets[0].momentum() + jets[1].momentum() + higgs.momentum()).pt();
0641           //VBF topology
0642           if (vbfTopo)
0643             return Category(GG2H_MJJ_350_700_PTHJJ_0_25 + vbfTopo - 1);
0644           //Njets >= 2jets without VBF topology (mjj<350)
0645           if (pTHjj < 25)
0646             return Category(GG2H_GE2J_MJJ_0_350_PTH_0_60_PTHJJ_0_25 + getBin(higgs.pt(), {0, 60, 120, 200}));
0647           else
0648             return Category(GG2H_GE2J_MJJ_0_350_PTH_0_60_PTHJJ_GT25 + getBin(higgs.pt(), {0, 60, 120, 200}));
0649         }
0650       }
0651 
0652       // 2. Electroweak qq->Hqq Stage 1.1 categories
0653       else if (prodMode == HTXS::VBF || (isVH(prodMode) && quarkDecay(V))) {
0654         if (std::abs(higgs.rapidity()) > 2.5)
0655           return QQ2HQQ_FWDH;
0656         int Njets = jets.size();
0657         if (Njets == 0)
0658           return QQ2HQQ_0J;
0659         else if (Njets == 1)
0660           return QQ2HQQ_1J;
0661         else if (Njets >= 2) {
0662           double mjj = (jets[0].mom() + jets[1].mom()).mass();
0663           double pTHjj = (jets[0].momentum() + jets[1].momentum() + higgs.momentum()).pt();
0664           if (mjj < 350) {
0665             if (pTHjj < 25)
0666               return Category(QQ2HQQ_MJJ_0_60_PTHJJ_0_25 + getBin(mjj, {0, 60, 120, 350}));
0667             else
0668               return Category(QQ2HQQ_MJJ_0_60_PTHJJ_GT25 + getBin(mjj, {0, 60, 120, 350}));
0669           } else {  //mjj>350 GeV
0670             if (higgs.pt() < 200) {
0671               return Category(QQ2HQQ_MJJ_350_700_PTHJJ_0_25 + vbfTopo - 1);
0672             } else {
0673               return Category(QQ2HQQ_PTH_GT200_MJJ_350_700_PTHJJ_0_25 + vbfTopo - 1);
0674             }
0675           }
0676         }
0677       }
0678       // 3. WH->Hlv categories
0679       else if (prodMode == HTXS::WH) {
0680         if (fwdHiggs)
0681           return QQ2HLNU_FWDH;
0682         int Njets = jets.size();
0683         if (Njets == 0)
0684           return Category(QQ2HLNU_PTV_0_75_0J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0685         if (Njets == 1)
0686           return Category(QQ2HLNU_PTV_0_75_1J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0687         return Category(QQ2HLNU_PTV_0_75_GE2J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0688       }
0689       // 4. qq->ZH->llH categories
0690       else if (prodMode == HTXS::QQ2ZH) {
0691         if (fwdHiggs)
0692           return QQ2HLL_FWDH;
0693         int Njets = jets.size();
0694         if (Njets == 0)
0695           return Category(QQ2HLL_PTV_0_75_0J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0696         if (Njets == 1)
0697           return Category(QQ2HLL_PTV_0_75_1J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0698         return Category(QQ2HLL_PTV_0_75_GE2J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0699       }
0700       // 5. gg->ZH->llH categories
0701       else if (prodMode == HTXS::GG2ZH) {
0702         if (fwdHiggs)
0703           return GG2HLL_FWDH;
0704         int Njets = jets.size();
0705         if (Njets == 0)
0706           return Category(GG2HLL_PTV_0_75_0J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0707         if (Njets == 1)
0708           return Category(GG2HLL_PTV_0_75_1J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0709         return Category(GG2HLL_PTV_0_75_GE2J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0710       }
0711       // 6.ttH,bbH,tH categories
0712       else if (prodMode == HTXS::TTH)
0713         return Category(TTH_FWDH + ctrlHiggs);
0714       else if (prodMode == HTXS::BBH)
0715         return Category(BBH_FWDH + ctrlHiggs);
0716       else if (prodMode == HTXS::TH)
0717         return Category(TH_FWDH + ctrlHiggs);
0718       return UNKNOWN;
0719     }
0720 
0721     /// @brief Stage-1.2 categorization
0722     HTXS::Stage1_2::Category getStage1_2_Category(const HTXS::HiggsProdMode prodMode,
0723                                                   const Particle &higgs,
0724                                                   const Jets &jets,
0725                                                   const Particle &V) {
0726       using namespace HTXS::Stage1_2;
0727       int Njets = jets.size(), ctrlHiggs = std::abs(higgs.rapidity()) < 2.5, fwdHiggs = !ctrlHiggs;
0728       int vbfTopo = vbfTopology_Stage1_X(jets, higgs);
0729 
0730       // 1. GGF Stage 1 categories
0731       //    Following YR4 write-up: XXXXX
0732       if (prodMode == HTXS::GGF || (prodMode == HTXS::GG2ZH && quarkDecay(V))) {
0733         if (fwdHiggs)
0734           return GG2H_FWDH;
0735         if (higgs.pt() > 200)
0736           return Category(GG2H_PTH_200_300 + getBin(higgs.pt(), {200, 300, 450, 650}));
0737         if (Njets == 0)
0738           return higgs.pt() < 10 ? GG2H_0J_PTH_0_10 : GG2H_0J_PTH_GT10;
0739         if (Njets == 1)
0740           return Category(GG2H_1J_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200}));
0741         if (Njets > 1) {
0742           //VBF topology
0743           if (vbfTopo)
0744             return Category(GG2H_GE2J_MJJ_350_700_PTH_0_200_PTHJJ_0_25 + vbfTopo - 1);
0745           //Njets >= 2jets without VBF topology (mjj<350)
0746           return Category(GG2H_GE2J_MJJ_0_350_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200}));
0747         }
0748       }
0749 
0750       // 2. Electroweak qq->Hqq Stage 1.2 categories
0751       else if (prodMode == HTXS::VBF || (isVH(prodMode) && quarkDecay(V))) {
0752         if (std::abs(higgs.rapidity()) > 2.5)
0753           return QQ2HQQ_FWDH;
0754         int Njets = jets.size();
0755         if (Njets == 0)
0756           return QQ2HQQ_0J;
0757         else if (Njets == 1)
0758           return QQ2HQQ_1J;
0759         else if (Njets >= 2) {
0760           double mjj = (jets[0].mom() + jets[1].mom()).mass();
0761           if (mjj < 60)
0762             return QQ2HQQ_GE2J_MJJ_0_60;
0763           else if (60 < mjj && mjj < 120)
0764             return QQ2HQQ_GE2J_MJJ_60_120;
0765           else if (120 < mjj && mjj < 350)
0766             return QQ2HQQ_GE2J_MJJ_120_350;
0767           else if (mjj > 350) {
0768             if (higgs.pt() > 200)
0769               return QQ2HQQ_GE2J_MJJ_GT350_PTH_GT200;
0770             if (vbfTopo)
0771               return Category(QQ2HQQ_GE2J_MJJ_GT350_PTH_GT200 + vbfTopo);
0772           }
0773         }
0774       }
0775       // 3. WH->Hlv categories
0776       else if (prodMode == HTXS::WH) {
0777         if (fwdHiggs)
0778           return QQ2HLNU_FWDH;
0779         else if (V.pt() < 75)
0780           return QQ2HLNU_PTV_0_75;
0781         else if (V.pt() < 150)
0782           return QQ2HLNU_PTV_75_150;
0783         else if (V.pt() > 250)
0784           return QQ2HLNU_PTV_GT250;
0785         // 150 < pTV/GeV < 250
0786         return jets.empty() ? QQ2HLNU_PTV_150_250_0J : QQ2HLNU_PTV_150_250_GE1J;
0787       }
0788       // 4. qq->ZH->llH categories
0789       else if (prodMode == HTXS::QQ2ZH) {
0790         if (fwdHiggs)
0791           return QQ2HLL_FWDH;
0792         else if (V.pt() < 75)
0793           return QQ2HLL_PTV_0_75;
0794         else if (V.pt() < 150)
0795           return QQ2HLL_PTV_75_150;
0796         else if (V.pt() > 250)
0797           return QQ2HLL_PTV_GT250;
0798         // 150 < pTV/GeV < 250
0799         return jets.empty() ? QQ2HLL_PTV_150_250_0J : QQ2HLL_PTV_150_250_GE1J;
0800       }
0801       // 5. gg->ZH->llH categories
0802       else if (prodMode == HTXS::GG2ZH) {
0803         if (fwdHiggs)
0804           return GG2HLL_FWDH;
0805         else if (V.pt() < 75)
0806           return GG2HLL_PTV_0_75;
0807         else if (V.pt() < 150)
0808           return GG2HLL_PTV_75_150;
0809         else if (V.pt() > 250)
0810           return GG2HLL_PTV_GT250;
0811         return jets.empty() ? GG2HLL_PTV_150_250_0J : GG2HLL_PTV_150_250_GE1J;
0812       }
0813       // 6.ttH,bbH,tH categories
0814       else if (prodMode == HTXS::TTH) {
0815         if (fwdHiggs)
0816           return TTH_FWDH;
0817         else
0818           return Category(TTH_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200, 300}));
0819       } else if (prodMode == HTXS::BBH)
0820         return Category(BBH_FWDH + ctrlHiggs);
0821       else if (prodMode == HTXS::TH)
0822         return Category(TH_FWDH + ctrlHiggs);
0823       return UNKNOWN;
0824     }
0825 
0826     /// @brief Stage-1.2 Fine categorization
0827     HTXS::Stage1_2_Fine::Category getStage1_2_Fine_Category(const HTXS::HiggsProdMode prodMode,
0828                                                             const Particle &higgs,
0829                                                             const Jets &jets,
0830                                                             const Particle &V) {
0831       using namespace HTXS::Stage1_2_Fine;
0832       int Njets = jets.size(), ctrlHiggs = std::abs(higgs.rapidity()) < 2.5, fwdHiggs = !ctrlHiggs;
0833       int vbfTopo = vbfTopology_Stage1_X_Fine(jets, higgs);
0834 
0835       // 1. GGF Stage 1.2 categories
0836       //    Following YR4 write-up: XXXXX
0837       if (prodMode == HTXS::GGF || (prodMode == HTXS::GG2ZH && quarkDecay(V))) {
0838         if (fwdHiggs)
0839           return GG2H_FWDH;
0840         if (higgs.pt() > 200) {
0841           if (Njets > 0) {
0842             double pTHj = (jets[0].momentum() + higgs.momentum()).pt();
0843             if (pTHj / higgs.pt() > 0.15)
0844               return Category(GG2H_PTH_200_300_PTHJoverPTH_GT15 + getBin(higgs.pt(), {200, 300, 450, 650}));
0845             else
0846               return Category(GG2H_PTH_200_300_PTHJoverPTH_0_15 + getBin(higgs.pt(), {200, 300, 450, 650}));
0847           } else
0848             return Category(GG2H_PTH_200_300_PTHJoverPTH_0_15 + getBin(higgs.pt(), {200, 300, 450, 650}));
0849         }
0850         if (Njets == 0)
0851           return higgs.pt() < 10 ? GG2H_0J_PTH_0_10 : GG2H_0J_PTH_GT10;
0852         if (Njets == 1)
0853           return Category(GG2H_1J_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200}));
0854         if (Njets > 1) {
0855           //double mjj = (jets[0].mom()+jets[1].mom()).mass();
0856           double pTHjj = (jets[0].momentum() + jets[1].momentum() + higgs.momentum()).pt();
0857           //VBF topology
0858           if (vbfTopo)
0859             return Category(GG2H_GE2J_MJJ_350_700_PTH_0_200_PTHJJ_0_25 + vbfTopo - 1);
0860           //Njets >= 2jets without VBF topology (mjj<350)
0861           if (pTHjj < 25)
0862             return Category(GG2H_GE2J_MJJ_0_350_PTH_0_60_PTHJJ_0_25 + getBin(higgs.pt(), {0, 60, 120, 200}));
0863           else
0864             return Category(GG2H_GE2J_MJJ_0_350_PTH_0_60_PTHJJ_GT25 + getBin(higgs.pt(), {0, 60, 120, 200}));
0865         }
0866       }
0867 
0868       // 2. Electroweak qq->Hqq Stage 1.2 categories
0869       else if (prodMode == HTXS::VBF || (isVH(prodMode) && quarkDecay(V))) {
0870         if (std::abs(higgs.rapidity()) > 2.5)
0871           return QQ2HQQ_FWDH;
0872         int Njets = jets.size();
0873         if (Njets == 0)
0874           return QQ2HQQ_0J;
0875         else if (Njets == 1)
0876           return QQ2HQQ_1J;
0877         else if (Njets >= 2) {
0878           double mjj = (jets[0].mom() + jets[1].mom()).mass();
0879           double pTHjj = (jets[0].momentum() + jets[1].momentum() + higgs.momentum()).pt();
0880           if (mjj < 350) {
0881             if (pTHjj < 25)
0882               return Category(QQ2HQQ_GE2J_MJJ_0_60_PTHJJ_0_25 + getBin(mjj, {0, 60, 120, 350}));
0883             else
0884               return Category(QQ2HQQ_GE2J_MJJ_0_60_PTHJJ_GT25 + getBin(mjj, {0, 60, 120, 350}));
0885           } else {  //mjj>350 GeV
0886             if (higgs.pt() < 200) {
0887               return Category(QQ2HQQ_GE2J_MJJ_350_700_PTH_0_200_PTHJJ_0_25 + vbfTopo - 1);
0888             } else {
0889               return Category(QQ2HQQ_GE2J_MJJ_350_700_PTH_GT200_PTHJJ_0_25 + vbfTopo - 1);
0890             }
0891           }
0892         }
0893       }
0894       // 3. WH->Hlv categories
0895       else if (prodMode == HTXS::WH) {
0896         if (fwdHiggs)
0897           return QQ2HLNU_FWDH;
0898         int Njets = jets.size();
0899         if (Njets == 0)
0900           return Category(QQ2HLNU_PTV_0_75_0J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0901         if (Njets == 1)
0902           return Category(QQ2HLNU_PTV_0_75_1J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0903         return Category(QQ2HLNU_PTV_0_75_GE2J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0904       }
0905       // 4. qq->ZH->llH categories
0906       else if (prodMode == HTXS::QQ2ZH) {
0907         if (fwdHiggs)
0908           return QQ2HLL_FWDH;
0909         int Njets = jets.size();
0910         if (Njets == 0)
0911           return Category(QQ2HLL_PTV_0_75_0J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0912         if (Njets == 1)
0913           return Category(QQ2HLL_PTV_0_75_1J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0914         return Category(QQ2HLL_PTV_0_75_GE2J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0915       }
0916       // 5. gg->ZH->llH categories
0917       else if (prodMode == HTXS::GG2ZH) {
0918         if (fwdHiggs)
0919           return GG2HLL_FWDH;
0920         int Njets = jets.size();
0921         if (Njets == 0)
0922           return Category(GG2HLL_PTV_0_75_0J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0923         if (Njets == 1)
0924           return Category(GG2HLL_PTV_0_75_1J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0925         return Category(GG2HLL_PTV_0_75_GE2J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0926       }
0927       // 6.ttH,bbH,tH categories
0928       else if (prodMode == HTXS::TTH) {
0929         if (fwdHiggs)
0930           return TTH_FWDH;
0931         else
0932           return Category(TTH_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200, 300, 450}));
0933       } else if (prodMode == HTXS::BBH)
0934         return Category(BBH_FWDH + ctrlHiggs);
0935       else if (prodMode == HTXS::TH)
0936         return Category(TH_FWDH + ctrlHiggs);
0937       return UNKNOWN;
0938     }
0939 
0940     /// @}
0941 
0942     /// @name Default Rivet analysis methods and steering methods
0943     /// @{
0944 
0945     /// @brief Sets the Higgs production mode
0946     void setHiggsProdMode(HTXS::HiggsProdMode prodMode) { m_HiggsProdMode = prodMode; }
0947 
0948     /// @brief default Rivet Analysis::init method
0949     /// Booking of histograms, initializing Rivet projection
0950     /// Extracts Higgs production mode from shell variable if not set manually using setHiggsProdMode
0951     void init() override {
0952       printf("==============================================================\n");
0953       printf("========     HiggsTemplateCrossSections Initialization     =========\n");
0954       printf("==============================================================\n");
0955       // check that the production mode has been set
0956       // if running in standalone Rivet the production mode is set through an env variable
0957       if (m_HiggsProdMode == HTXS::UNKNOWN) {
0958         char *pm_env = std::getenv("HIGGSPRODMODE");
0959         string pm(pm_env == nullptr ? "" : pm_env);
0960         if (pm == "GGF")
0961           m_HiggsProdMode = HTXS::GGF;
0962         else if (pm == "VBF")
0963           m_HiggsProdMode = HTXS::VBF;
0964         else if (pm == "WH")
0965           m_HiggsProdMode = HTXS::WH;
0966         else if (pm == "ZH")
0967           m_HiggsProdMode = HTXS::QQ2ZH;
0968         else if (pm == "QQ2ZH")
0969           m_HiggsProdMode = HTXS::QQ2ZH;
0970         else if (pm == "GG2ZH")
0971           m_HiggsProdMode = HTXS::GG2ZH;
0972         else if (pm == "TTH")
0973           m_HiggsProdMode = HTXS::TTH;
0974         else if (pm == "BBH")
0975           m_HiggsProdMode = HTXS::BBH;
0976         else if (pm == "TH")
0977           m_HiggsProdMode = HTXS::TH;
0978         else {
0979           MSG_WARNING("No HIGGSPRODMODE shell variable found. Needed when running Rivet stand-alone.");
0980         }
0981       }
0982 
0983       // Projections for final state particles
0984       const FinalState FS;
0985       declare(FS, "FS");
0986 
0987       // initialize the histograms with for each of the stages
0988       initializeHistos();
0989       m_sumw = 0.0;
0990       m_sumevents = 0;
0991       printf("==============================================================\n");
0992       printf("========             Higgs prod mode %d              =========\n", m_HiggsProdMode);
0993       printf("========          Sucessful Initialization           =========\n");
0994       printf("==============================================================\n");
0995     }
0996 
0997     // Perform the per-event analysis
0998     void analyze(const Event &event) override {
0999       // get the classification
1000       HiggsClassification cat = classifyEvent(event, m_HiggsProdMode);
1001 
1002       // Fill histograms: categorization --> linerize the categories
1003       const double weight = 1.0;
1004       m_sumw += weight;
1005 
1006       int F = cat.stage0_cat % 10, P = cat.stage1_cat_pTjet30GeV / 100;
1007       hist_stage0->fill(cat.stage0_cat / 10 * 2 + F, weight);
1008 
1009       // Stage 1 enum offsets for each production mode: GGF=12, VBF=6, WH= 5, QQ2ZH=5, GG2ZH=4, TTH=2, BBH=2, TH=2
1010       vector<int> offset({0, 1, 13, 19, 24, 29, 33, 35, 37, 39});
1011       int off = offset[P];
1012       hist_stage1_pTjet25->fill(cat.stage1_cat_pTjet25GeV % 100 + off, weight);
1013       hist_stage1_pTjet30->fill(cat.stage1_cat_pTjet30GeV % 100 + off, weight);
1014 
1015       // Stage 1_2 enum offsets for each production mode: GGF=17, VBF=11, WH= 6, QQ2ZH=6, GG2ZH=6, TTH=6, BBH=2, TH=2
1016       static vector<int> offset1_2({0, 1, 18, 29, 35, 41, 47, 53, 55, 57});
1017       int off1_2 = offset1_2[P];
1018       // Stage 1_2 Fine enum offsets for each production mode: GGF=28, VBF=25, WH= 16, QQ2ZH=16, GG2ZH=16, TTH=7, BBH=2, TH=2
1019       static vector<int> offset1_2f({0, 1, 29, 54, 70, 86, 102, 109, 111, 113});
1020       int off1_2f = offset1_2f[P];
1021       hist_stage1_2_pTjet25->fill(cat.stage1_2_cat_pTjet25GeV % 100 + off1_2, weight);
1022       hist_stage1_2_pTjet30->fill(cat.stage1_2_cat_pTjet30GeV % 100 + off1_2, weight);
1023       hist_stage1_2_fine_pTjet25->fill(cat.stage1_2_fine_cat_pTjet25GeV % 100 + off1_2f, weight);
1024       hist_stage1_2_fine_pTjet30->fill(cat.stage1_2_fine_cat_pTjet30GeV % 100 + off1_2f, weight);
1025 
1026       // Fill histograms: variables used in the categorization
1027       hist_pT_Higgs->fill(cat.higgs.pT(), weight);
1028       hist_y_Higgs->fill(cat.higgs.rapidity(), weight);
1029       hist_pT_V->fill(cat.V.pT(), weight);
1030 
1031       hist_Njets25->fill(cat.jets25.size(), weight);
1032       hist_Njets30->fill(cat.jets30.size(), weight);
1033 
1034       hist_isZ2vv->fill(cat.isZ2vvDecay, weight);
1035 
1036       // Jet variables. Use jet collection with pT threshold at 30 GeV
1037       if (!cat.jets30.empty())
1038         hist_pT_jet1->fill(cat.jets30[0].pt(), weight);
1039       if (cat.jets30.size() >= 2) {
1040         const FourMomentum &j1 = cat.jets30[0].momentum(), &j2 = cat.jets30[1].momentum();
1041         hist_deltay_jj->fill(std::abs(j1.rapidity() - j2.rapidity()), weight);
1042         hist_dijet_mass->fill((j1 + j2).mass(), weight);
1043         hist_pT_Hjj->fill((j1 + j2 + cat.higgs.momentum()).pt(), weight);
1044       }
1045     }
1046 
1047     void printClassificationSummary() {
1048       MSG_INFO(" ====================================================== ");
1049       MSG_INFO("      Higgs Template X-Sec Categorization Tool          ");
1050       MSG_INFO("                Status Code Summary                     ");
1051       MSG_INFO(" ====================================================== ");
1052       bool allSuccess = (m_sumevents == m_errorCount[HTXS::SUCCESS]);
1053       if (allSuccess)
1054         MSG_INFO("     >>>> All " << m_errorCount[HTXS::SUCCESS] << " events successfully categorized!");
1055       else {
1056         MSG_INFO("     >>>> " << m_errorCount[HTXS::SUCCESS] << " events successfully categorized");
1057         MSG_INFO("     >>>> --> the following errors occured:");
1058         MSG_INFO("     >>>> " << m_errorCount[HTXS::PRODMODE_DEFINED] << " had an undefined Higgs production mode.");
1059         MSG_INFO("     >>>> " << m_errorCount[HTXS::MOMENTUM_CONSERVATION] << " failed momentum conservation.");
1060         MSG_INFO("     >>>> " << m_errorCount[HTXS::HIGGS_IDENTIFICATION]
1061                               << " failed to identify a valid Higgs boson.");
1062         MSG_INFO("     >>>> " << m_errorCount[HTXS::HS_VTX_IDENTIFICATION]
1063                               << " failed to identify the hard scatter vertex.");
1064         MSG_INFO("     >>>> " << m_errorCount[HTXS::VH_IDENTIFICATION] << " VH: to identify a valid V-boson.");
1065         MSG_INFO("     >>>> " << m_errorCount[HTXS::TOP_W_IDENTIFICATION]
1066                               << " failed to identify valid Ws from top decay.");
1067       }
1068       MSG_INFO(" ====================================================== ");
1069       MSG_INFO(" ====================================================== ");
1070     }
1071 
1072     void finalize() override {
1073       printClassificationSummary();
1074       double sf = m_sumw > 0 ? 1.0 / m_sumw : 1.0;
1075       scale({hist_stage0,
1076              hist_stage1_pTjet25,
1077              hist_stage1_pTjet30,
1078              hist_stage1_2_pTjet25,
1079              hist_stage1_2_pTjet30,
1080              hist_stage1_2_fine_pTjet25,
1081              hist_stage1_2_fine_pTjet30,
1082              hist_Njets25,
1083              hist_Njets30,
1084              hist_pT_Higgs,
1085              hist_y_Higgs,
1086              hist_pT_V,
1087              hist_pT_jet1,
1088              hist_deltay_jj,
1089              hist_dijet_mass,
1090              hist_pT_Hjj},
1091             sf);
1092     }
1093 
1094     /*
1095      *  initialize histograms
1096      */
1097 
1098     void initializeHistos() {
1099       book(hist_stage0, "HTXS_stage0", 20, 0, 20);
1100       book(hist_stage1_pTjet25, "HTXS_stage1_pTjet25", 40, 0, 40);
1101       book(hist_stage1_pTjet30, "HTXS_stage1_pTjet30", 40, 0, 40);
1102       book(hist_stage1_2_pTjet25, "HTXS_stage1_2_pTjet25", 57, 0, 57);
1103       book(hist_stage1_2_pTjet30, "HTXS_stage1_2_pTjet30", 57, 0, 57);
1104       book(hist_stage1_2_fine_pTjet25, "HTXS_stage1_2_fine_pTjet25", 113, 0, 113);
1105       book(hist_stage1_2_fine_pTjet30, "HTXS_stage1_2_fine_pTjet30", 113, 0, 113);
1106       book(hist_pT_Higgs, "pT_Higgs", 80, 0, 400);
1107       book(hist_y_Higgs, "y_Higgs", 80, -4, 4);
1108       book(hist_pT_V, "pT_V", 80, 0, 400);
1109       book(hist_pT_jet1, "pT_jet1", 80, 0, 400);
1110       book(hist_deltay_jj, "deltay_jj", 50, 0, 10);
1111       book(hist_dijet_mass, "m_jj", 50, 0, 2000);
1112       book(hist_pT_Hjj, "pT_Hjj", 50, 0, 250);
1113       book(hist_Njets25, "Njets25", 10, 0, 10);
1114       book(hist_Njets30, "Njets30", 10, 0, 10);
1115       book(hist_isZ2vv, "isZ2vv", 2, 0, 2);
1116     }
1117     /// @}
1118 
1119     /*
1120      *    initialize private members used in the classification procedure
1121      */
1122 
1123   private:
1124     double m_sumw;
1125     size_t m_sumevents;
1126     HTXS::HiggsProdMode m_HiggsProdMode;
1127     std::map<HTXS::ErrorCode, size_t> m_errorCount;
1128     Histo1DPtr hist_stage0;
1129     Histo1DPtr hist_stage1_pTjet25, hist_stage1_pTjet30;
1130     Histo1DPtr hist_stage1_2_pTjet25, hist_stage1_2_pTjet30;
1131     Histo1DPtr hist_stage1_2_fine_pTjet25, hist_stage1_2_fine_pTjet30;
1132     Histo1DPtr hist_pT_Higgs, hist_y_Higgs;
1133     Histo1DPtr hist_pT_V, hist_pT_jet1;
1134     Histo1DPtr hist_deltay_jj, hist_dijet_mass, hist_pT_Hjj;
1135     Histo1DPtr hist_Njets25, hist_Njets30;
1136     Histo1DPtr hist_isZ2vv;
1137   };
1138 
1139   // the PLUGIN only needs to be decleared when running standalone Rivet
1140   // and causes compilation / linking issues if included in Athena / RootCore
1141   //check for Rivet environment variable RIVET_ANALYSIS_PATH
1142 #ifdef RIVET_ANALYSIS_PATH
1143   // The hook for the plugin system
1144   DECLARE_RIVET_PLUGIN(HiggsTemplateCrossSections);
1145 #endif
1146 
1147 }  // namespace Rivet
1148 
1149 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!