HiggsTemplateCrossSections.cc

Go to the documentation of this file.

#ifndef TRUTHRIVETTOOLS_HIGGSTEMPLATECROSSSECTIONS_CC
#define TRUTHRIVETTOOLS_HIGGSTEMPLATECROSSSECTIONS_CC

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Particle.hh"
#include "Rivet/Projections/FastJets.hh"

// Definition of the StatusCode and Category enums
//#include "HiggsTemplateCrossSections.h"
#include "SimDataFormats/HTXS/interface/HiggsTemplateCrossSections.h" // 

#include <atomic>

namespace Rivet {
  
  class HiggsTemplateCrossSections : public Analysis {
  public:
    // Constructor
    HiggsTemplateCrossSections()
      : Analysis("HiggsTemplateCrossSections"),
    m_HiggsProdMode(HTXS::UNKNOWN) {}

  public:


    Particle getLastInstance(Particle ptcl) {
      if ( ptcl.genParticle()->end_vertex() ) {
    if ( !hasChild(ptcl.genParticle(),ptcl.pdgId()) ) return ptcl;
    else return getLastInstance(ptcl.children()[0]);
      }
      return ptcl;
    }
    
    bool originateFrom(const Particle& p, const Particles& ptcls ) {
      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 originateFrom(const Particle& p, const Particle& p2 ) { 
      Particles ptcls = {p2}; return originateFrom(p,ptcls);
    }
    
    bool hasChild(const GenParticle *ptcl, int pdgID) {
      for (auto child:Particle(*ptcl).children())
        if (child.pdgId()==pdgID) return true;
      return false;
    }
    
    bool hasParent(const GenParticle *ptcl, int pdgID) {
      for (auto parent:particles(ptcl->production_vertex(),HepMC::parents))
        if (parent->pdg_id()==pdgID) return true;
      return false;
    }

    bool quarkDecay(const Particle &p) {
      for (auto child:p.children())
        if (PID::isQuark(child.pdgId())) return true;
      return false;
    }
    
    HiggsClassification error(HiggsClassification &cat, HTXS::ErrorCode err, 
                  std::string msg="", int NmaxWarnings=20) {
      // Set the error, and keep statistics
      cat.errorCode = err;
      ++m_errorCount[err];

      // Print warning message to the screen/log
      static std::atomic<int> Nwarnings{0};
      if ( !msg.empty() && ++Nwarnings < NmaxWarnings )
      MSG_WARNING(msg);

      return cat;
    }

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

      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;
      cat.stage1_1_cat_pTjet25GeV = HTXS::Stage1_1::UNKNOWN;
      cat.stage1_1_cat_pTjet30GeV = HTXS::Stage1_1::UNKNOWN;
      cat.stage1_1_fine_cat_pTjet25GeV = HTXS::Stage1_1_Fine::UNKNOWN;
      cat.stage1_1_fine_cat_pTjet30GeV = HTXS::Stage1_1_Fine::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.stage1_1_cat_pTjet25GeV = getStage1_1_Category(prodMode,cat.higgs,cat.jets25,cat.V);
      cat.stage1_1_cat_pTjet30GeV = getStage1_1_Category(prodMode,cat.higgs,cat.jets30,cat.V);
      cat.stage1_1_fine_cat_pTjet25GeV = getStage1_1_Fine_Category(prodMode,cat.higgs,cat.jets25,cat.V);
      cat.stage1_1_fine_cat_pTjet30GeV = getStage1_1_Fine_Category(prodMode,cat.higgs,cat.jets30,cat.V);

      cat.errorCode = HTXS::SUCCESS; ++m_errorCount[HTXS::SUCCESS];

      return cat;
    }    

    
    int getBin(double x, std::vector<double> bins) {
      for (size_t i=1;i<bins.size();++i)
        if (x<bins[i]) return i-1;
      return bins.size()-1;
    }
    
    int vbfTopology(const Jets &jets, const Particle &higgs) {
      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;
    }

    int vbfTopology_Stage1_1(const Jets &jets, const Particle &higgs) {
      if (jets.size()<2) return 0;
      const FourMomentum &j1=jets[0].momentum(), &j2=jets[1].momentum();
      double mjj = (j1+j2).mass();
      if(mjj>350 && mjj<=700) return (j1+j2+higgs.momentum()).pt()<25 ? 1 : 2;
      else if(mjj>700) return (j1+j2+higgs.momentum()).pt()<25 ? 3 : 4;
      else return 0;
    }

    int vbfTopology_Stage1_1_Fine(const Jets &jets, const Particle &higgs) {
        if (jets.size()<2) return 0;
        const FourMomentum &j1=jets[0].momentum(), &j2=jets[1].momentum();
        double mjj = (j1+j2).mass();
        if(mjj>350 && mjj<=700) return (j1+j2+higgs.momentum()).pt()<25 ? 1 : 2;
        else if(mjj>700 && mjj<=1000) return (j1+j2+higgs.momentum()).pt()<25 ? 3 : 4;
        else if(mjj>1000 && mjj<=1500) return (j1+j2+higgs.momentum()).pt()<25 ? 5 : 6;
        else if(mjj>1500) return (j1+j2+higgs.momentum()).pt()<25 ? 7 : 8;
        else return 0;
    }

    bool isVH(HTXS::HiggsProdMode p) { return p==HTXS::WH || p==HTXS::QQ2ZH || p==HTXS::GG2ZH; }
    
    HTXS::Stage0::Category getStage0Category(const HTXS::HiggsProdMode prodMode,
                         const Particle &higgs,
                         const Particle &V) {
      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 getStage1Category(const HTXS::HiggsProdMode prodMode,
                         const Particle &higgs,
                         const Jets &jets,
                         const Particle &V) {
      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;
    }

    HTXS::Stage1_1::Category getStage1_1_Category(const HTXS::HiggsProdMode prodMode,
                          const Particle &higgs,
                          const Jets &jets,
                          const Particle &V) {
      using namespace HTXS::Stage1_1;
      int Njets=jets.size(), ctrlHiggs = std::abs(higgs.rapidity())<2.5, fwdHiggs = !ctrlHiggs;
      int vbfTopo = vbfTopology_Stage1_1(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 ( higgs.pt()>200 ) return GG2H_PTH_GT200;
    if (Njets==0)  return higgs.pt()<10 ? GG2H_0J_PTH_0_10 : GG2H_0J_PTH_GT10;
    if (Njets==1)  return Category(GG2H_1J_PTH_0_60+getBin(higgs.pt(),{0,60,120,200}));
    if (Njets>1){
      //VBF topology
      if(vbfTopo) return Category(GG2H_MJJ_350_700_PTHJJ_0_25+vbfTopo-1);
      //Njets >= 2jets without VBF topology (mjj<350)
      return Category(GG2H_GE2J_MJJ_0_350_PTH_0_60+getBin(higgs.pt(),{0,60,120,200}));
    }
      }

      // 2. Electroweak qq->Hqq Stage 1.1 categories
      else if (prodMode==HTXS::VBF || ( isVH(prodMode) && quarkDecay(V)) ) {
    if (std::abs(higgs.rapidity())>2.5) return QQ2HQQ_FWDH;
    int Njets=jets.size();
    if (Njets==0)        return QQ2HQQ_0J;
    else if (Njets==1)   return QQ2HQQ_1J;
    else if (Njets>=2) {
      double mjj = (jets[0].mom()+jets[1].mom()).mass();
      if ( mjj < 60 )      return QQ2HQQ_MJJ_0_60;
      else if ( 60 < mjj && mjj < 120 ) return QQ2HQQ_MJJ_60_120;
      else if ( 120 < mjj && mjj < 350 ) return QQ2HQQ_MJJ_120_350;
      else if (  mjj > 350 ) {
            if (higgs.pt()>200) return QQ2HQQ_MJJ_GT350_PTH_GT200;
            if(vbfTopo) return Category(QQ2HQQ_MJJ_GT350_PTH_GT200+vbfTopo);
      }
    }
      }
      // 3. WH->Hlv categories
      else if (prodMode==HTXS::WH) {
        if (fwdHiggs) return QQ2HLNU_FWDH;
        else if (V.pt()<75) return QQ2HLNU_PTV_0_75;
        else if (V.pt()<150) return QQ2HLNU_PTV_75_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()<75) return QQ2HLL_PTV_0_75;
        else if (V.pt()<150) return QQ2HLL_PTV_75_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;
        else if (V.pt()<75) return GG2HLL_PTV_0_75;
        else if (V.pt()<150) return GG2HLL_PTV_75_150;
        else if (V.pt()>250) return GG2HLL_PTV_GT250;
        return jets.empty() ? GG2HLL_PTV_150_250_0J : GG2HLL_PTV_150_250_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;
    }

    HTXS::Stage1_1_Fine::Category getStage1_1_Fine_Category(const HTXS::HiggsProdMode prodMode,
                                const Particle &higgs,
                                const Jets &jets,
                                const Particle &V) {
      using namespace HTXS::Stage1_1_Fine;
      int Njets=jets.size(), ctrlHiggs = std::abs(higgs.rapidity())<2.5, fwdHiggs = !ctrlHiggs;
      int vbfTopo = vbfTopology_Stage1_1_Fine(jets,higgs);

      // 1. GGF Stage 1.1 categories
      //    Following YR4 write-up: XXXXX
      if (prodMode==HTXS::GGF || (prodMode==HTXS::GG2ZH && quarkDecay(V)) ) {
    if (fwdHiggs)        return GG2H_FWDH;
    if ( higgs.pt()>200 ) return GG2H_PTH_GT200;
    if (Njets==0)  return higgs.pt()<10 ? GG2H_0J_PTH_0_10 : GG2H_0J_PTH_GT10;
    if (Njets==1)  return Category(GG2H_1J_PTH_0_60+getBin(higgs.pt(),{0,60,120,200}));
    if (Njets>1){
      //double mjj = (jets[0].mom()+jets[1].mom()).mass();
      double pTHjj = (jets[0].momentum()+jets[1].momentum()+higgs.momentum()).pt();
      //VBF topology
      if(vbfTopo) return Category(GG2H_MJJ_350_700_PTHJJ_0_25+vbfTopo-1);
      //Njets >= 2jets without VBF topology (mjj<350)
      if (pTHjj<25) return Category(GG2H_GE2J_MJJ_0_350_PTH_0_60_PTHJJ_0_25+getBin(higgs.pt(),{0,60,120,200}));
      else return Category(GG2H_GE2J_MJJ_0_350_PTH_0_60_PTHJJ_GT25+getBin(higgs.pt(),{0,60,120,200}));
    }
      }

      // 2. Electroweak qq->Hqq Stage 1.1 categories
      else if (prodMode==HTXS::VBF || ( isVH(prodMode) && quarkDecay(V)) ) {
    if (std::abs(higgs.rapidity())>2.5) return QQ2HQQ_FWDH;
    int Njets=jets.size();
    if (Njets==0)        return QQ2HQQ_0J;
    else if (Njets==1)   return QQ2HQQ_1J;
    else if (Njets>=2) {
      double mjj = (jets[0].mom()+jets[1].mom()).mass();
      double pTHjj = (jets[0].momentum()+jets[1].momentum()+higgs.momentum()).pt();
      if (mjj<350){
            if (pTHjj<25) return Category(QQ2HQQ_MJJ_0_60_PTHJJ_0_25+getBin(mjj,{0,60,120,350}));
            else return Category(QQ2HQQ_MJJ_0_60_PTHJJ_GT25+getBin(mjj,{0,60,120,350}));
      } else { //mjj>350 GeV
            if (higgs.pt()<200){
          return Category(QQ2HQQ_MJJ_350_700_PTHJJ_0_25+vbfTopo-1);
            } else {
          return Category(QQ2HQQ_PTH_GT200_MJJ_350_700_PTHJJ_0_25+vbfTopo-1);
            }
      }
    }
      }
      // 3. WH->Hlv categories
      else if (prodMode==HTXS::WH) {
        if (fwdHiggs) return QQ2HLNU_FWDH;
        int Njets=jets.size();
        if (Njets==0) return Category(QQ2HLNU_PTV_0_75_0J+getBin(V.pt(),{0,75,150,250,400}));
        if (Njets==1) return Category(QQ2HLNU_PTV_0_75_1J+getBin(V.pt(),{0,75,150,250,400}));
        return Category(QQ2HLNU_PTV_0_75_GE2J+getBin(V.pt(),{0,75,150,250,400}));
      }
      // 4. qq->ZH->llH categories
      else if (prodMode==HTXS::QQ2ZH) {
        if (fwdHiggs) return QQ2HLL_FWDH;
        int Njets=jets.size();
        if (Njets==0) return Category(QQ2HLL_PTV_0_75_0J+getBin(V.pt(),{0,75,150,250,400}));
        if (Njets==1) return Category(QQ2HLL_PTV_0_75_1J+getBin(V.pt(),{0,75,150,250,400}));
        return Category(QQ2HLL_PTV_0_75_GE2J+getBin(V.pt(),{0,75,150,250,400}));
      }
      // 5. gg->ZH->llH categories
      else if (prodMode==HTXS::GG2ZH ) {
        if (fwdHiggs) return GG2HLL_FWDH;
        int Njets=jets.size();
        if (Njets==0) return Category(GG2HLL_PTV_0_75_0J+getBin(V.pt(),{0,75,150,250,400}));
        if (Njets==1) return Category(GG2HLL_PTV_0_75_1J+getBin(V.pt(),{0,75,150,250,400}));
        return Category(GG2HLL_PTV_0_75_GE2J+getBin(V.pt(),{0,75,150,250,400}));
      }
      // 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;
    }





    void setHiggsProdMode( HTXS::HiggsProdMode prodMode ){ m_HiggsProdMode = prodMode; }

    void init() override {
      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");
    }
        
    // Perform the per-event analysis
    void analyze(const Event& event) override {

      // 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);
      }
    }
    
    void printClassificationSummary(){
      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 (" ====================================================== ");
    }


    void finalize() override {
      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);
    }
    
    /*
     *  initialize histograms
     */

    void initializeHistos(){
      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);
    }

    /*
     *    initialize private members used in the classification procedure
     */
    
  private:
    double m_sumw;
    HTXS::HiggsProdMode m_HiggsProdMode;
    std::map<HTXS::ErrorCode,size_t> m_errorCount;
    Histo1DPtr hist_stage0;
    Histo1DPtr hist_stage1_pTjet25, hist_stage1_pTjet30;
    Histo1DPtr hist_pT_Higgs, hist_y_Higgs;
    Histo1DPtr hist_pT_V, hist_pT_jet1;
    Histo1DPtr hist_deltay_jj, hist_dijet_mass, hist_pT_Hjj;
    Histo1DPtr hist_Njets25, hist_Njets30;
  };

  // the PLUGIN only needs to be decleared when running standalone Rivet
  // and causes compilation / linking issues if included in Athena / RootCore
  //check for Rivet environment variable RIVET_ANALYSIS_PATH
#ifdef RIVET_ANALYSIS_PATH
  // The hook for the plugin system
  DECLARE_RIVET_PLUGIN(HiggsTemplateCrossSections);
#endif

}

#endif