Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-01-09 23:33:36 

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 (const 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 (const 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 (const 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 (const 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       // Temporary fix: add variables to perform STXS 1.3 classification with nanoAOD on-the-fly
0306       // Vector-boson pt for VH production modes
0307       if (isVH(prodMode)) {
0308         cat.V_pt = cat.V.pt();
0309       } else {
0310         cat.V_pt = -999;
0311       }
0312       // Dijet variables using jets30 collection
0313       if (cat.jets30.size() >= 2) {
0314         cat.Mjj = (cat.jets30[0].mom() + cat.jets30[1].mom()).mass();
0315         cat.ptHjj = (cat.jets30[0].mom() + cat.jets30[1].mom() + cat.higgs.momentum()).pt();
0316         cat.dPhijj = cat.jets30[0].mom().phi() - cat.jets30[1].mom().phi();
0317         // Return phi angle in the interval [-PI,PI)
0318         if (cat.dPhijj >= Rivet::pi)
0319           cat.dPhijj -= 2 * Rivet::pi;
0320         else if (cat.dPhijj < -1 * Rivet::pi)
0321           cat.dPhijj += 2 * Rivet::pi;
0322       } else {
0323         cat.Mjj = -999;
0324         cat.ptHjj = -999;
0325         cat.dPhijj = -999;
0326       }
0327 
0328       // check that four mometum sum of all stable particles satisfies momentum consevation
0329       /*
0330       if ( sum.pt()>0.1 )
0331     return error(cat,HTXS::MOMENTUM_CONSERVATION,"Four vector sum does not amount to pT=0, m=E=sqrt(s), but pT="+
0332              std::to_string(sum.pt())+" GeV and m = "+std::to_string(sum.mass())+" GeV");
0333 */
0334       // check if V-boson was not included in the event record but decay particles were
0335       // EFT contact interaction: return UNKNOWN for category but set all event/particle kinematics
0336       if (is_uncatdV)
0337         return error(cat, HTXS::VH_IDENTIFICATION, "Failed to identify associated V-boson!");
0338 
0339       /*****
0340        * Step 4.
0341        *   Classify and save output
0342        */
0343 
0344       // Apply the categorization categorization
0345       cat.isZ2vvDecay = false;
0346       if ((prodMode == HTXS::GG2ZH || prodMode == HTXS::QQ2ZH) && !quarkDecay(cat.V) && !ChLeptonDecay(cat.V))
0347         cat.isZ2vvDecay = true;
0348       cat.stage0_cat = getStage0Category(prodMode, cat.higgs, cat.V);
0349       cat.stage1_cat_pTjet25GeV = getStage1Category(prodMode, cat.higgs, cat.jets25, cat.V);
0350       cat.stage1_cat_pTjet30GeV = getStage1Category(prodMode, cat.higgs, cat.jets30, cat.V);
0351       cat.stage1_1_cat_pTjet25GeV = getStage1_1_Category(prodMode, cat.higgs, cat.jets25, cat.V);
0352       cat.stage1_1_cat_pTjet30GeV = getStage1_1_Category(prodMode, cat.higgs, cat.jets30, cat.V);
0353       cat.stage1_1_fine_cat_pTjet25GeV = getStage1_1_Fine_Category(prodMode, cat.higgs, cat.jets25, cat.V);
0354       cat.stage1_1_fine_cat_pTjet30GeV = getStage1_1_Fine_Category(prodMode, cat.higgs, cat.jets30, cat.V);
0355       cat.stage1_2_cat_pTjet25GeV = getStage1_2_Category(prodMode, cat.higgs, cat.jets25, cat.V);
0356       cat.stage1_2_cat_pTjet30GeV = getStage1_2_Category(prodMode, cat.higgs, cat.jets30, cat.V);
0357       cat.stage1_2_fine_cat_pTjet25GeV = getStage1_2_Fine_Category(prodMode, cat.higgs, cat.jets25, cat.V);
0358       cat.stage1_2_fine_cat_pTjet30GeV = getStage1_2_Fine_Category(prodMode, cat.higgs, cat.jets30, cat.V);
0359 
0360       cat.errorCode = HTXS::SUCCESS;
0361       ++m_errorCount[HTXS::SUCCESS];
0362       ++m_sumevents;
0363 
0364       return cat;
0365     }
0366 
0367     /// @name Categorization methods
0368     /// Methods to assign the truth category based
0369     /// on the identified Higgs boson and associated
0370     /// vector bosons and/or reconstructed jets
0371     /// @{
0372 
0373     /// @brief Return bin index of x given the provided bin edges. 0=first bin, -1=underflow bin.
0374     int getBin(double x, std::vector<double> bins) {
0375       for (size_t i = 1; i < bins.size(); ++i)
0376         if (x < bins[i])
0377           return i - 1;
0378       return bins.size() - 1;
0379     }
0380 
0381     /// @brief VBF topolog selection
0382     /// 0 = fail loose selction: m_jj > 400 GeV and Dy_jj > 2.8
0383     /// 1 pass loose, but fail additional cut pT(Hjj)<25. 2 pass tight selection
0384     int vbfTopology(const Jets &jets, const Particle &higgs) {
0385       if (jets.size() < 2)
0386         return 0;
0387       const FourMomentum &j1 = jets[0].momentum(), &j2 = jets[1].momentum();
0388       bool VBFtopo = (j1 + j2).mass() > 400.0 && std::abs(j1.rapidity() - j2.rapidity()) > 2.8;
0389       return VBFtopo ? (j1 + j2 + higgs.momentum()).pt() < 25 ? 2 : 1 : 0;
0390     }
0391 
0392     /// @brief VBF topology selection Stage1.1 and Stage1.2
0393     /// 0 = fail loose selection: m_jj > 350 GeV
0394     /// 1 pass loose, but fail additional cut pT(Hjj)<25. 2 pass pT(Hjj)>25 selection
0395     /// 3 pass tight (m_jj>700 GeV), but fail additional cut pT(Hjj)<25. 4 pass pT(Hjj)>25 selection
0396     int vbfTopology_Stage1_X(const Jets &jets, const Particle &higgs) {
0397       if (jets.size() < 2)
0398         return 0;
0399       const FourMomentum &j1 = jets[0].momentum(), &j2 = jets[1].momentum();
0400       double mjj = (j1 + j2).mass();
0401       if (mjj > 350 && mjj <= 700)
0402         return (j1 + j2 + higgs.momentum()).pt() < 25 ? 1 : 2;
0403       else if (mjj > 700)
0404         return (j1 + j2 + higgs.momentum()).pt() < 25 ? 3 : 4;
0405       else
0406         return 0;
0407     }
0408 
0409     /// @brief VBF topology selection for Stage1.1 and Stage 1.2 Fine
0410     /// 0 = fail loose selection: m_jj > 350 GeV
0411     /// 1 pass loose, but fail additional cut pT(Hjj)<25. 2 pass pT(Hjj)>25 selection
0412     /// 3 pass 700<m_jj<1000 GeV, but fail additional cut pT(Hjj)<25. 4 pass pT(Hjj)>25 selection
0413     /// 5 pass 1000<m_jj<1500 GeV, but fail additional cut pT(Hjj)<25. 6 pass pT(Hjj)>25 selection
0414     /// 7 pass m_jj>1500 GeV, but fail additional cut pT(Hjj)<25. 8 pass pT(Hjj)>25 selection
0415     int vbfTopology_Stage1_X_Fine(const Jets &jets, const Particle &higgs) {
0416       if (jets.size() < 2)
0417         return 0;
0418       const FourMomentum &j1 = jets[0].momentum(), &j2 = jets[1].momentum();
0419       double mjj = (j1 + j2).mass();
0420       if (mjj > 350 && mjj <= 700)
0421         return (j1 + j2 + higgs.momentum()).pt() < 25 ? 1 : 2;
0422       else if (mjj > 700 && mjj <= 1000)
0423         return (j1 + j2 + higgs.momentum()).pt() < 25 ? 3 : 4;
0424       else if (mjj > 1000 && mjj <= 1500)
0425         return (j1 + j2 + higgs.momentum()).pt() < 25 ? 5 : 6;
0426       else if (mjj > 1500)
0427         return (j1 + j2 + higgs.momentum()).pt() < 25 ? 7 : 8;
0428       else
0429         return 0;
0430     }
0431 
0432     /// @brief Whether the Higgs is produced in association with a vector boson (VH)
0433     bool isVH(HTXS::HiggsProdMode p) { return p == HTXS::WH || p == HTXS::QQ2ZH || p == HTXS::GG2ZH; }
0434 
0435     /// @brief Stage-0 HTXS categorization
0436     HTXS::Stage0::Category getStage0Category(const HTXS::HiggsProdMode prodMode,
0437                                              const Particle &higgs,
0438                                              const Particle &V) {
0439       using namespace HTXS::Stage0;
0440       int ctrlHiggs = std::abs(higgs.rapidity()) < 2.5;
0441       // Special cases first, qq→Hqq
0442       if ((prodMode == HTXS::WH || prodMode == HTXS::QQ2ZH) && quarkDecay(V)) {
0443         return ctrlHiggs ? VH2HQQ : VH2HQQ_FWDH;
0444       } else if (prodMode == HTXS::GG2ZH && quarkDecay(V)) {
0445         return Category(HTXS::GGF * 10 + ctrlHiggs);
0446       }
0447       // General case after
0448       return Category(prodMode * 10 + ctrlHiggs);
0449     }
0450 
0451     /// @brief Stage-1 categorization
0452     HTXS::Stage1::Category getStage1Category(const HTXS::HiggsProdMode prodMode,
0453                                              const Particle &higgs,
0454                                              const Jets &jets,
0455                                              const Particle &V) {
0456       using namespace HTXS::Stage1;
0457       int Njets = jets.size(), ctrlHiggs = std::abs(higgs.rapidity()) < 2.5, fwdHiggs = !ctrlHiggs;
0458       double pTj1 = !jets.empty() ? jets[0].momentum().pt() : 0;
0459       int vbfTopo = vbfTopology(jets, higgs);
0460 
0461       // 1. GGF Stage 1 categories
0462       //    Following YR4 write-up: XXXXX
0463       if (prodMode == HTXS::GGF || (prodMode == HTXS::GG2ZH && quarkDecay(V))) {
0464         if (fwdHiggs)
0465           return GG2H_FWDH;
0466         if (Njets == 0)
0467           return GG2H_0J;
0468         else if (Njets == 1)
0469           return Category(GG2H_1J_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200}));
0470         else if (Njets >= 2) {
0471           // events with pT_H>200 get priority over VBF cuts
0472           if (higgs.pt() <= 200) {
0473             if (vbfTopo == 2)
0474               return GG2H_VBFTOPO_JET3VETO;
0475             else if (vbfTopo == 1)
0476               return GG2H_VBFTOPO_JET3;
0477           }
0478           // Njets >= 2jets without VBF topology
0479           return Category(GG2H_GE2J_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200}));
0480         }
0481       }
0482       // 2. Electroweak qq->Hqq Stage 1 categories
0483       else if (prodMode == HTXS::VBF || (isVH(prodMode) && quarkDecay(V))) {
0484         if (std::abs(higgs.rapidity()) > 2.5)
0485           return QQ2HQQ_FWDH;
0486         if (pTj1 > 200)
0487           return QQ2HQQ_PTJET1_GT200;
0488         if (vbfTopo == 2)
0489           return QQ2HQQ_VBFTOPO_JET3VETO;
0490         if (vbfTopo == 1)
0491           return QQ2HQQ_VBFTOPO_JET3;
0492         double mjj = jets.size() > 1 ? (jets[0].mom() + jets[1].mom()).mass() : 0;
0493         if (60 < mjj && mjj < 120)
0494           return QQ2HQQ_VH2JET;
0495         return QQ2HQQ_REST;
0496       }
0497       // 3. WH->Hlv categories
0498       else if (prodMode == HTXS::WH) {
0499         if (fwdHiggs)
0500           return QQ2HLNU_FWDH;
0501         else if (V.pt() < 150)
0502           return QQ2HLNU_PTV_0_150;
0503         else if (V.pt() > 250)
0504           return QQ2HLNU_PTV_GT250;
0505         // 150 < pTV/GeV < 250
0506         return jets.empty() ? QQ2HLNU_PTV_150_250_0J : QQ2HLNU_PTV_150_250_GE1J;
0507       }
0508       // 4. qq->ZH->llH categories
0509       else if (prodMode == HTXS::QQ2ZH) {
0510         if (fwdHiggs)
0511           return QQ2HLL_FWDH;
0512         else if (V.pt() < 150)
0513           return QQ2HLL_PTV_0_150;
0514         else if (V.pt() > 250)
0515           return QQ2HLL_PTV_GT250;
0516         // 150 < pTV/GeV < 250
0517         return jets.empty() ? QQ2HLL_PTV_150_250_0J : QQ2HLL_PTV_150_250_GE1J;
0518       }
0519       // 5. gg->ZH->llH categories
0520       else if (prodMode == HTXS::GG2ZH) {
0521         if (fwdHiggs)
0522           return GG2HLL_FWDH;
0523         if (V.pt() < 150)
0524           return GG2HLL_PTV_0_150;
0525         else if (jets.empty())
0526           return GG2HLL_PTV_GT150_0J;
0527         return GG2HLL_PTV_GT150_GE1J;
0528       }
0529       // 6.ttH,bbH,tH categories
0530       else if (prodMode == HTXS::TTH)
0531         return Category(TTH_FWDH + ctrlHiggs);
0532       else if (prodMode == HTXS::BBH)
0533         return Category(BBH_FWDH + ctrlHiggs);
0534       else if (prodMode == HTXS::TH)
0535         return Category(TH_FWDH + ctrlHiggs);
0536       return UNKNOWN;
0537     }
0538 
0539     /// @brief Stage-1.1 categorization
0540     HTXS::Stage1_1::Category getStage1_1_Category(const HTXS::HiggsProdMode prodMode,
0541                                                   const Particle &higgs,
0542                                                   const Jets &jets,
0543                                                   const Particle &V) {
0544       using namespace HTXS::Stage1_1;
0545       int Njets = jets.size(), ctrlHiggs = std::abs(higgs.rapidity()) < 2.5, fwdHiggs = !ctrlHiggs;
0546       int vbfTopo = vbfTopology_Stage1_X(jets, higgs);
0547 
0548       // 1. GGF Stage 1 categories
0549       //    Following YR4 write-up: XXXXX
0550       if (prodMode == HTXS::GGF || (prodMode == HTXS::GG2ZH && quarkDecay(V))) {
0551         if (fwdHiggs)
0552           return GG2H_FWDH;
0553         if (higgs.pt() > 200)
0554           return GG2H_PTH_GT200;
0555         if (Njets == 0)
0556           return higgs.pt() < 10 ? GG2H_0J_PTH_0_10 : GG2H_0J_PTH_GT10;
0557         if (Njets == 1)
0558           return Category(GG2H_1J_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200}));
0559         if (Njets > 1) {
0560           //VBF topology
0561           if (vbfTopo)
0562             return Category(GG2H_MJJ_350_700_PTHJJ_0_25 + vbfTopo - 1);
0563           //Njets >= 2jets without VBF topology (mjj<350)
0564           return Category(GG2H_GE2J_MJJ_0_350_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200}));
0565         }
0566       }
0567 
0568       // 2. Electroweak qq->Hqq Stage 1.1 categories
0569       else if (prodMode == HTXS::VBF || (isVH(prodMode) && quarkDecay(V))) {
0570         if (std::abs(higgs.rapidity()) > 2.5)
0571           return QQ2HQQ_FWDH;
0572         int Njets = jets.size();
0573         if (Njets == 0)
0574           return QQ2HQQ_0J;
0575         else if (Njets == 1)
0576           return QQ2HQQ_1J;
0577         else if (Njets >= 2) {
0578           double mjj = (jets[0].mom() + jets[1].mom()).mass();
0579           if (mjj < 60)
0580             return QQ2HQQ_MJJ_0_60;
0581           else if (60 < mjj && mjj < 120)
0582             return QQ2HQQ_MJJ_60_120;
0583           else if (120 < mjj && mjj < 350)
0584             return QQ2HQQ_MJJ_120_350;
0585           else if (mjj > 350) {
0586             if (higgs.pt() > 200)
0587               return QQ2HQQ_MJJ_GT350_PTH_GT200;
0588             if (vbfTopo)
0589               return Category(QQ2HQQ_MJJ_GT350_PTH_GT200 + vbfTopo);
0590           }
0591         }
0592       }
0593       // 3. WH->Hlv categories
0594       else if (prodMode == HTXS::WH) {
0595         if (fwdHiggs)
0596           return QQ2HLNU_FWDH;
0597         else if (V.pt() < 75)
0598           return QQ2HLNU_PTV_0_75;
0599         else if (V.pt() < 150)
0600           return QQ2HLNU_PTV_75_150;
0601         else if (V.pt() > 250)
0602           return QQ2HLNU_PTV_GT250;
0603         // 150 < pTV/GeV < 250
0604         return jets.empty() ? QQ2HLNU_PTV_150_250_0J : QQ2HLNU_PTV_150_250_GE1J;
0605       }
0606       // 4. qq->ZH->llH categories
0607       else if (prodMode == HTXS::QQ2ZH) {
0608         if (fwdHiggs)
0609           return QQ2HLL_FWDH;
0610         else if (V.pt() < 75)
0611           return QQ2HLL_PTV_0_75;
0612         else if (V.pt() < 150)
0613           return QQ2HLL_PTV_75_150;
0614         else if (V.pt() > 250)
0615           return QQ2HLL_PTV_GT250;
0616         // 150 < pTV/GeV < 250
0617         return jets.empty() ? QQ2HLL_PTV_150_250_0J : QQ2HLL_PTV_150_250_GE1J;
0618       }
0619       // 5. gg->ZH->llH categories
0620       else if (prodMode == HTXS::GG2ZH) {
0621         if (fwdHiggs)
0622           return GG2HLL_FWDH;
0623         else if (V.pt() < 75)
0624           return GG2HLL_PTV_0_75;
0625         else if (V.pt() < 150)
0626           return GG2HLL_PTV_75_150;
0627         else if (V.pt() > 250)
0628           return GG2HLL_PTV_GT250;
0629         return jets.empty() ? GG2HLL_PTV_150_250_0J : GG2HLL_PTV_150_250_GE1J;
0630       }
0631       // 6.ttH,bbH,tH categories
0632       else if (prodMode == HTXS::TTH)
0633         return Category(TTH_FWDH + ctrlHiggs);
0634       else if (prodMode == HTXS::BBH)
0635         return Category(BBH_FWDH + ctrlHiggs);
0636       else if (prodMode == HTXS::TH)
0637         return Category(TH_FWDH + ctrlHiggs);
0638       return UNKNOWN;
0639     }
0640 
0641     /// @brief Stage-1_1 categorization
0642     HTXS::Stage1_1_Fine::Category getStage1_1_Fine_Category(const HTXS::HiggsProdMode prodMode,
0643                                                             const Particle &higgs,
0644                                                             const Jets &jets,
0645                                                             const Particle &V) {
0646       using namespace HTXS::Stage1_1_Fine;
0647       int Njets = jets.size(), ctrlHiggs = std::abs(higgs.rapidity()) < 2.5, fwdHiggs = !ctrlHiggs;
0648       int vbfTopo = vbfTopology_Stage1_X_Fine(jets, higgs);
0649 
0650       // 1. GGF Stage 1.1 categories
0651       //    Following YR4 write-up: XXXXX
0652       if (prodMode == HTXS::GGF || (prodMode == HTXS::GG2ZH && quarkDecay(V))) {
0653         if (fwdHiggs)
0654           return GG2H_FWDH;
0655         if (higgs.pt() > 200)
0656           return GG2H_PTH_GT200;
0657         if (Njets == 0)
0658           return higgs.pt() < 10 ? GG2H_0J_PTH_0_10 : GG2H_0J_PTH_GT10;
0659         if (Njets == 1)
0660           return Category(GG2H_1J_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200}));
0661         if (Njets > 1) {
0662           //double mjj = (jets[0].mom()+jets[1].mom()).mass();
0663           double pTHjj = (jets[0].momentum() + jets[1].momentum() + higgs.momentum()).pt();
0664           //VBF topology
0665           if (vbfTopo)
0666             return Category(GG2H_MJJ_350_700_PTHJJ_0_25 + vbfTopo - 1);
0667           //Njets >= 2jets without VBF topology (mjj<350)
0668           if (pTHjj < 25)
0669             return Category(GG2H_GE2J_MJJ_0_350_PTH_0_60_PTHJJ_0_25 + getBin(higgs.pt(), {0, 60, 120, 200}));
0670           else
0671             return Category(GG2H_GE2J_MJJ_0_350_PTH_0_60_PTHJJ_GT25 + getBin(higgs.pt(), {0, 60, 120, 200}));
0672         }
0673       }
0674 
0675       // 2. Electroweak qq->Hqq Stage 1.1 categories
0676       else if (prodMode == HTXS::VBF || (isVH(prodMode) && quarkDecay(V))) {
0677         if (std::abs(higgs.rapidity()) > 2.5)
0678           return QQ2HQQ_FWDH;
0679         int Njets = jets.size();
0680         if (Njets == 0)
0681           return QQ2HQQ_0J;
0682         else if (Njets == 1)
0683           return QQ2HQQ_1J;
0684         else if (Njets >= 2) {
0685           double mjj = (jets[0].mom() + jets[1].mom()).mass();
0686           double pTHjj = (jets[0].momentum() + jets[1].momentum() + higgs.momentum()).pt();
0687           if (mjj < 350) {
0688             if (pTHjj < 25)
0689               return Category(QQ2HQQ_MJJ_0_60_PTHJJ_0_25 + getBin(mjj, {0, 60, 120, 350}));
0690             else
0691               return Category(QQ2HQQ_MJJ_0_60_PTHJJ_GT25 + getBin(mjj, {0, 60, 120, 350}));
0692           } else {  //mjj>350 GeV
0693             if (higgs.pt() < 200) {
0694               return Category(QQ2HQQ_MJJ_350_700_PTHJJ_0_25 + vbfTopo - 1);
0695             } else {
0696               return Category(QQ2HQQ_PTH_GT200_MJJ_350_700_PTHJJ_0_25 + vbfTopo - 1);
0697             }
0698           }
0699         }
0700       }
0701       // 3. WH->Hlv categories
0702       else if (prodMode == HTXS::WH) {
0703         if (fwdHiggs)
0704           return QQ2HLNU_FWDH;
0705         int Njets = jets.size();
0706         if (Njets == 0)
0707           return Category(QQ2HLNU_PTV_0_75_0J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0708         if (Njets == 1)
0709           return Category(QQ2HLNU_PTV_0_75_1J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0710         return Category(QQ2HLNU_PTV_0_75_GE2J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0711       }
0712       // 4. qq->ZH->llH categories
0713       else if (prodMode == HTXS::QQ2ZH) {
0714         if (fwdHiggs)
0715           return QQ2HLL_FWDH;
0716         int Njets = jets.size();
0717         if (Njets == 0)
0718           return Category(QQ2HLL_PTV_0_75_0J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0719         if (Njets == 1)
0720           return Category(QQ2HLL_PTV_0_75_1J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0721         return Category(QQ2HLL_PTV_0_75_GE2J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0722       }
0723       // 5. gg->ZH->llH categories
0724       else if (prodMode == HTXS::GG2ZH) {
0725         if (fwdHiggs)
0726           return GG2HLL_FWDH;
0727         int Njets = jets.size();
0728         if (Njets == 0)
0729           return Category(GG2HLL_PTV_0_75_0J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0730         if (Njets == 1)
0731           return Category(GG2HLL_PTV_0_75_1J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0732         return Category(GG2HLL_PTV_0_75_GE2J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0733       }
0734       // 6.ttH,bbH,tH categories
0735       else if (prodMode == HTXS::TTH)
0736         return Category(TTH_FWDH + ctrlHiggs);
0737       else if (prodMode == HTXS::BBH)
0738         return Category(BBH_FWDH + ctrlHiggs);
0739       else if (prodMode == HTXS::TH)
0740         return Category(TH_FWDH + ctrlHiggs);
0741       return UNKNOWN;
0742     }
0743 
0744     /// @brief Stage-1.2 categorization
0745     HTXS::Stage1_2::Category getStage1_2_Category(const HTXS::HiggsProdMode prodMode,
0746                                                   const Particle &higgs,
0747                                                   const Jets &jets,
0748                                                   const Particle &V) {
0749       using namespace HTXS::Stage1_2;
0750       int Njets = jets.size(), ctrlHiggs = std::abs(higgs.rapidity()) < 2.5, fwdHiggs = !ctrlHiggs;
0751       int vbfTopo = vbfTopology_Stage1_X(jets, higgs);
0752 
0753       // 1. GGF Stage 1 categories
0754       //    Following YR4 write-up: XXXXX
0755       if (prodMode == HTXS::GGF || (prodMode == HTXS::GG2ZH && quarkDecay(V))) {
0756         if (fwdHiggs)
0757           return GG2H_FWDH;
0758         if (higgs.pt() > 200)
0759           return Category(GG2H_PTH_200_300 + getBin(higgs.pt(), {200, 300, 450, 650}));
0760         if (Njets == 0)
0761           return higgs.pt() < 10 ? GG2H_0J_PTH_0_10 : GG2H_0J_PTH_GT10;
0762         if (Njets == 1)
0763           return Category(GG2H_1J_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200}));
0764         if (Njets > 1) {
0765           //VBF topology
0766           if (vbfTopo)
0767             return Category(GG2H_GE2J_MJJ_350_700_PTH_0_200_PTHJJ_0_25 + vbfTopo - 1);
0768           //Njets >= 2jets without VBF topology (mjj<350)
0769           return Category(GG2H_GE2J_MJJ_0_350_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200}));
0770         }
0771       }
0772 
0773       // 2. Electroweak qq->Hqq Stage 1.2 categories
0774       else if (prodMode == HTXS::VBF || (isVH(prodMode) && quarkDecay(V))) {
0775         if (std::abs(higgs.rapidity()) > 2.5)
0776           return QQ2HQQ_FWDH;
0777         int Njets = jets.size();
0778         if (Njets == 0)
0779           return QQ2HQQ_0J;
0780         else if (Njets == 1)
0781           return QQ2HQQ_1J;
0782         else if (Njets >= 2) {
0783           double mjj = (jets[0].mom() + jets[1].mom()).mass();
0784           if (mjj < 60)
0785             return QQ2HQQ_GE2J_MJJ_0_60;
0786           else if (60 < mjj && mjj < 120)
0787             return QQ2HQQ_GE2J_MJJ_60_120;
0788           else if (120 < mjj && mjj < 350)
0789             return QQ2HQQ_GE2J_MJJ_120_350;
0790           else if (mjj > 350) {
0791             if (higgs.pt() > 200)
0792               return QQ2HQQ_GE2J_MJJ_GT350_PTH_GT200;
0793             if (vbfTopo)
0794               return Category(QQ2HQQ_GE2J_MJJ_GT350_PTH_GT200 + vbfTopo);
0795           }
0796         }
0797       }
0798       // 3. WH->Hlv categories
0799       else if (prodMode == HTXS::WH) {
0800         if (fwdHiggs)
0801           return QQ2HLNU_FWDH;
0802         else if (V.pt() < 75)
0803           return QQ2HLNU_PTV_0_75;
0804         else if (V.pt() < 150)
0805           return QQ2HLNU_PTV_75_150;
0806         else if (V.pt() > 250)
0807           return QQ2HLNU_PTV_GT250;
0808         // 150 < pTV/GeV < 250
0809         return jets.empty() ? QQ2HLNU_PTV_150_250_0J : QQ2HLNU_PTV_150_250_GE1J;
0810       }
0811       // 4. qq->ZH->llH categories
0812       else if (prodMode == HTXS::QQ2ZH) {
0813         if (fwdHiggs)
0814           return QQ2HLL_FWDH;
0815         else if (V.pt() < 75)
0816           return QQ2HLL_PTV_0_75;
0817         else if (V.pt() < 150)
0818           return QQ2HLL_PTV_75_150;
0819         else if (V.pt() > 250)
0820           return QQ2HLL_PTV_GT250;
0821         // 150 < pTV/GeV < 250
0822         return jets.empty() ? QQ2HLL_PTV_150_250_0J : QQ2HLL_PTV_150_250_GE1J;
0823       }
0824       // 5. gg->ZH->llH categories
0825       else if (prodMode == HTXS::GG2ZH) {
0826         if (fwdHiggs)
0827           return GG2HLL_FWDH;
0828         else if (V.pt() < 75)
0829           return GG2HLL_PTV_0_75;
0830         else if (V.pt() < 150)
0831           return GG2HLL_PTV_75_150;
0832         else if (V.pt() > 250)
0833           return GG2HLL_PTV_GT250;
0834         return jets.empty() ? GG2HLL_PTV_150_250_0J : GG2HLL_PTV_150_250_GE1J;
0835       }
0836       // 6.ttH,bbH,tH categories
0837       else if (prodMode == HTXS::TTH) {
0838         if (fwdHiggs)
0839           return TTH_FWDH;
0840         else
0841           return Category(TTH_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200, 300}));
0842       } else if (prodMode == HTXS::BBH)
0843         return Category(BBH_FWDH + ctrlHiggs);
0844       else if (prodMode == HTXS::TH)
0845         return Category(TH_FWDH + ctrlHiggs);
0846       return UNKNOWN;
0847     }
0848 
0849     /// @brief Stage-1.2 Fine categorization
0850     HTXS::Stage1_2_Fine::Category getStage1_2_Fine_Category(const HTXS::HiggsProdMode prodMode,
0851                                                             const Particle &higgs,
0852                                                             const Jets &jets,
0853                                                             const Particle &V) {
0854       using namespace HTXS::Stage1_2_Fine;
0855       int Njets = jets.size(), ctrlHiggs = std::abs(higgs.rapidity()) < 2.5, fwdHiggs = !ctrlHiggs;
0856       int vbfTopo = vbfTopology_Stage1_X_Fine(jets, higgs);
0857 
0858       // 1. GGF Stage 1.2 categories
0859       //    Following YR4 write-up: XXXXX
0860       if (prodMode == HTXS::GGF || (prodMode == HTXS::GG2ZH && quarkDecay(V))) {
0861         if (fwdHiggs)
0862           return GG2H_FWDH;
0863         if (higgs.pt() > 200) {
0864           if (Njets > 0) {
0865             double pTHj = (jets[0].momentum() + higgs.momentum()).pt();
0866             if (pTHj / higgs.pt() > 0.15)
0867               return Category(GG2H_PTH_200_300_PTHJoverPTH_GT15 + getBin(higgs.pt(), {200, 300, 450, 650}));
0868             else
0869               return Category(GG2H_PTH_200_300_PTHJoverPTH_0_15 + getBin(higgs.pt(), {200, 300, 450, 650}));
0870           } else
0871             return Category(GG2H_PTH_200_300_PTHJoverPTH_0_15 + getBin(higgs.pt(), {200, 300, 450, 650}));
0872         }
0873         if (Njets == 0)
0874           return higgs.pt() < 10 ? GG2H_0J_PTH_0_10 : GG2H_0J_PTH_GT10;
0875         if (Njets == 1)
0876           return Category(GG2H_1J_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200}));
0877         if (Njets > 1) {
0878           //double mjj = (jets[0].mom()+jets[1].mom()).mass();
0879           double pTHjj = (jets[0].momentum() + jets[1].momentum() + higgs.momentum()).pt();
0880           //VBF topology
0881           if (vbfTopo)
0882             return Category(GG2H_GE2J_MJJ_350_700_PTH_0_200_PTHJJ_0_25 + vbfTopo - 1);
0883           //Njets >= 2jets without VBF topology (mjj<350)
0884           if (pTHjj < 25)
0885             return Category(GG2H_GE2J_MJJ_0_350_PTH_0_60_PTHJJ_0_25 + getBin(higgs.pt(), {0, 60, 120, 200}));
0886           else
0887             return Category(GG2H_GE2J_MJJ_0_350_PTH_0_60_PTHJJ_GT25 + getBin(higgs.pt(), {0, 60, 120, 200}));
0888         }
0889       }
0890 
0891       // 2. Electroweak qq->Hqq Stage 1.2 categories
0892       else if (prodMode == HTXS::VBF || (isVH(prodMode) && quarkDecay(V))) {
0893         if (std::abs(higgs.rapidity()) > 2.5)
0894           return QQ2HQQ_FWDH;
0895         int Njets = jets.size();
0896         if (Njets == 0)
0897           return QQ2HQQ_0J;
0898         else if (Njets == 1)
0899           return QQ2HQQ_1J;
0900         else if (Njets >= 2) {
0901           double mjj = (jets[0].mom() + jets[1].mom()).mass();
0902           double pTHjj = (jets[0].momentum() + jets[1].momentum() + higgs.momentum()).pt();
0903           if (mjj < 350) {
0904             if (pTHjj < 25)
0905               return Category(QQ2HQQ_GE2J_MJJ_0_60_PTHJJ_0_25 + getBin(mjj, {0, 60, 120, 350}));
0906             else
0907               return Category(QQ2HQQ_GE2J_MJJ_0_60_PTHJJ_GT25 + getBin(mjj, {0, 60, 120, 350}));
0908           } else {  //mjj>350 GeV
0909             if (higgs.pt() < 200) {
0910               return Category(QQ2HQQ_GE2J_MJJ_350_700_PTH_0_200_PTHJJ_0_25 + vbfTopo - 1);
0911             } else {
0912               return Category(QQ2HQQ_GE2J_MJJ_350_700_PTH_GT200_PTHJJ_0_25 + vbfTopo - 1);
0913             }
0914           }
0915         }
0916       }
0917       // 3. WH->Hlv categories
0918       else if (prodMode == HTXS::WH) {
0919         if (fwdHiggs)
0920           return QQ2HLNU_FWDH;
0921         int Njets = jets.size();
0922         if (Njets == 0)
0923           return Category(QQ2HLNU_PTV_0_75_0J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0924         if (Njets == 1)
0925           return Category(QQ2HLNU_PTV_0_75_1J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0926         return Category(QQ2HLNU_PTV_0_75_GE2J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0927       }
0928       // 4. qq->ZH->llH categories
0929       else if (prodMode == HTXS::QQ2ZH) {
0930         if (fwdHiggs)
0931           return QQ2HLL_FWDH;
0932         int Njets = jets.size();
0933         if (Njets == 0)
0934           return Category(QQ2HLL_PTV_0_75_0J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0935         if (Njets == 1)
0936           return Category(QQ2HLL_PTV_0_75_1J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0937         return Category(QQ2HLL_PTV_0_75_GE2J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0938       }
0939       // 5. gg->ZH->llH categories
0940       else if (prodMode == HTXS::GG2ZH) {
0941         if (fwdHiggs)
0942           return GG2HLL_FWDH;
0943         int Njets = jets.size();
0944         if (Njets == 0)
0945           return Category(GG2HLL_PTV_0_75_0J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0946         if (Njets == 1)
0947           return Category(GG2HLL_PTV_0_75_1J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0948         return Category(GG2HLL_PTV_0_75_GE2J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0949       }
0950       // 6.ttH,bbH,tH categories
0951       else if (prodMode == HTXS::TTH) {
0952         if (fwdHiggs)
0953           return TTH_FWDH;
0954         else
0955           return Category(TTH_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200, 300, 450}));
0956       } else if (prodMode == HTXS::BBH)
0957         return Category(BBH_FWDH + ctrlHiggs);
0958       else if (prodMode == HTXS::TH)
0959         return Category(TH_FWDH + ctrlHiggs);
0960       return UNKNOWN;
0961     }
0962 
0963     /// @}
0964 
0965     /// @name Default Rivet analysis methods and steering methods
0966     /// @{
0967 
0968     /// @brief Sets the Higgs production mode
0969     void setHiggsProdMode(HTXS::HiggsProdMode prodMode) { m_HiggsProdMode = prodMode; }
0970 
0971     /// @brief default Rivet Analysis::init method
0972     /// Booking of histograms, initializing Rivet projection
0973     /// Extracts Higgs production mode from shell variable if not set manually using setHiggsProdMode
0974     void init() override {
0975       printf("==============================================================\n");
0976       printf("========     HiggsTemplateCrossSections Initialization     =========\n");
0977       printf("==============================================================\n");
0978       // check that the production mode has been set
0979       // if running in standalone Rivet the production mode is set through an env variable
0980       if (m_HiggsProdMode == HTXS::UNKNOWN) {
0981         char *pm_env = std::getenv("HIGGSPRODMODE");
0982         string pm(pm_env == nullptr ? "" : pm_env);
0983         if (pm == "GGF")
0984           m_HiggsProdMode = HTXS::GGF;
0985         else if (pm == "VBF")
0986           m_HiggsProdMode = HTXS::VBF;
0987         else if (pm == "WH")
0988           m_HiggsProdMode = HTXS::WH;
0989         else if (pm == "ZH")
0990           m_HiggsProdMode = HTXS::QQ2ZH;
0991         else if (pm == "QQ2ZH")
0992           m_HiggsProdMode = HTXS::QQ2ZH;
0993         else if (pm == "GG2ZH")
0994           m_HiggsProdMode = HTXS::GG2ZH;
0995         else if (pm == "TTH")
0996           m_HiggsProdMode = HTXS::TTH;
0997         else if (pm == "BBH")
0998           m_HiggsProdMode = HTXS::BBH;
0999         else if (pm == "TH")
1000           m_HiggsProdMode = HTXS::TH;
1001         else {
1002           MSG_WARNING("No HIGGSPRODMODE shell variable found. Needed when running Rivet stand-alone.");
1003         }
1004       }
1005 
1006       // Projections for final state particles
1007       const FinalState FS;
1008       declare(FS, "FS");
1009 
1010       // initialize the histograms with for each of the stages
1011       initializeHistos();
1012       m_sumw = 0.0;
1013       m_sumevents = 0;
1014       printf("==============================================================\n");
1015       printf("========             Higgs prod mode %d              =========\n", m_HiggsProdMode);
1016       printf("========          Sucessful Initialization           =========\n");
1017       printf("==============================================================\n");
1018     }
1019 
1020     // Perform the per-event analysis
1021     void analyze(const Event &event) override {
1022       // get the classification
1023       HiggsClassification cat = classifyEvent(event, m_HiggsProdMode);
1024 
1025       // Fill histograms: categorization --> linerize the categories
1026       const double weight = 1.0;
1027       m_sumw += weight;
1028 
1029       int F = cat.stage0_cat % 10, P = cat.stage1_cat_pTjet30GeV / 100;
1030       hist_stage0->fill(cat.stage0_cat / 10 * 2 + F, weight);
1031 
1032       // Stage 1 enum offsets for each production mode: GGF=12, VBF=6, WH= 5, QQ2ZH=5, GG2ZH=4, TTH=2, BBH=2, TH=2
1033       vector<int> offset({0, 1, 13, 19, 24, 29, 33, 35, 37, 39});
1034       int off = offset[P];
1035       hist_stage1_pTjet25->fill(cat.stage1_cat_pTjet25GeV % 100 + off, weight);
1036       hist_stage1_pTjet30->fill(cat.stage1_cat_pTjet30GeV % 100 + off, weight);
1037 
1038       // 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
1039       static vector<int> offset1_2({0, 1, 18, 29, 35, 41, 47, 53, 55, 57});
1040       int off1_2 = offset1_2[P];
1041       // 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
1042       static vector<int> offset1_2f({0, 1, 29, 54, 70, 86, 102, 109, 111, 113});
1043       int off1_2f = offset1_2f[P];
1044       hist_stage1_2_pTjet25->fill(cat.stage1_2_cat_pTjet25GeV % 100 + off1_2, weight);
1045       hist_stage1_2_pTjet30->fill(cat.stage1_2_cat_pTjet30GeV % 100 + off1_2, weight);
1046       hist_stage1_2_fine_pTjet25->fill(cat.stage1_2_fine_cat_pTjet25GeV % 100 + off1_2f, weight);
1047       hist_stage1_2_fine_pTjet30->fill(cat.stage1_2_fine_cat_pTjet30GeV % 100 + off1_2f, weight);
1048 
1049       // Fill histograms: variables used in the categorization
1050       hist_pT_Higgs->fill(cat.higgs.pT(), weight);
1051       hist_y_Higgs->fill(cat.higgs.rapidity(), weight);
1052       hist_pT_V->fill(cat.V.pT(), weight);
1053 
1054       hist_Njets25->fill(cat.jets25.size(), weight);
1055       hist_Njets30->fill(cat.jets30.size(), weight);
1056 
1057       hist_isZ2vv->fill(cat.isZ2vvDecay, weight);
1058 
1059       // Jet variables. Use jet collection with pT threshold at 30 GeV
1060       if (!cat.jets30.empty())
1061         hist_pT_jet1->fill(cat.jets30[0].pt(), weight);
1062       if (cat.jets30.size() >= 2) {
1063         const FourMomentum &j1 = cat.jets30[0].momentum(), &j2 = cat.jets30[1].momentum();
1064         hist_deltay_jj->fill(std::abs(j1.rapidity() - j2.rapidity()), weight);
1065         hist_dijet_mass->fill((j1 + j2).mass(), weight);
1066         hist_pT_Hjj->fill((j1 + j2 + cat.higgs.momentum()).pt(), weight);
1067       }
1068     }
1069 
1070     void printClassificationSummary() {
1071       MSG_INFO(" ====================================================== ");
1072       MSG_INFO("      Higgs Template X-Sec Categorization Tool          ");
1073       MSG_INFO("                Status Code Summary                     ");
1074       MSG_INFO(" ====================================================== ");
1075       bool allSuccess = (m_sumevents == m_errorCount[HTXS::SUCCESS]);
1076       if (allSuccess)
1077         MSG_INFO("     >>>> All " << m_errorCount[HTXS::SUCCESS] << " events successfully categorized!");
1078       else {
1079         MSG_INFO("     >>>> " << m_errorCount[HTXS::SUCCESS] << " events successfully categorized");
1080         MSG_INFO("     >>>> --> the following errors occured:");
1081         MSG_INFO("     >>>> " << m_errorCount[HTXS::PRODMODE_DEFINED] << " had an undefined Higgs production mode.");
1082         MSG_INFO("     >>>> " << m_errorCount[HTXS::MOMENTUM_CONSERVATION] << " failed momentum conservation.");
1083         MSG_INFO("     >>>> " << m_errorCount[HTXS::HIGGS_IDENTIFICATION]
1084                               << " failed to identify a valid Higgs boson.");
1085         MSG_INFO("     >>>> " << m_errorCount[HTXS::HS_VTX_IDENTIFICATION]
1086                               << " failed to identify the hard scatter vertex.");
1087         MSG_INFO("     >>>> " << m_errorCount[HTXS::VH_IDENTIFICATION] << " VH: to identify a valid V-boson.");
1088         MSG_INFO("     >>>> " << m_errorCount[HTXS::TOP_W_IDENTIFICATION]
1089                               << " failed to identify valid Ws from top decay.");
1090       }
1091       MSG_INFO(" ====================================================== ");
1092       MSG_INFO(" ====================================================== ");
1093     }
1094 
1095     void finalize() override {
1096       printClassificationSummary();
1097       double sf = m_sumw > 0 ? 1.0 / m_sumw : 1.0;
1098       scale({hist_stage0,
1099              hist_stage1_pTjet25,
1100              hist_stage1_pTjet30,
1101              hist_stage1_2_pTjet25,
1102              hist_stage1_2_pTjet30,
1103              hist_stage1_2_fine_pTjet25,
1104              hist_stage1_2_fine_pTjet30,
1105              hist_Njets25,
1106              hist_Njets30,
1107              hist_pT_Higgs,
1108              hist_y_Higgs,
1109              hist_pT_V,
1110              hist_pT_jet1,
1111              hist_deltay_jj,
1112              hist_dijet_mass,
1113              hist_pT_Hjj},
1114             sf);
1115     }
1116 
1117     /*
1118      *  initialize histograms
1119      */
1120 
1121     void initializeHistos() {
1122       book(hist_stage0, "HTXS_stage0", 20, 0, 20);
1123       book(hist_stage1_pTjet25, "HTXS_stage1_pTjet25", 40, 0, 40);
1124       book(hist_stage1_pTjet30, "HTXS_stage1_pTjet30", 40, 0, 40);
1125       book(hist_stage1_2_pTjet25, "HTXS_stage1_2_pTjet25", 57, 0, 57);
1126       book(hist_stage1_2_pTjet30, "HTXS_stage1_2_pTjet30", 57, 0, 57);
1127       book(hist_stage1_2_fine_pTjet25, "HTXS_stage1_2_fine_pTjet25", 113, 0, 113);
1128       book(hist_stage1_2_fine_pTjet30, "HTXS_stage1_2_fine_pTjet30", 113, 0, 113);
1129       book(hist_pT_Higgs, "pT_Higgs", 80, 0, 400);
1130       book(hist_y_Higgs, "y_Higgs", 80, -4, 4);
1131       book(hist_pT_V, "pT_V", 80, 0, 400);
1132       book(hist_pT_jet1, "pT_jet1", 80, 0, 400);
1133       book(hist_deltay_jj, "deltay_jj", 50, 0, 10);
1134       book(hist_dijet_mass, "m_jj", 50, 0, 2000);
1135       book(hist_pT_Hjj, "pT_Hjj", 50, 0, 250);
1136       book(hist_Njets25, "Njets25", 10, 0, 10);
1137       book(hist_Njets30, "Njets30", 10, 0, 10);
1138       book(hist_isZ2vv, "isZ2vv", 2, 0, 2);
1139     }
1140     /// @}
1141 
1142     /*
1143      *    initialize private members used in the classification procedure
1144      */
1145 
1146   private:
1147     double m_sumw;
1148     size_t m_sumevents;
1149     HTXS::HiggsProdMode m_HiggsProdMode;
1150     std::map<HTXS::ErrorCode, size_t> m_errorCount;
1151     Histo1DPtr hist_stage0;
1152     Histo1DPtr hist_stage1_pTjet25, hist_stage1_pTjet30;
1153     Histo1DPtr hist_stage1_2_pTjet25, hist_stage1_2_pTjet30;
1154     Histo1DPtr hist_stage1_2_fine_pTjet25, hist_stage1_2_fine_pTjet30;
1155     Histo1DPtr hist_pT_Higgs, hist_y_Higgs;
1156     Histo1DPtr hist_pT_V, hist_pT_jet1;
1157     Histo1DPtr hist_deltay_jj, hist_dijet_mass, hist_pT_Hjj;
1158     Histo1DPtr hist_Njets25, hist_Njets30;
1159     Histo1DPtr hist_isZ2vv;
1160   };
1161 
1162   // the PLUGIN only needs to be decleared when running standalone Rivet
1163   // and causes compilation / linking issues if included in Athena / RootCore
1164   //check for Rivet environment variable RIVET_ANALYSIS_PATH
1165 #ifdef RIVET_ANALYSIS_PATH
1166   // The hook for the plugin system
1167   DECLARE_RIVET_PLUGIN(HiggsTemplateCrossSections);
1168 #endif
1169 
1170 }  // namespace Rivet
1171 
1172 #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!