Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-02-14 13:07:31 

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