Rivet::HiggsTemplateCrossSections Class Reference

Rivet routine for classifying MC events according to the Higgs template cross section categories. More...

Inheritance diagram for Rivet::HiggsTemplateCrossSections:

Public Member Functions
HiggsClassification	classifyEvent (const Event &event, const HTXS::HiggsProdMode prodMode)
	Main classificaion method. More...
	HiggsTemplateCrossSections ()
Utility methods
Methods to identify the Higgs boson and associated vector boson and to build jets
Particle	getLastInstance (Particle ptcl)
	follow a "propagating" particle and return its last instance More...
bool	originateFrom (const Particle &p, const Particles &ptcls)
	Whether particle p originate from any of the ptcls. More...
bool	originateFrom (const Particle &p, const Particle &p2)
	Whether particle p originates from p2. More...
bool	hasChild (const GenParticle *ptcl, int pdgID)
	Checks whether the input particle has a child with a given PDGID. More...
bool	hasParent (const GenParticle *ptcl, int pdgID)
	Checks whether the input particle has a parent with a given PDGID. More...
bool	quarkDecay (const Particle &p)
	Return true is particle decays to quarks. More...
HiggsClassification	error (HiggsClassification &cat, HTXS::ErrorCode err, std::string msg="", int NmaxWarnings=20)
	Returns the classification object with the error code set. Prints an warning message, and keeps track of number of errors. More...
Categorization methods
Methods to assign the truth category based on the identified Higgs boson and associated vector bosons and/or reconstructed jets
int	getBin (double x, std::vector< double > bins)
	Return bin index of x given the provided bin edges. 0=first bin, -1=underflow bin. More...
int	vbfTopology (const Jets &jets, const Particle &higgs)
	VBF topolog selection 0 = fail loose selction: m_jj > 400 GeV and Dy_jj > 2.8 1 pass loose, but fail additional cut pT(Hjj)<25. 2 pass tight selection. More...
bool	isVH (HTXS::HiggsProdMode p)
	Whether the Higgs is produced in association with a vector boson (VH) More...
HTXS::Stage0::Category	getStage0Category (const HTXS::HiggsProdMode prodMode, const Particle &higgs, const Particle &V)
	Stage-0 HTXS categorization. More...
HTXS::Stage1::Category	getStage1Category (const HTXS::HiggsProdMode prodMode, const Particle &higgs, const Jets &jets, const Particle &V)
	Stage-1 categorization. More...
Default Rivet analysis methods and steering methods
void	setHiggsProdMode (HTXS::HiggsProdMode prodMode)
	Sets the Higgs production mode. More...
void	init () override
	default Rivet Analysis::init method Booking of histograms, initializing Rivet projection Extracts Higgs production mode from shell variable if not set manually using setHiggsProdMode More...
void	analyze (const Event &event) override
void	printClassificationSummary ()
void	finalize () override
void	initializeHistos ()

Private Attributes
Histo1DPtr	hist_deltay_jj
Histo1DPtr	hist_dijet_mass
Histo1DPtr	hist_Njets25
Histo1DPtr	hist_Njets30
Histo1DPtr	hist_pT_Higgs
Histo1DPtr	hist_pT_Hjj
Histo1DPtr	hist_pT_jet1
Histo1DPtr	hist_pT_V
Histo1DPtr	hist_stage0
Histo1DPtr	hist_stage1_pTjet25
Histo1DPtr	hist_stage1_pTjet30
Histo1DPtr	hist_y_Higgs
std::map< HTXS::ErrorCode, size_t >	m_errorCount
HTXS::HiggsProdMode	m_HiggsProdMode
double	m_sumw

Detailed Description

Rivet routine for classifying MC events according to the Higgs template cross section categories.

Author

Jim Lacey (DESY) <james.nosp@m..lac.nosp@m.ey@ce.nosp@m.rn.c.nosp@m.h,jlace.nosp@m.y@de.nosp@m.sy.de>

Dag Gillberg (Carleton University) dag.g.nosp@m.illb.nosp@m.erg@c.nosp@m.ern..nosp@m.ch

Definition at line 19 of file HiggsTemplateCrossSections.cc.

Constructor & Destructor Documentation

Rivet::HiggsTemplateCrossSections::HiggsTemplateCrossSections ( )

inline

Definition at line 22 of file HiggsTemplateCrossSections.cc.

      : Analysis("HiggsTemplateCrossSections"),
    m_HiggsProdMode(HTXS::UNKNOWN) {}

Member Function Documentation

void Rivet::HiggsTemplateCrossSections::analyze ( const Event &  event )

inlineoverride

Definition at line 435 of file HiggsTemplateCrossSections.cc.

References funct::abs(), eostools::cat(), classifyEvent(), F(), hist_deltay_jj, hist_dijet_mass, hist_Njets25, hist_Njets30, hist_pT_Higgs, hist_pT_Hjj, hist_pT_jet1, hist_pT_V, hist_stage0, hist_stage1_pTjet25, hist_stage1_pTjet30, hist_y_Higgs, m_HiggsProdMode, m_sumw, ResonanceBuilder::mass, PFRecoTauDiscriminationByIsolation_cfi::offset, EnergyCorrector::pt, and mps_merge::weight.

                                              {

      // get the classification
      HiggsClassification cat = classifyEvent(event,m_HiggsProdMode);

      // Fill histograms: categorization --> linerize the categories
      const double weight = event.weight();
      m_sumw += weight;

      int F=cat.stage0_cat%10, P=cat.stage1_cat_pTjet30GeV/100;
      hist_stage0->fill( cat.stage0_cat/10*2+F, weight );
      
      // Stage 1 enum offsets for each production mode: GGF=12, VBF=6, WH= 5, QQ2ZH=5, GG2ZH=4, TTH=2, BBH=2, TH=2
      vector<int> offset({0,1,13,19,24,29,33,35,37,39});
      int off = offset[P];
      hist_stage1_pTjet25->fill(cat.stage1_cat_pTjet25GeV%100 + off, weight);
      hist_stage1_pTjet30->fill(cat.stage1_cat_pTjet30GeV%100 + off, weight);

      // Fill histograms: variables used in the categorization
      hist_pT_Higgs->fill(cat.higgs.pT(),weight);
      hist_y_Higgs->fill(cat.higgs.rapidity(),weight);
      hist_pT_V->fill(cat.V.pT(),weight);

      hist_Njets25->fill(cat.jets25.size(),weight);
      hist_Njets30->fill(cat.jets30.size(),weight);

      // Jet variables. Use jet collection with pT threshold at 30 GeV
      if (!cat.jets30.empty()) hist_pT_jet1->fill(cat.jets30[0].pt(),weight);
      if (cat.jets30.size()>=2) {
    const FourMomentum &j1 = cat.jets30[0].momentum(), &j2 = cat.jets30[1].momentum();
    hist_deltay_jj->fill(std::abs(j1.rapidity()-j2.rapidity()),weight);
    hist_dijet_mass->fill((j1+j2).mass(),weight);
    hist_pT_Hjj->fill((j1+j2+cat.higgs.momentum()).pt(),weight);
      }
    }

HiggsClassification Rivet::HiggsTemplateCrossSections::classifyEvent ( const Event &  event, const HTXS::HiggsProdMode  prodMode  )

inline

Main classificaion method.

Definition at line 99 of file HiggsTemplateCrossSections.cc.

References eostools::cat(), class-composition::child, class-composition::children, error(), event(), getLastInstance(), getStage0Category(), getStage1Category(), HTXS::GG2ZH, hasChild(), hasParent(), HTXS::HIGGS_DECAY_IDENTIFICATION, HTXS::HIGGS_IDENTIFICATION, HTXS::HS_VTX_IDENTIFICATION, isVH(), fwrapper::jets, m_errorCount, m_HiggsProdMode, originateFrom(), AlCaHLTBitMon_ParallelJobs::p, HadronAndPartonSelector_cfi::particles, HTXS::PRODMODE_DEFINED, PVValHelper::pT, HTXS::QQ2ZH, quarkDecay(), HTXS::SUCCESS, HTXS::TH, HTXS::TOP_W_IDENTIFICATION, HTXS::TTH, HTXS::UNDEFINED, HTXS::UNKNOWN, HTXS::Stage0::UNKNOWN, HTXS::Stage1::UNKNOWN, HTXS::VH_DECAY_IDENTIFICATION, HTXS::VH_IDENTIFICATION, and HTXS::WH.

Referenced by analyze().

                                                                                             {

      if (m_HiggsProdMode==HTXS::UNKNOWN) m_HiggsProdMode = prodMode; 
      
      // the classification object
      HiggsClassification cat;
      cat.prodMode = prodMode;
      cat.errorCode = HTXS::UNDEFINED;
      cat.stage0_cat = HTXS::Stage0::UNKNOWN;
      cat.stage1_cat_pTjet25GeV = HTXS::Stage1::UNKNOWN;
      cat.stage1_cat_pTjet30GeV = HTXS::Stage1::UNKNOWN;

      if (prodMode == HTXS::UNKNOWN) 
    return error(cat,HTXS::PRODMODE_DEFINED,
             "Unkown Higgs production mechanism. Cannot classify event."                                                                       
             " Classification for all events will most likely fail.");

      /*****
       * Step 1. 
       *  Idenfify the Higgs boson and the hard scatter vertex
       *  There should be only one of each.
       */

      GenVertex *HSvtx = event.genEvent()->signal_process_vertex();
      int Nhiggs=0;
      for ( const GenParticle *ptcl : Rivet::particles(event.genEvent()) ) {

        // a) Reject all non-Higgs particles
        if ( !PID::isHiggs(ptcl->pdg_id()) ) continue;
        // b) select only the final Higgs boson copy, prior to decay
        if ( ptcl->end_vertex() && !hasChild(ptcl,PID::HIGGS) ) {
          cat.higgs = Particle(ptcl); ++Nhiggs;
        }
        // c) if HepMC::signal_proces_vertex is missing
        //    set hard-scatter vertex based on first Higgs boson
        if ( HSvtx==nullptr && ptcl->production_vertex() && !hasParent(ptcl,PID::HIGGS) )
      HSvtx = ptcl->production_vertex();
      }

      // Make sure things are in order so far
      if (Nhiggs!=1) 
    return error(cat,HTXS::HIGGS_IDENTIFICATION,
             "Current event has "+std::to_string(Nhiggs)+" Higgs bosons. There must be only one.");
      if (cat.higgs.children().size()<2) 
    return error(cat,HTXS::HIGGS_DECAY_IDENTIFICATION,
             "Could not identify Higgs boson decay products.");

      if (HSvtx == nullptr) 
    return error(cat,HTXS::HS_VTX_IDENTIFICATION,"Cannot find hard-scatter vertex of current event.");

      /*****
       * Step 2. 
       *   Identify associated vector bosons
       */
 
      // Find associated vector bosons
      bool is_uncatdV = false;
      Particles uncatV_decays;
      FourMomentum uncatV_p4(0,0,0,0);
      FourVector uncatV_v4(0,0,0,0);
      int nWs=0, nZs=0, nTop=0;
      if ( isVH(prodMode) ) {
    for (auto ptcl:particles(HSvtx,HepMC::children)) {
      if (PID::isW(ptcl->pdg_id())) { ++nWs; cat.V=Particle(ptcl); }
      if (PID::isZ(ptcl->pdg_id())) { ++nZs; cat.V=Particle(ptcl); }
    }
    if(nWs+nZs>0) cat.V = getLastInstance(cat.V);
    else {
      for (auto ptcl:particles(HSvtx,HepMC::children)) {
        if (!PID::isHiggs(ptcl->pdg_id())) {
          uncatV_decays += Particle(ptcl);
          uncatV_p4 += Particle(ptcl).momentum();
          // uncatV_v4 += Particle(ptcl).origin();
        }
      }
      // is_uncatdV = true; cat.V = Particle(24,uncatV_p4,uncatV_v4);
      is_uncatdV = true; cat.V = Particle(24,uncatV_p4);
    }
      }
      
      if ( !is_uncatdV ){

    if ( isVH(prodMode) && !cat.V.genParticle()->end_vertex() )
      return error(cat,HTXS::VH_DECAY_IDENTIFICATION,"Vector boson does not decay!");

    if ( isVH(prodMode) && cat.V.children().size()<2 )
      return error(cat,HTXS::VH_DECAY_IDENTIFICATION,"Vector boson does not decay!");
    
    if ( ( prodMode==HTXS::WH && (nZs>0||nWs!=1) ) ||
         ( (prodMode==HTXS::QQ2ZH||prodMode==HTXS::GG2ZH) && (nZs!=1||nWs>0) ) ) 
      return error(cat,HTXS::VH_IDENTIFICATION,"Found "+std::to_string(nWs)+" W-bosons and "+
               std::to_string(nZs)+" Z-bosons. Inconsitent with VH expectation.");
      }
      
      // Find and store the W-bosons from ttH->WbWbH
      Particles Ws;
      if ( prodMode==HTXS::TTH || prodMode==HTXS::TH ){
    // loop over particles produced in hard-scatter vertex
        for ( auto ptcl : particles(HSvtx,HepMC::children) ) {
          if ( !PID::isTop(ptcl->pdg_id()) ) continue;
      ++nTop;
      Particle top = getLastInstance(Particle(ptcl));
      if ( top.genParticle()->end_vertex() ) 
        for (auto child:top.children())
          if ( PID::isW(child.pdgId()) ) Ws += child;
    }
      }

      // Make sure result make sense
      if ( (prodMode==HTXS::TTH && Ws.size()<2) || (prodMode==HTXS::TH && Ws.empty() ) )
    return error(cat,HTXS::TOP_W_IDENTIFICATION,"Failed to identify W-boson(s) from t-decay!");

      /*****
       * Step 3.
       *   Build jets
       *   Make sure all stable particles are present
       */

      // Create a list of the vector bosons that decay leptonically
      // Either the vector boson produced in association with the Higgs boson,
      // or the ones produced from decays of top quarks produced with the Higgs
      Particles leptonicVs;
      if ( !is_uncatdV ){
    if ( isVH(prodMode) && !quarkDecay(cat.V) ) leptonicVs += cat.V;
      }else leptonicVs = uncatV_decays;
      for ( auto W:Ws ) if ( W.genParticle()->end_vertex() && !quarkDecay(W) ) leptonicVs += W;

      // Obtain all stable, final-state particles
      const ParticleVector FS = applyProjection<FinalState>(event, "FS").particles();
      Particles hadrons;
      FourMomentum sum(0,0,0,0), vSum(0,0,0,0), hSum(0,0,0,0);
      for ( const Particle &p : FS ) {
        // Add up the four momenta of all stable particles as a cross check
    sum += p.momentum();
    // ignore particles from the Higgs boson
    if ( originateFrom(p,cat.higgs) ) { hSum += p.momentum(); continue; }
    // Cross-check the V decay products for VH
    if ( isVH(prodMode) && !is_uncatdV && originateFrom(p,Ws) ) vSum += p.momentum();
    // ignore final state particles from leptonic V decays
    if ( !leptonicVs.empty() && originateFrom(p,leptonicVs) ) continue;
    // All particles reaching here are considered hadrons and will be used to build jets
    hadrons += p;
      }
      
      cat.p4decay_higgs = hSum;
      cat.p4decay_V = vSum;

      FinalState fps_temp;
      FastJets jets(fps_temp, FastJets::ANTIKT, 0.4 );
      jets.calc(hadrons);

      cat.jets25 = jets.jetsByPt( Cuts::pT > 25.0 );
      cat.jets30 = jets.jetsByPt( Cuts::pT > 30.0 );
 
      // check that four mometum sum of all stable particles satisfies momentum consevation
/*
      if ( sum.pt()>0.1 )
    return error(cat,HTXS::MOMENTUM_CONSERVATION,"Four vector sum does not amount to pT=0, m=E=sqrt(s), but pT="+
             std::to_string(sum.pt())+" GeV and m = "+std::to_string(sum.mass())+" GeV");
*/      
      // check if V-boson was not included in the event record but decay particles were
      // EFT contact interaction: return UNKNOWN for category but set all event/particle kinematics
      if(is_uncatdV) 
    return error(cat,HTXS::VH_IDENTIFICATION,"Failed to identify associated V-boson!");
       
      /*****
       * Step 4.
       *   Classify and save output
       */
      
      // Apply the categorization categorization
      cat.stage0_cat = getStage0Category(prodMode,cat.higgs,cat.V);
      cat.stage1_cat_pTjet25GeV = getStage1Category(prodMode,cat.higgs,cat.jets25,cat.V);
      cat.stage1_cat_pTjet30GeV = getStage1Category(prodMode,cat.higgs,cat.jets30,cat.V);      
      cat.errorCode = HTXS::SUCCESS; ++m_errorCount[HTXS::SUCCESS];

      return cat;
    }    

HiggsClassification Rivet::HiggsTemplateCrossSections::error ( HiggsClassification &  cat, HTXS::ErrorCode  err, std::string  msg = "", int  NmaxWarnings = 20  )

inline

Returns the classification object with the error code set. Prints an warning message, and keeps track of number of errors.

Definition at line 83 of file HiggsTemplateCrossSections.cc.

References eostools::cat(), m_errorCount, and mps_check::msg.

Referenced by argparse.ArgumentParser::_get_option_tuples(), python.rootplot.argparse.ArgumentParser::_get_option_tuples(), argparse.ArgumentParser::_parse_known_args(), python.rootplot.argparse.ArgumentParser::_parse_known_args(), argparse.ArgumentParser::_parse_optional(), python.rootplot.argparse.ArgumentParser::_parse_optional(), argparse.ArgumentParser::_read_args_from_files(), python.rootplot.argparse.ArgumentParser::_read_args_from_files(), argparse.ArgumentParser::add_subparsers(), python.rootplot.argparse.ArgumentParser::add_subparsers(), Page1Parser.Page1Parser::check_for_whole_start_tag(), classifyEvent(), argparse.ArgumentParser::parse_args(), python.rootplot.argparse.ArgumentParser::parse_args(), argparse.ArgumentParser::parse_known_args(), and python.rootplot.argparse.ArgumentParser::parse_known_args().

                                                         {
      // Set the error, and keep statistics
      cat.errorCode = err;
      ++m_errorCount[err];

      // Print warning message to the screen/log
      static int Nwarnings = 0;
      if ( msg!="" && ++Nwarnings < NmaxWarnings )
      MSG_WARNING(msg);

      return cat;
    }

void Rivet::HiggsTemplateCrossSections::finalize ( void  )

inlineoverride

Definition at line 493 of file HiggsTemplateCrossSections.cc.

References create_public_lumi_plots::hist, hist_deltay_jj, hist_dijet_mass, hist_Njets25, hist_Njets30, hist_pT_Higgs, hist_pT_Hjj, hist_pT_jet1, hist_pT_V, hist_stage0, hist_stage1_pTjet25, hist_stage1_pTjet30, hist_y_Higgs, m_sumw, printClassificationSummary(), and Scenarios_cff::scale.

                             {
      printClassificationSummary();
      double sf = m_sumw>0?1.0/m_sumw:1.0;
      for (auto hist:{hist_stage0,hist_stage1_pTjet25,hist_stage1_pTjet30,hist_Njets25,hist_Njets30,
        hist_pT_Higgs,hist_y_Higgs,hist_pT_V,hist_pT_jet1,hist_deltay_jj,hist_dijet_mass,hist_pT_Hjj})
    scale(hist, sf);
    }

int Rivet::HiggsTemplateCrossSections::getBin ( double  x, std::vector< double >  bins  )

inline

Return bin index of x given the provided bin edges. 0=first bin, -1=underflow bin.

Definition at line 285 of file HiggsTemplateCrossSections.cc.

References mps_fire::i.

Referenced by BTagWeightCalculator.BTagWeightCalculator::calcJetWeightImpl(), and getStage1Category().

                                                 {
      for (size_t i=1;i<bins.size();++i)
        if (x<bins[i]) return i-1;
      return bins.size()-1;
    }

Particle Rivet::HiggsTemplateCrossSections::getLastInstance ( Particle  ptcl )

inline

follow a "propagating" particle and return its last instance

Definition at line 34 of file HiggsTemplateCrossSections.cc.

References hasChild().

Referenced by classifyEvent().

                                            {
      if ( ptcl.genParticle()->end_vertex() ) {
    if ( !hasChild(ptcl.genParticle(),ptcl.pdgId()) ) return ptcl;
    else return getLastInstance(ptcl.children()[0]);
      }
      return ptcl;
    }

HTXS::Stage0::Category Rivet::HiggsTemplateCrossSections::getStage0Category ( const HTXS::HiggsProdMode  prodMode, const Particle &  higgs, const Particle &  V  )

inline

Stage-0 HTXS categorization.

Definition at line 305 of file HiggsTemplateCrossSections.cc.

References funct::abs(), HTXS::GG2ZH, HTXS::GGF, HTXS::QQ2ZH, quarkDecay(), HTXS::Stage0::VH2HQQ, HTXS::Stage0::VH2HQQ_FWDH, and HTXS::WH.

Referenced by classifyEvent().

                                            {
      using namespace HTXS::Stage0;
      int ctrlHiggs = std::abs(higgs.rapidity())<2.5;
      // Special cases first, qq→Hqq
      if ( (prodMode==HTXS::WH||prodMode==HTXS::QQ2ZH) && quarkDecay(V) ) {
    return ctrlHiggs ? VH2HQQ : VH2HQQ_FWDH;
      } else if ( prodMode==HTXS::GG2ZH && quarkDecay(V) ) {
    return Category(HTXS::GGF*10 + ctrlHiggs);
      }
      // General case after
      return  Category(prodMode*10 + ctrlHiggs);
    }

HTXS::Stage1::Category Rivet::HiggsTemplateCrossSections::getStage1Category ( const HTXS::HiggsProdMode  prodMode, const Particle &  higgs, const Jets &  jets, const Particle &  V  )

inline

Stage-1 categorization.

Definition at line 321 of file HiggsTemplateCrossSections.cc.

References funct::abs(), HTXS::BBH, HTXS::Stage0::BBH_FWDH, getBin(), HTXS::Stage1::GG2H_0J, HTXS::Stage1::GG2H_1J_PTH_0_60, HTXS::Stage0::GG2H_FWDH, HTXS::Stage1::GG2H_GE2J_PTH_0_60, HTXS::Stage1::GG2H_VBFTOPO_JET3, HTXS::Stage1::GG2H_VBFTOPO_JET3VETO, HTXS::Stage0::GG2HLL_FWDH, HTXS::Stage1::GG2HLL_PTV_0_150, HTXS::Stage1::GG2HLL_PTV_GT150_0J, HTXS::Stage1::GG2HLL_PTV_GT150_GE1J, HTXS::GG2ZH, HTXS::GGF, isVH(), ResonanceBuilder::mass, ECF_cff::Njets, HTXS::Stage0::QQ2HLL_FWDH, HTXS::Stage1::QQ2HLL_PTV_0_150, HTXS::Stage1::QQ2HLL_PTV_150_250_0J, HTXS::Stage1::QQ2HLL_PTV_150_250_GE1J, HTXS::Stage1::QQ2HLL_PTV_GT250, HTXS::Stage0::QQ2HLNU_FWDH, HTXS::Stage1::QQ2HLNU_PTV_0_150, HTXS::Stage1::QQ2HLNU_PTV_150_250_0J, HTXS::Stage1::QQ2HLNU_PTV_150_250_GE1J, HTXS::Stage1::QQ2HLNU_PTV_GT250, HTXS::Stage1::QQ2HQQ_FWDH, HTXS::Stage1::QQ2HQQ_PTJET1_GT200, HTXS::Stage1::QQ2HQQ_REST, HTXS::Stage1::QQ2HQQ_VBFTOPO_JET3, HTXS::Stage1::QQ2HQQ_VBFTOPO_JET3VETO, HTXS::Stage1::QQ2HQQ_VH2JET, HTXS::QQ2ZH, quarkDecay(), HTXS::TH, HTXS::Stage0::TH_FWDH, HTXS::TTH, HTXS::Stage0::TTH_FWDH, pat::UNKNOWN, HTXS::VBF, vbfTopology(), and HTXS::WH.

Referenced by classifyEvent().

                                            {
      using namespace HTXS::Stage1;
      int Njets=jets.size(), ctrlHiggs = std::abs(higgs.rapidity())<2.5, fwdHiggs = !ctrlHiggs;
      double pTj1 = !jets.empty() ? jets[0].momentum().pt() : 0;
      int vbfTopo = vbfTopology(jets,higgs);

      // 1. GGF Stage 1 categories
      //    Following YR4 write-up: XXXXX
      if (prodMode==HTXS::GGF || (prodMode==HTXS::GG2ZH && quarkDecay(V)) ) {
    if (fwdHiggs)        return GG2H_FWDH;
    if (Njets==0)        return GG2H_0J;
    else if (Njets==1)   return Category(GG2H_1J_PTH_0_60+getBin(higgs.pt(),{0,60,120,200}));
    else if (Njets>=2) {
      // events with pT_H>200 get priority over VBF cuts
      if(higgs.pt()<=200){
        if      (vbfTopo==2) return GG2H_VBFTOPO_JET3VETO;
        else if (vbfTopo==1) return GG2H_VBFTOPO_JET3;
      }
      // Njets >= 2jets without VBF topology
      return Category(GG2H_GE2J_PTH_0_60+getBin(higgs.pt(),{0,60,120,200}));
    }
      }
      // 2. Electroweak qq->Hqq Stage 1 categories
      else if (prodMode==HTXS::VBF || ( isVH(prodMode) && quarkDecay(V)) ) {
    if (std::abs(higgs.rapidity())>2.5) return QQ2HQQ_FWDH;
        if (pTj1>200) return QQ2HQQ_PTJET1_GT200;
    if (vbfTopo==2) return QQ2HQQ_VBFTOPO_JET3VETO;
    if (vbfTopo==1) return QQ2HQQ_VBFTOPO_JET3;
    double mjj = jets.size()>1 ? (jets[0].mom()+jets[1].mom()).mass():0;
    if ( 60 < mjj && mjj < 120 ) return QQ2HQQ_VH2JET;
        return QQ2HQQ_REST;
      }
      // 3. WH->Hlv categories
      else if (prodMode==HTXS::WH) {
    if (fwdHiggs) return QQ2HLNU_FWDH;
    else if (V.pt()<150) return QQ2HLNU_PTV_0_150;
        else if (V.pt()>250) return QQ2HLNU_PTV_GT250;
        // 150 < pTV/GeV < 250
        return jets.empty() ? QQ2HLNU_PTV_150_250_0J : QQ2HLNU_PTV_150_250_GE1J;
      }
      // 4. qq->ZH->llH categories
      else if (prodMode==HTXS::QQ2ZH) {
    if (fwdHiggs) return QQ2HLL_FWDH;
        else if (V.pt()<150) return QQ2HLL_PTV_0_150;
        else if (V.pt()>250) return QQ2HLL_PTV_GT250;
        // 150 < pTV/GeV < 250
        return jets.empty() ? QQ2HLL_PTV_150_250_0J : QQ2HLL_PTV_150_250_GE1J;
      }
      // 5. gg->ZH->llH categories
      else if (prodMode==HTXS::GG2ZH ) {
    if (fwdHiggs) return GG2HLL_FWDH;
    if      (V.pt()<150) return GG2HLL_PTV_0_150;
    else if (jets.empty()) return GG2HLL_PTV_GT150_0J;
    return GG2HLL_PTV_GT150_GE1J;
      }
      // 6.ttH,bbH,tH categories
      else if (prodMode==HTXS::TTH) return Category(TTH_FWDH+ctrlHiggs);
      else if (prodMode==HTXS::BBH) return Category(BBH_FWDH+ctrlHiggs);
      else if (prodMode==HTXS::TH ) return Category(TH_FWDH+ctrlHiggs);
      return UNKNOWN;
    }

bool Rivet::HiggsTemplateCrossSections::hasChild ( const GenParticle *  ptcl, int  pdgID  )

inline

Checks whether the input particle has a child with a given PDGID.

Definition at line 61 of file HiggsTemplateCrossSections.cc.

Referenced by classifyEvent(), and getLastInstance().

                                                      {
      for (auto child:Particle(*ptcl).children())
        if (child.pdgId()==pdgID) return true;
      return false;
    }

bool Rivet::HiggsTemplateCrossSections::hasParent ( const GenParticle *  ptcl, int  pdgID  )

inline

Checks whether the input particle has a parent with a given PDGID.

Definition at line 68 of file HiggsTemplateCrossSections.cc.

References class-composition::parent, parents, and HadronAndPartonSelector_cfi::particles.

Referenced by classifyEvent().

                                                       {
      for (auto parent:particles(ptcl->production_vertex(),HepMC::parents))
        if (parent->pdg_id()==pdgID) return true;
      return false;
    }

void Rivet::HiggsTemplateCrossSections::init ( void  )

inlineoverride

default Rivet Analysis::init method Booking of histograms, initializing Rivet projection Extracts Higgs production mode from shell variable if not set manually using setHiggsProdMode

Definition at line 398 of file HiggsTemplateCrossSections.cc.

References HTXS::BBH, HTXS::GG2ZH, HTXS::GGF, initializeHistos(), m_HiggsProdMode, m_sumw, HTXS::QQ2ZH, HTXS::TH, HTXS::TTH, HTXS::UNKNOWN, HTXS::VBF, and HTXS::WH.

                         {
      printf("==============================================================\n");
      printf("========     HiggsTemplateCrossSections Initialization     =========\n");
      printf("==============================================================\n");
      // check that the production mode has been set
      // if running in standalone Rivet the production mode is set through an env variable
      if (m_HiggsProdMode==HTXS::UNKNOWN) {
    char *pm_env = getenv("HIGGSPRODMODE");
    string pm(pm_env==nullptr?"":pm_env);
        if      ( pm == "GGF"   ) m_HiggsProdMode = HTXS::GGF;
        else if ( pm == "VBF"   ) m_HiggsProdMode = HTXS::VBF;
        else if ( pm == "WH"    ) m_HiggsProdMode = HTXS::WH;
        else if ( pm == "ZH"    ) m_HiggsProdMode = HTXS::QQ2ZH;
        else if ( pm == "QQ2ZH" ) m_HiggsProdMode = HTXS::QQ2ZH;
        else if ( pm == "GG2ZH" ) m_HiggsProdMode = HTXS::GG2ZH;
        else if ( pm == "TTH"   ) m_HiggsProdMode = HTXS::TTH;
        else if ( pm == "BBH"   ) m_HiggsProdMode = HTXS::BBH;
        else if ( pm == "TH"    ) m_HiggsProdMode = HTXS::TH;
        else {
          MSG_WARNING("No HIGGSPRODMODE shell variable found. Needed when running Rivet stand-alone.");
        }
      }

      // Projections for final state particles
      const FinalState FS; 
      addProjection(FS,"FS");

      // initialize the histograms with for each of the stages
      initializeHistos();
      m_sumw = 0.0;
      printf("==============================================================\n");
      printf("========             Higgs prod mode %d              =========\n",m_HiggsProdMode);
      printf("========          Sucessful Initialization           =========\n");
      printf("==============================================================\n");
    }

void Rivet::HiggsTemplateCrossSections::initializeHistos ( )

inline

Definition at line 505 of file HiggsTemplateCrossSections.cc.

References hist_deltay_jj, hist_dijet_mass, hist_Njets25, hist_Njets30, hist_pT_Higgs, hist_pT_Hjj, hist_pT_jet1, hist_pT_V, hist_stage0, hist_stage1_pTjet25, hist_stage1_pTjet30, and hist_y_Higgs.

Referenced by init().

                           {
      hist_stage0          = bookHisto1D("HTXS_stage0",20,0,20);
      hist_stage1_pTjet25  = bookHisto1D("HTXS_stage1_pTjet25",40,0,40);
      hist_stage1_pTjet30  = bookHisto1D("HTXS_stage1_pTjet30",40,0,40);
      hist_pT_Higgs    = bookHisto1D("pT_Higgs",80,0,400);
      hist_y_Higgs     = bookHisto1D("y_Higgs",80,-4,4);
      hist_pT_V        = bookHisto1D("pT_V",80,0,400);
      hist_pT_jet1     = bookHisto1D("pT_jet1",80,0,400);
      hist_deltay_jj   = bookHisto1D("deltay_jj",50,0,10);
      hist_dijet_mass  = bookHisto1D("m_jj",50,0,2000);
      hist_pT_Hjj      = bookHisto1D("pT_Hjj",50,0,250);
      hist_Njets25     = bookHisto1D("Njets25",10,0,10);
      hist_Njets30     = bookHisto1D("Njets30",10,0,10);
    }

bool Rivet::HiggsTemplateCrossSections::isVH ( HTXS::HiggsProdMode  p )

inline

Whether the Higgs is produced in association with a vector boson (VH)

Definition at line 302 of file HiggsTemplateCrossSections.cc.

References HTXS::GG2ZH, HTXS::QQ2ZH, and HTXS::WH.

Referenced by classifyEvent(), and getStage1Category().

{ return p==HTXS::WH || p==HTXS::QQ2ZH || p==HTXS::GG2ZH; }

bool Rivet::HiggsTemplateCrossSections::originateFrom ( const Particle &  p, const Particles &  ptcls  )

inline

Whether particle p originate from any of the ptcls.

Definition at line 43 of file HiggsTemplateCrossSections.cc.

References HadronAndPartonSelector_cfi::particles.

Referenced by classifyEvent(), and originateFrom().

                                                                   {
      const GenVertex* prodVtx = p.genParticle()->production_vertex();
      if (prodVtx == nullptr) return false;
      // for each ancestor, check if it matches any of the input particles
      for (const auto & ancestor:particles(prodVtx, HepMC::ancestors)){ 
    for ( auto part:ptcls ) 
      if ( ancestor==part.genParticle() ) return true;
      }
      // if we get here, no ancetor matched any input particle
      return false; 
    }

bool Rivet::HiggsTemplateCrossSections::originateFrom ( const Particle &  p, const Particle &  p2  )

inline

Whether particle p originates from p2.

Definition at line 56 of file HiggsTemplateCrossSections.cc.

References originateFrom().

                                                               { 
      Particles ptcls = {p2}; return originateFrom(p,ptcls);
    }

void Rivet::HiggsTemplateCrossSections::printClassificationSummary ( )

inline

Definition at line 471 of file HiggsTemplateCrossSections.cc.

References HTXS::HIGGS_IDENTIFICATION, HTXS::HS_VTX_IDENTIFICATION, m_errorCount, HTXS::MOMENTUM_CONSERVATION, simpleEdmComparison::numEvents, HTXS::PRODMODE_DEFINED, HTXS::SUCCESS, HTXS::TOP_W_IDENTIFICATION, and HTXS::VH_IDENTIFICATION.

Referenced by finalize().

                                     {
      MSG_INFO (" ====================================================== ");
      MSG_INFO ("      Higgs Template X-Sec Categorization Tool          ");
      MSG_INFO ("                Status Code Summary                     ");
      MSG_INFO (" ====================================================== ");
      bool allSuccess = (numEvents()==m_errorCount[HTXS::SUCCESS]);
      if ( allSuccess ) MSG_INFO ("     >>>> All "<< m_errorCount[HTXS::SUCCESS] <<" events successfully categorized!");
      else{
    MSG_INFO ("     >>>> "<< m_errorCount[HTXS::SUCCESS] <<" events successfully categorized");
    MSG_INFO ("     >>>> --> the following errors occured:");
    MSG_INFO ("     >>>> "<< m_errorCount[HTXS::PRODMODE_DEFINED] <<" had an undefined Higgs production mode.");
    MSG_INFO ("     >>>> "<< m_errorCount[HTXS::MOMENTUM_CONSERVATION] <<" failed momentum conservation.");
    MSG_INFO ("     >>>> "<< m_errorCount[HTXS::HIGGS_IDENTIFICATION] <<" failed to identify a valid Higgs boson.");
    MSG_INFO ("     >>>> "<< m_errorCount[HTXS::HS_VTX_IDENTIFICATION] <<" failed to identify the hard scatter vertex.");
    MSG_INFO ("     >>>> "<< m_errorCount[HTXS::VH_IDENTIFICATION] <<" VH: to identify a valid V-boson.");
    MSG_INFO ("     >>>> "<< m_errorCount[HTXS::TOP_W_IDENTIFICATION] <<" failed to identify valid Ws from top decay.");
      }
      MSG_INFO (" ====================================================== ");
      MSG_INFO (" ====================================================== ");
    }

bool Rivet::HiggsTemplateCrossSections::quarkDecay ( const Particle &  p )

inline

Return true is particle decays to quarks.

Definition at line 75 of file HiggsTemplateCrossSections.cc.

Referenced by classifyEvent(), getStage0Category(), and getStage1Category().

                                       {
      for (auto child:p.children())
        if (PID::isQuark(child.pdgId())) return true;
      return false;
    }

void Rivet::HiggsTemplateCrossSections::setHiggsProdMode ( HTXS::HiggsProdMode  prodMode )

inline

Sets the Higgs production mode.

Definition at line 393 of file HiggsTemplateCrossSections.cc.

References m_HiggsProdMode.

{ m_HiggsProdMode = prodMode; }

int Rivet::HiggsTemplateCrossSections::vbfTopology ( const Jets &  jets, const Particle &  higgs  )

inline

VBF topolog selection 0 = fail loose selction: m_jj > 400 GeV and Dy_jj > 2.8 1 pass loose, but fail additional cut pT(Hjj)<25. 2 pass tight selection.

Definition at line 294 of file HiggsTemplateCrossSections.cc.

References funct::abs(), ResonanceBuilder::mass, and EnergyCorrector::pt.

Referenced by getStage1Category().

                                                             {
      if (jets.size()<2) return 0;
      const FourMomentum &j1=jets[0].momentum(), &j2=jets[1].momentum();
      bool VBFtopo = (j1+j2).mass() > 400.0 && std::abs(j1.rapidity()-j2.rapidity()) > 2.8;
      return VBFtopo ? (j1+j2+higgs.momentum()).pt()<25 ? 2 : 1 : 0;
    }

Member Data Documentation

Histo1DPtr Rivet::HiggsTemplateCrossSections::hist_deltay_jj

private

Definition at line 533 of file HiggsTemplateCrossSections.cc.

Referenced by analyze(), finalize(), and initializeHistos().

Histo1DPtr Rivet::HiggsTemplateCrossSections::hist_dijet_mass

private

Definition at line 533 of file HiggsTemplateCrossSections.cc.

Referenced by analyze(), finalize(), and initializeHistos().

Histo1DPtr Rivet::HiggsTemplateCrossSections::hist_Njets25

private

Definition at line 534 of file HiggsTemplateCrossSections.cc.

Referenced by analyze(), finalize(), and initializeHistos().

Histo1DPtr Rivet::HiggsTemplateCrossSections::hist_Njets30

private

Definition at line 534 of file HiggsTemplateCrossSections.cc.

Referenced by analyze(), finalize(), and initializeHistos().

Histo1DPtr Rivet::HiggsTemplateCrossSections::hist_pT_Higgs

private

Definition at line 531 of file HiggsTemplateCrossSections.cc.

Referenced by analyze(), finalize(), and initializeHistos().

Histo1DPtr Rivet::HiggsTemplateCrossSections::hist_pT_Hjj

private

Definition at line 533 of file HiggsTemplateCrossSections.cc.

Referenced by analyze(), finalize(), and initializeHistos().

Histo1DPtr Rivet::HiggsTemplateCrossSections::hist_pT_jet1

private

Definition at line 532 of file HiggsTemplateCrossSections.cc.

Referenced by analyze(), finalize(), and initializeHistos().

Histo1DPtr Rivet::HiggsTemplateCrossSections::hist_pT_V

private

Definition at line 532 of file HiggsTemplateCrossSections.cc.

Referenced by analyze(), finalize(), and initializeHistos().

Histo1DPtr Rivet::HiggsTemplateCrossSections::hist_stage0

private

Definition at line 529 of file HiggsTemplateCrossSections.cc.

Referenced by analyze(), finalize(), and initializeHistos().

Histo1DPtr Rivet::HiggsTemplateCrossSections::hist_stage1_pTjet25

private

Definition at line 530 of file HiggsTemplateCrossSections.cc.

Referenced by analyze(), finalize(), and initializeHistos().

Histo1DPtr Rivet::HiggsTemplateCrossSections::hist_stage1_pTjet30

private

Definition at line 530 of file HiggsTemplateCrossSections.cc.

Referenced by analyze(), finalize(), and initializeHistos().

Histo1DPtr Rivet::HiggsTemplateCrossSections::hist_y_Higgs

private

Definition at line 531 of file HiggsTemplateCrossSections.cc.

Referenced by analyze(), finalize(), and initializeHistos().

std::map<HTXS::ErrorCode,size_t> Rivet::HiggsTemplateCrossSections::m_errorCount

private

Definition at line 528 of file HiggsTemplateCrossSections.cc.

Referenced by classifyEvent(), error(), and printClassificationSummary().

HTXS::HiggsProdMode Rivet::HiggsTemplateCrossSections::m_HiggsProdMode

private

Definition at line 527 of file HiggsTemplateCrossSections.cc.

Referenced by analyze(), classifyEvent(), init(), and setHiggsProdMode().

double Rivet::HiggsTemplateCrossSections::m_sumw

private

Definition at line 526 of file HiggsTemplateCrossSections.cc.

Referenced by analyze(), finalize(), and init().