diff --git a/Generators/GenInterfaces/GenInterfaces/IHiggsTruthCategoryTool.h b/Generators/GenInterfaces/GenInterfaces/IHiggsTruthCategoryTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..97285a8c8f6a5526aa8fe377ed87f1c9991ced47
--- /dev/null
+++ b/Generators/GenInterfaces/GenInterfaces/IHiggsTruthCategoryTool.h
@@ -0,0 +1,31 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+/*
+ * Dual-use tool interface for Rivet routine for classifying MC events according to the Higgs template cross section categories
+ * Authors: Jim Lacey (Carleton University)
+ * <james.lacey@cern.ch,jlacey@physics.carleton.ca>
+ */
+
+#ifndef GENINTERFACES_IHIGGSTRUTHCATEGORYTOOL_H
+#define GENINTERFACES_IHIGGSTRUTHCATEGORYTOOL_H 1
+
+#include "AsgTools/IAsgTool.h"
+#include "HepMC/GenEvent.h"
+
+namespace HTXS {
+ class HiggsClassification;
+}
+
+class IHiggsTruthCategoryTool : public virtual asg::IAsgTool {
+ public:
+ ASG_TOOL_INTERFACE( IHiggsTruthCategoryTool ) //declares the interface to athena
+ virtual ~IHiggsTruthCategoryTool() {};
+ public:
+ virtual StatusCode initialize() = 0;
+ virtual StatusCode finalize () = 0;
+ virtual HTXS::HiggsClassification* getHiggsTruthCategoryObject(const HepMC::GenEvent& HepMCEvent, const HTXS::HiggsProdMode prodMode)=0;
+};
+
+#endif //> !GENINTERFACES_IHIGGSTRUTHCATEGORYTOOL_H
diff --git a/Generators/TruthRivetTools/CMakeLists.txt b/Generators/TruthRivetTools/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..5ffd61b530ec3f72fbf013ec2f6f3d74d80c85a1
--- /dev/null
+++ b/Generators/TruthRivetTools/CMakeLists.txt
@@ -0,0 +1,44 @@
+################################################################################
+# Package: TruthRivetTools
+################################################################################
+
+# Declare the name of the package:
+atlas_subdir( TruthRivetTools )
+
+# Declare the dependencies of the package:
+atlas_depends_on_subdirs(
+ PUBLIC
+ Control/AthToolSupport/AsgTools
+ Generators/AtlasHepMC
+ PRIVATE
+ GaudiKernel
+ GenInterfaces )
+
+# External dependencies:
+find_package( HepMC )
+find_package( Rivet )
+find_package( YODA )
+find_package( FastJet )
+find_package( ROOT COMPONENTS Core Physics )
+find_package( GSL )
+
+# Remove the --as-needed linker flags:
+atlas_disable_as_needed()
+
+# Component(s) in the package:
+atlas_add_library( TruthRivetToolsLib
+ Root/*.cxx TruthRivetTools/*.h
+ SHARED
+ PUBLIC_HEADERS TruthRivetTools
+ INCLUDE_DIRS ${HEPMC_INCLUDE_DIRS} ${RIVET_INCLUDE_DIRS} ${YODA_INCLUDE_DIRS}
+ ${FASTJET_INCLUDE_DIRS} ${ROOT_INCLUDE_DIRS} ${GSL_INCLUDE_DIRS}
+ LINK_LIBRARIES ${HEPMC_LIBRARIES} ${RIVET_LIBRARIES} ${YODA_LIBRARIES}
+ ${FASTJET_LIBRARIES} ${ROOT_LIBRARIES} ${GSL_LIBRARIES} AsgTools AtlasHepMCLib )
+
+atlas_add_component( TruthRivetTools
+ src/components/*.cxx
+ LINK_LIBRARIES TruthRivetToolsLib GaudiKernel )
+
+atlas_add_dictionary( TruthRivetToolsDict
+ TruthRivetTools/TruthRivetToolsDict.h TruthRivetTools/selection.xml
+ LINK_LIBRARIES TruthRivetToolsLib )
diff --git a/Generators/TruthRivetTools/Root/HiggsTruthCategoryTool.cxx b/Generators/TruthRivetTools/Root/HiggsTruthCategoryTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..e7c854ff32e480328be3ca7d57142b074d799b35
--- /dev/null
+++ b/Generators/TruthRivetTools/Root/HiggsTruthCategoryTool.cxx
@@ -0,0 +1,59 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// TruthRivetTools includes
+#include "TruthRivetTools/HiggsTruthCategoryTool.h"
+
+HiggsTruthCategoryTool::HiggsTruthCategoryTool( const std::string& name)
+: asg::AsgTool( name ),
+ rivetAnaHandler(nullptr),
+ higgsTemplateCrossSections(nullptr)
+{
+ // cannot be set to true until the issue with the beam protons in the truth event record is fixed..
+ // see JIRA ticket: https://its.cern.ch/jira/browse/ATLASRECTS-3072?filter=-2
+ m_outHistos = false;
+ m_isInitialized = false;
+}
+
+
+StatusCode HiggsTruthCategoryTool::initialize() {
+ ATH_MSG_INFO ("Initializing " << name() << "...");
+ // Rivet analysis :: Higgs truth event classifier class
+ higgsTemplateCrossSections = new Rivet::HiggsTemplateCrossSections();
+ // crreate an instance of the Rivet analysis handler
+ rivetAnaHandler = new Rivet::AnalysisHandler();
+ // Add the Higgs truth classifier class to the handler
+ rivetAnaHandler->addAnalysis(&(*higgsTemplateCrossSections));
+ return StatusCode::SUCCESS;
+}
+
+StatusCode HiggsTruthCategoryTool :: finalize () {
+ ATH_MSG_INFO (" ====================================================== ");
+ ATH_MSG_INFO (" ---- Finalizing" << name() << "...");
+ ATH_MSG_INFO (" ====================================================== ");
+ if ( !m_outHistos )
+ higgsTemplateCrossSections->printClassificationSummary( );
+ else{
+ // TODO:: update the tool properly deal with output files/paths
+ rivetAnaHandler->finalize();
+ rivetAnaHandler->writeData("HiggsTruthCategoryTool.yoda");
+ }
+ ATH_MSG_INFO (" ====================================================== ");
+ return StatusCode::SUCCESS;
+}
+
+HTXS::HiggsClassification* HiggsTruthCategoryTool :: getHiggsTruthCategoryObject (const HepMC::GenEvent& HepMCEvent, const HTXS::HiggsProdMode prodMode){
+ if ( !m_isInitialized ) {
+ higgsTemplateCrossSections->setHiggsProdMode(prodMode);
+ rivetAnaHandler->init(HepMCEvent);
+ m_isInitialized=true;
+ }
+ // fill histos if flag is specified
+ if ( m_outHistos ) rivetAnaHandler->analyze(HepMCEvent);
+ // get the category output object containing the template cross section category,
+ // and Higgs, V-boson, jets 4-vectors
+ const Rivet::HiggsClassification htxs_cat_rivet = higgsTemplateCrossSections->classifyEvent(HepMCEvent,prodMode);
+ return HTXS::Rivet2Root(htxs_cat_rivet);
+}
+
diff --git a/Generators/TruthRivetTools/Root/LinkDef.h b/Generators/TruthRivetTools/Root/LinkDef.h
new file mode 100644
index 0000000000000000000000000000000000000000..e5ed76f140e9866849dbf1b6c460ca5a852be5cd
--- /dev/null
+++ b/Generators/TruthRivetTools/Root/LinkDef.h
@@ -0,0 +1,16 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include <TruthRivetTools/HiggsTruthCategoryTool.h>
+
+#ifdef __CINT__
+#pragma link off all globals;
+#pragma link off all classes;
+#pragma link off all functions;
+#pragma link C++ nestedclass;
+#endif
+
+#ifdef __CINT__
+#pragma link C++ class HiggsTruthCategoryTool+;
+#endif
diff --git a/Generators/TruthRivetTools/TruthRivetTools/HiggsTemplateCrossSections.h b/Generators/TruthRivetTools/TruthRivetTools/HiggsTemplateCrossSections.h
new file mode 100644
index 0000000000000000000000000000000000000000..b875bcdf982fefb9adcd358d649faec2fe7db676
--- /dev/null
+++ b/Generators/TruthRivetTools/TruthRivetTools/HiggsTemplateCrossSections.h
@@ -0,0 +1,791 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef TRUTHRIVETTOOLS_HIGGSTEMPLATECROSSSECTIONS_H
+#define TRUTHRIVETTOOLS_HIGGSTEMPLATECROSSSECTIONS_H
+
+// -*- C++ -*-
+#include "Rivet/Analysis.hh"
+#include "Rivet/Particle.hh"
+#include "Rivet/Projections/FastJets.hh"
+
+// Definition of the StatusCode and Category enums
+// Note: the Template XSec Defs *depends* on having included
+// the TLorentzVector header *before* it is included -- it
+// uses the include guard from TLorentzVector to decide
+// what is available
+#include "TLorentzVector.h"
+#include "TruthRivetTools/HiggsTemplateCrossSectionsDefs.h"
+#include "AtlasHepMC/Relatives.h"
+#include "AtlasHepMC/GenVertex.h"
+#include "AtlasHepMC/GenParticle.h"
+#include "AtlasHepMC/GenRanges.h"
+
+namespace Rivet {
+
+ /// @class HiggsTemplateCrossSections
+ /// @brief Rivet routine for classifying MC events according to the Higgs template cross section categories
+ /// @author Jim Lacey (DESY) <james.lacey@cern.ch,jlacey@desy.de>
+ /// @author Dag Gillberg (Carleton University) <dag.gillberg@cern.ch>
+ class HiggsTemplateCrossSections : public Analysis {
+ public:
+ // Constructor
+ HiggsTemplateCrossSections()
+ : Analysis("HiggsTemplateCrossSections"),
+ m_HiggsProdMode(HTXS::UNKNOWN) {}
+
+ public:
+
+ /// @name Utility methods
+ /// Methods to identify the Higgs boson and
+ /// associated vector boson and to build jets
+ /// @{
+
+ /// follow a "propagating" particle and return its last instance
+ Particle getLastInstance(Particle ptcl) {
+ if ( ptcl.genParticle()->end_vertex() ) {
+ if ( !hasChild(ptcl.genParticle(),ptcl.pid()) ) return ptcl;
+ else return getLastInstance(ptcl.children()[0]);
+ }
+ return ptcl;
+ }
+
+ /// @brief Whether particle p originate from any of the ptcls
+ bool originateFrom(const Particle& p, const Particles& ptcls ) {
+ const HepMC::GenVertex* prodVtx = p.genParticle()->production_vertex();
+ if (prodVtx == nullptr) return false;
+ // for each ancestor, check if it matches any of the input particles
+ for (auto ancestor:Rivet::HepMCUtils::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;
+ }
+
+ /// @brief Whether particle p originates from p2
+ bool originateFrom(const Particle& p, const Particle& p2 ) {
+ Particles ptcls = {p2}; return originateFrom(p,ptcls);
+ }
+
+ /// @brief Checks whether the input particle has a child with a given PDGID
+ bool hasChild(const HepMC::GenParticle *ptcl, int pdgID) {
+ for (const Particle& child:Particle(*ptcl).children())
+ if (child.pid()==pdgID) return true;
+ return false;
+ }
+
+ /// @brief Checks whether the input particle has a parent with a given PDGID
+ bool hasParent(const HepMC::GenParticle *ptcl, int pdgID) {
+ for (auto parent:Rivet::HepMCUtils::particles(ptcl->production_vertex(),HepMC::parents))
+ if (parent->pdg_id()==pdgID) return true;
+ return false;
+ }
+
+ /// @brief Return true is particle decays to quarks
+ bool quarkDecay(const Particle &p) {
+ for (const Particle& child:p.children())
+ if (PID::isQuark(child.pid())) return true;
+ return false;
+ }
+
+ /// @brief Return true if particle decays to charged leptons.
+ bool ChLeptonDecay(const Particle &p) {
+ for (const Particle& child:p.children())
+ if (PID::isChLepton(child.pid())) return true;
+ return false;
+ }
+
+ /// @brief Returns the classification object with the error code set.
+ /// Prints an warning message, and keeps track of number of errors
+ 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 int Nwarnings = 0;
+ if ( msg!="" && ++Nwarnings < NmaxWarnings )
+ MSG_WARNING(msg);
+
+ return cat;
+ }
+ /// @}
+
+ /// @brief Main classificaion method.
+ 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_2_cat_pTjet25GeV = HTXS::Stage1_2::UNKNOWN;
+ cat.stage1_2_cat_pTjet30GeV = HTXS::Stage1_2::UNKNOWN;
+ cat.stage1_2_fine_cat_pTjet25GeV = HTXS::Stage1_2_Fine::UNKNOWN;
+ cat.stage1_2_fine_cat_pTjet30GeV = HTXS::Stage1_2_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.
+ */
+
+ HepMC::GenVertex *HSvtx = event.genEvent()->signal_process_vertex();
+ int Nhiggs=0;
+ for ( const HepMC::GenParticle *ptcl : Rivet::HepMCUtils::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:Rivet::HepMCUtils::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:Rivet::HepMCUtils::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);
+ }
+ }
+
+ 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 : Rivet::HepMCUtils::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.pid()) ) Ws += getLastInstance(child);
+ }
+ }
+
+ // Make sure result make sense
+ if ( (prodMode==HTXS::TTH && Ws.size()<2) || (prodMode==HTXS::TH && Ws.size()<1 ) )
+ 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 Particles FS = apply<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.size() && 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.isZ2vvDecay = false;
+ if( (prodMode==HTXS::GG2ZH || prodMode==HTXS::QQ2ZH) && !quarkDecay(cat.V) && !ChLeptonDecay(cat.V) ) cat.isZ2vvDecay = true;
+ 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_2_cat_pTjet25GeV = getStage1_2_Category(prodMode,cat.higgs,cat.jets25,cat.V);
+ cat.stage1_2_cat_pTjet30GeV = getStage1_2_Category(prodMode,cat.higgs,cat.jets30,cat.V);
+ cat.stage1_2_fine_cat_pTjet25GeV = getStage1_2_Fine_Category(prodMode,cat.higgs,cat.jets25,cat.V);
+ cat.stage1_2_fine_cat_pTjet30GeV = getStage1_2_Fine_Category(prodMode,cat.higgs,cat.jets30,cat.V);
+ cat.errorCode = HTXS::SUCCESS; ++m_errorCount[HTXS::SUCCESS];
+
+ return cat;
+ }
+
+ /// @name Categorization methods
+ /// Methods to assign the truth category based
+ /// on the identified Higgs boson and associated
+ /// vector bosons and/or reconstructed jets
+ /// @{
+
+ /// @brief Return bin index of x given the provided bin edges. 0=first bin, -1=underflow bin.
+ int getBin(double x, std::vector<double> bins) {
+ if (bins.size()==0||x<bins[0]) return -1; // should not happen!
+ for (size_t i=1;i<bins.size();++i)
+ if (x<bins[i]) return i-1;
+ return bins.size()-1;
+ }
+
+ /// @brief 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
+ 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;
+ }
+ /// @brief VBF topology selection
+ /// 0 = fail loose selection: m_jj > 350 GeV
+ /// 1 pass loose, but fail additional cut pT(Hjj)<25. 2 pass pT(Hjj)>25 selection
+ /// 3 pass tight (m_jj>700 GeV), but fail additional cut pT(Hjj)<25. 4 pass pT(Hjj)>25 selection
+ int vbfTopology_Stage1_2(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;
+ }
+ /// @brief VBF topology selection for Stage1_2
+ /// 0 = fail loose selection: m_jj > 350 GeV
+ /// 1 pass loose, but fail additional cut pT(Hjj)<25. 2 pass pT(Hjj)>25 selection
+ /// 3 pass 700<m_jj<1000 GeV, but fail additional cut pT(Hjj)<25. 4 pass pT(Hjj)>25 selection
+ /// 5 pass 1000<m_jj<1500 GeV, but fail additional cut pT(Hjj)<25. 6 pass pT(Hjj)>25 selection
+ /// 7 pass m_jj>1500 GeV, but fail additional cut pT(Hjj)<25. 8 pass pT(Hjj)>25 selection
+ int vbfTopology_Stage1_2_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;
+ }
+
+ /// @brief Whether the Higgs is produced in association with a vector boson (VH)
+ bool isVH(HTXS::HiggsProdMode p) { return p==HTXS::WH || p==HTXS::QQ2ZH || p==HTXS::GG2ZH; }
+
+ /// @brief Stage-0 HTXS categorization
+ 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);
+ }
+
+ /// @brief Stage-1 categorization
+ 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.size() ? 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.size()==0 ? 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.size()==0 ? 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.size()==0) 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;
+ }
+
+ /// @brief Stage-1.2 categorization
+ HTXS::Stage1_2::Category getStage1_2_Category(const HTXS::HiggsProdMode prodMode,
+ const Particle &higgs,
+ const Jets &jets,
+ const Particle &V) {
+ using namespace HTXS::Stage1_2;
+ int Njets=jets.size(), ctrlHiggs = std::abs(higgs.rapidity())<2.5, fwdHiggs = !ctrlHiggs;
+ int vbfTopo = vbfTopology_Stage1_2(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 Category(GG2H_PTH_200_300+getBin(higgs.pt(),{200,300,450,650}));
+ 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_GE2J_MJJ_350_700_PTH_0_200_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.2 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_GE2J_MJJ_0_60;
+ else if ( 60 < mjj && mjj < 120 ) return QQ2HQQ_GE2J_MJJ_60_120;
+ else if ( 120 < mjj && mjj < 350 ) return QQ2HQQ_GE2J_MJJ_120_350;
+ else if ( mjj > 350 ) {
+ if (higgs.pt()>200) return QQ2HQQ_GE2J_MJJ_GT350_PTH_GT200;
+ if(vbfTopo) return Category(QQ2HQQ_GE2J_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.size()==0 ? 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.size()==0 ? 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.size()==0 ? GG2HLL_PTV_150_250_0J : GG2HLL_PTV_150_250_GE1J;
+ }
+ // 6.ttH,bbH,tH categories
+ else if (prodMode==HTXS::TTH) {
+ if (fwdHiggs) return TTH_FWDH;
+ else return Category(TTH_PTH_0_60+getBin(higgs.pt(),{0,60,120,200,300}));
+ }
+ else if (prodMode==HTXS::BBH) return Category(BBH_FWDH+ctrlHiggs);
+ else if (prodMode==HTXS::TH ) return Category(TH_FWDH+ctrlHiggs);
+ return UNKNOWN;
+ }
+
+ /// @brief Stage-1.2_Fine categorization
+ HTXS::Stage1_2_Fine::Category getStage1_2_Fine_Category(const HTXS::HiggsProdMode prodMode,
+ const Particle &higgs,
+ const Jets &jets,
+ const Particle &V) {
+ using namespace HTXS::Stage1_2_Fine;
+ int Njets=jets.size(), ctrlHiggs = std::abs(higgs.rapidity())<2.5, fwdHiggs = !ctrlHiggs;
+ int vbfTopo = vbfTopology_Stage1_2_Fine(jets,higgs);
+
+ // 1. GGF Stage 1.2 categories
+ // Following YR4 write-up: XXXXX
+ if (prodMode==HTXS::GGF || (prodMode==HTXS::GG2ZH && quarkDecay(V)) ) {
+ if (fwdHiggs) return GG2H_FWDH;
+ if ( higgs.pt()>200 ){
+ if (Njets>0){
+ double pTHj = (jets[0].momentum()+higgs.momentum()).pt();
+ if( pTHj/higgs.pt()>0.15 ) return Category(GG2H_PTH_200_300_PTHJoverPTH_GT15+getBin(higgs.pt(),{200,300,450,650}));
+ else return Category(GG2H_PTH_200_300_PTHJoverPTH_0_15+getBin(higgs.pt(),{200,300,450,650}));
+ }
+ else return Category(GG2H_PTH_200_300_PTHJoverPTH_0_15+getBin(higgs.pt(),{200,300,450,650}));
+ }
+ 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_GE2J_MJJ_350_700_PTH_0_200_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.2 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_GE2J_MJJ_0_60_PTHJJ_0_25+getBin(mjj,{0,60,120,350}));
+ else return Category(QQ2HQQ_GE2J_MJJ_0_60_PTHJJ_GT25+getBin(mjj,{0,60,120,350}));
+ } else { //mjj>350 GeV
+ if (higgs.pt()<200){
+ return Category(QQ2HQQ_GE2J_MJJ_350_700_PTH_0_200_PTHJJ_0_25+vbfTopo-1);
+ } else {
+ return Category(QQ2HQQ_GE2J_MJJ_350_700_PTH_GT200_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) {
+ if (fwdHiggs) return TTH_FWDH;
+ else return Category(TTH_PTH_0_60+getBin(higgs.pt(),{0,60,120,200,300,450}));
+ }
+ else if (prodMode==HTXS::BBH) return Category(BBH_FWDH+ctrlHiggs);
+ else if (prodMode==HTXS::TH ) return Category(TH_FWDH+ctrlHiggs);
+ return UNKNOWN;
+ }
+
+ /// @}
+
+
+ /// @name Default Rivet analysis methods and steering methods
+ /// @{
+
+ /// @brief Sets the Higgs production mode
+ void setHiggsProdMode( HTXS::HiggsProdMode prodMode ){ m_HiggsProdMode = prodMode; }
+
+ /// @brief default Rivet Analysis::init method
+ /// Booking of histograms, initializing Rivet projection
+ /// Extracts Higgs production mode from shell variable if not set manually using setHiggsProdMode
+ void init() {
+ 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;
+ declare(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) {
+
+ // get the classification
+ HiggsClassification cat = classifyEvent(event,m_HiggsProdMode);
+
+ // Fill histograms: categorization --> linerize the categories
+ const double weight = 1.; // Event weights are now all 1 in Rivet
+ m_sumw += weight;
+
+ int F=cat.stage0_cat%10, P=cat.stage1_cat_pTjet30GeV/100;
+ m_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
+ static vector<int> offset({0,1,13,19,24,29,33,35,37,39});
+ int off = offset[P];
+ // 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
+ static vector<int> offset1_2({0,1,18,29,35,41,47,53,55,57});
+ int off1_2 = offset1_2[P];
+ // 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
+ static vector<int> offset1_2_Fine({0,1,29,54,70,86,102,109,111,113});
+ int off1_2_Fine = offset1_2_Fine[P];
+
+ m_hist_stage1_pTjet25->fill(cat.stage1_cat_pTjet25GeV%100 + off, weight);
+ m_hist_stage1_pTjet30->fill(cat.stage1_cat_pTjet30GeV%100 + off, weight);
+ m_hist_stage1_2_pTjet25->fill(cat.stage1_2_cat_pTjet25GeV%100 + off1_2, weight);
+ m_hist_stage1_2_pTjet30->fill(cat.stage1_2_cat_pTjet30GeV%100 + off1_2, weight);
+ m_hist_stage1_2_fine_pTjet25->fill(cat.stage1_2_fine_cat_pTjet25GeV%100 + off1_2_Fine, weight);
+ m_hist_stage1_2_fine_pTjet30->fill(cat.stage1_2_fine_cat_pTjet30GeV%100 + off1_2_Fine, weight);
+
+ // Fill histograms: variables used in the categorization
+ m_hist_pT_Higgs->fill(cat.higgs.pT(),weight);
+ m_hist_y_Higgs->fill(cat.higgs.rapidity(),weight);
+ m_hist_pT_V->fill(cat.V.pT(),weight);
+
+ m_hist_Njets25->fill(cat.jets25.size(),weight);
+ m_hist_Njets30->fill(cat.jets30.size(),weight);
+
+ m_hist_isZ2vv->fill(cat.isZ2vvDecay, weight);
+
+ // Jet variables. Use jet collection with pT threshold at 30 GeV
+ if (cat.jets30.size()) m_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();
+ m_hist_deltay_jj->fill(std::abs(j1.rapidity()-j2.rapidity()),weight);
+ m_hist_dijet_mass->fill((j1+j2).mass(),weight);
+ m_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() {
+ printClassificationSummary();
+ double sf = m_sumw>0?1.0/m_sumw:1.0;
+ for (auto hist:{m_hist_stage0,m_hist_stage1_pTjet25,m_hist_stage1_pTjet30,m_hist_stage1_2_pTjet25,m_hist_stage1_2_pTjet30,m_hist_stage1_2_fine_pTjet25,m_hist_stage1_2_fine_pTjet30,
+ m_hist_Njets25,m_hist_Njets30,m_hist_pT_Higgs,m_hist_y_Higgs,m_hist_pT_V,m_hist_pT_jet1,m_hist_deltay_jj,m_hist_dijet_mass,m_hist_pT_Hjj,m_hist_isZ2vv})
+ scale(hist, sf);
+ }
+
+ /*
+ * initialize histograms
+ */
+
+ void initializeHistos(){
+ book(m_hist_stage0,"HTXS_stage0",20,0,20);
+ book(m_hist_stage1_pTjet25,"HTXS_stage1_pTjet25",40,0,40);
+ book(m_hist_stage1_pTjet30,"HTXS_stage1_pTjet30",40,0,40);
+ book(m_hist_stage1_2_pTjet25,"HTXS_stage1_2_pTjet25",57,0,57);
+ book(m_hist_stage1_2_pTjet30,"HTXS_stage1_2_pTjet30",57,0,57);
+ book(m_hist_stage1_2_fine_pTjet25,"HTXS_stage1_2_fine_pTjet25",113,0,113);
+ book(m_hist_stage1_2_fine_pTjet30,"HTXS_stage1_2_fine_pTjet30",113,0,113);
+ book(m_hist_pT_Higgs,"pT_Higgs",80,0,400);
+ book(m_hist_y_Higgs,"y_Higgs",80,-4,4);
+ book(m_hist_pT_V,"pT_V",80,0,400);
+ book(m_hist_pT_jet1,"pT_jet1",80,0,400);
+ book(m_hist_deltay_jj ,"deltay_jj",50,0,10);
+ book(m_hist_dijet_mass,"m_jj",50,0,2000);
+ book(m_hist_pT_Hjj,"pT_Hjj",50,0,250);
+ book(m_hist_Njets25,"Njets25",10,0,10);
+ book(m_hist_Njets30,"Njets30",10,0,10);
+ book(m_hist_isZ2vv,"isZ2vv",2,0,2);
+ }
+ /// @}
+
+ /*
+ * 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 m_hist_stage0;
+ Histo1DPtr m_hist_stage1_pTjet25, m_hist_stage1_pTjet30;
+ Histo1DPtr m_hist_stage1_2_pTjet25, m_hist_stage1_2_pTjet30;
+ Histo1DPtr m_hist_stage1_2_fine_pTjet25, m_hist_stage1_2_fine_pTjet30;
+ Histo1DPtr m_hist_pT_Higgs, m_hist_y_Higgs;
+ Histo1DPtr m_hist_pT_V, m_hist_pT_jet1;
+ Histo1DPtr m_hist_deltay_jj, m_hist_dijet_mass, m_hist_pT_Hjj;
+ Histo1DPtr m_hist_Njets25, m_hist_Njets30;
+ Histo1DPtr m_hist_isZ2vv;
+ };
+
+ // 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
+
diff --git a/Generators/TruthRivetTools/TruthRivetTools/HiggsTemplateCrossSectionsDefs.h b/Generators/TruthRivetTools/TruthRivetTools/HiggsTemplateCrossSectionsDefs.h
new file mode 100644
index 0000000000000000000000000000000000000000..10492c7aab971b9887d91e8a2f74a4982efec655
--- /dev/null
+++ b/Generators/TruthRivetTools/TruthRivetTools/HiggsTemplateCrossSectionsDefs.h
@@ -0,0 +1,549 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef TRUTHRIVETTOOLS_HIGGSTEMPLATECROSSSECTIONSDEFS_H
+#define TRUTHRIVETTOOLS_HIGGSTEMPLATECROSSSECTIONSDEFS_H 1
+
+/// Higgs Template Cross Section namespace
+namespace HTXS {
+
+ /// Error code: whether the classification was successful or failed
+ enum ErrorCode {
+ UNDEFINED=-99,
+ SUCCESS = 0, ///< successful classification
+ PRODMODE_DEFINED = 1, ///< production mode not defined
+ MOMENTUM_CONSERVATION = 2, ///< failed momentum conservation
+ HIGGS_IDENTIFICATION = 3, ///< failed to identify Higgs boson
+ HIGGS_DECAY_IDENTIFICATION = 4, ///< failed to identify Higgs boson decay products
+ HS_VTX_IDENTIFICATION = 5, ///< failed to identify hard scatter vertex
+ VH_IDENTIFICATION = 6, ///< failed to identify associated vector boson
+ VH_DECAY_IDENTIFICATION = 7, ///< failed to identify associated vector boson decay products
+ TOP_W_IDENTIFICATION = 8 ///< failed to identify top decay
+ };
+
+ /// Higgs production modes, corresponding to input sample
+ enum HiggsProdMode {
+ UNKNOWN = 0,
+ GGF = 1, VBF = 2, WH = 3, QQ2ZH = 4, GG2ZH = 5,
+ TTH = 6, BBH = 7, TH = 8
+ };
+
+ /// Additional identifier flag for TH production modes
+ enum tH_type { noTH=0, THQB=1, TWH=2 };
+
+ /// Two digit number of format PF
+ /// P is digit for the physics process
+ /// and F is 0 for |yH|>2.5 and 11 for |yH|<2.5 ("in fiducial")
+
+ /// Namespace for Stage0 categorization
+ namespace Stage0 {
+ /// @enum Stage-0 ategorization: Two-digit number of format PF, with P for process and F being 0 for |yH|>2.5 and 1 for |yH|<2.5
+ enum Category {
+ UNKNOWN = 0, GG2H_FWDH = 10, GG2H = 11, VBF_FWDH = 20, VBF = 21, VH2HQQ_FWDH = 22, VH2HQQ = 23,
+ QQ2HLNU_FWDH = 30, QQ2HLNU = 31, QQ2HLL_FWDH = 40, QQ2HLL = 41, GG2HLL_FWDH = 50, GG2HLL = 51,
+ TTH_FWDH = 60, TTH = 61, BBH_FWDH = 70, BBH = 71, TH_FWDH = 80, TH = 81 };
+ }
+
+ /// Categorization Stage 1:
+ /// Three digit integer of format PF
+ /// Where P is a digit representing the process
+ /// F is a unique integer ( F < 99 ) corresponding to each Stage1 phase-space region (bin)
+ namespace Stage1 {
+ enum Category {
+ UNKNOWN = 0,
+ // Gluon fusion
+ GG2H_FWDH = 100,
+ GG2H_VBFTOPO_JET3VETO = 101, GG2H_VBFTOPO_JET3 = 102,
+ GG2H_0J = 103,
+ GG2H_1J_PTH_0_60 = 104, GG2H_1J_PTH_60_120 = 105,
+ GG2H_1J_PTH_120_200 = 106, GG2H_1J_PTH_GT200 = 107,
+ GG2H_GE2J_PTH_0_60 = 108, GG2H_GE2J_PTH_60_120 = 109,
+ GG2H_GE2J_PTH_120_200 = 110, GG2H_GE2J_PTH_GT200 = 111,
+ // "VBF"
+ QQ2HQQ_FWDH = 200,
+ QQ2HQQ_VBFTOPO_JET3VETO = 201, QQ2HQQ_VBFTOPO_JET3 = 202,
+ QQ2HQQ_VH2JET = 203, QQ2HQQ_REST = 204, QQ2HQQ_PTJET1_GT200 = 205,
+ // qq -> WH
+ QQ2HLNU_FWDH = 300,
+ QQ2HLNU_PTV_0_150 = 301,
+ QQ2HLNU_PTV_150_250_0J = 302,
+ QQ2HLNU_PTV_150_250_GE1J = 303,
+ QQ2HLNU_PTV_GT250 = 304,
+ // qq -> ZH
+ QQ2HLL_FWDH = 400,
+ QQ2HLL_PTV_0_150 = 401,
+ QQ2HLL_PTV_150_250_0J = 402,
+ QQ2HLL_PTV_150_250_GE1J = 403,
+ QQ2HLL_PTV_GT250 = 404,
+ // gg -> ZH
+ GG2HLL_FWDH = 500,
+ GG2HLL_PTV_0_150 = 501,
+ GG2HLL_PTV_GT150_0J = 502,
+ GG2HLL_PTV_GT150_GE1J = 503,
+ // ttH
+ TTH_FWDH = 600, TTH = 601,
+ // bbH
+ BBH_FWDH = 700, BBH = 701,
+ // tH
+ TH_FWDH = 800, TH = 801
+ };
+ } // namespace Stage1
+
+ /// Categorization Stage 1.2:
+ /// Three digit integer of format PF
+ /// Where P is a digit representing the process
+ /// F is a unique integer ( F < 99 ) corresponding to each Stage1_2 phase-space region (bin)
+ namespace Stage1_2 {
+ enum Category {
+ UNKNOWN = 0,
+ // Gluon fusion
+ GG2H_FWDH = 100,
+ GG2H_PTH_200_300 = 101,
+ GG2H_PTH_300_450 = 102,
+ GG2H_PTH_450_650 = 103,
+ GG2H_PTH_GT650 = 104,
+ GG2H_0J_PTH_0_10 = 105,
+ GG2H_0J_PTH_GT10 = 106,
+ GG2H_1J_PTH_0_60 = 107,
+ GG2H_1J_PTH_60_120 = 108,
+ GG2H_1J_PTH_120_200 = 109,
+ GG2H_GE2J_MJJ_0_350_PTH_0_60 = 110,
+ GG2H_GE2J_MJJ_0_350_PTH_60_120 = 111,
+ GG2H_GE2J_MJJ_0_350_PTH_120_200 = 112,
+ GG2H_GE2J_MJJ_350_700_PTH_0_200_PTHJJ_0_25 = 113,
+ GG2H_GE2J_MJJ_350_700_PTH_0_200_PTHJJ_GT25 = 114,
+ GG2H_GE2J_MJJ_GT700_PTH_0_200_PTHJJ_0_25 = 115,
+ GG2H_GE2J_MJJ_GT700_PTH_0_200_PTHJJ_GT25 = 116,
+ // "VBF"
+ QQ2HQQ_FWDH = 200,
+ QQ2HQQ_0J = 201,
+ QQ2HQQ_1J = 202,
+ QQ2HQQ_GE2J_MJJ_0_60 = 203,
+ QQ2HQQ_GE2J_MJJ_60_120 = 204,
+ QQ2HQQ_GE2J_MJJ_120_350 = 205,
+ QQ2HQQ_GE2J_MJJ_GT350_PTH_GT200 = 206,
+ QQ2HQQ_GE2J_MJJ_350_700_PTH_0_200_PTHJJ_0_25 = 207,
+ QQ2HQQ_GE2J_MJJ_350_700_PTH_0_200_PTHJJ_GT25 = 208,
+ QQ2HQQ_GE2J_MJJ_GT700_PTH_0_200_PTHJJ_0_25 = 209,
+ QQ2HQQ_GE2J_MJJ_GT700_PTH_0_200_PTHJJ_GT25 = 210,
+ // qq -> WH
+ QQ2HLNU_FWDH = 300,
+ QQ2HLNU_PTV_0_75 = 301,
+ QQ2HLNU_PTV_75_150 = 302,
+ QQ2HLNU_PTV_150_250_0J = 303,
+ QQ2HLNU_PTV_150_250_GE1J = 304,
+ QQ2HLNU_PTV_GT250 = 305,
+ // qq -> ZH
+ QQ2HLL_FWDH = 400,
+ QQ2HLL_PTV_0_75 = 401,
+ QQ2HLL_PTV_75_150 = 402,
+ QQ2HLL_PTV_150_250_0J = 403,
+ QQ2HLL_PTV_150_250_GE1J = 404,
+ QQ2HLL_PTV_GT250 = 405,
+ // gg -> ZH
+ GG2HLL_FWDH = 500,
+ GG2HLL_PTV_0_75 = 501,
+ GG2HLL_PTV_75_150 = 502,
+ GG2HLL_PTV_150_250_0J = 503,
+ GG2HLL_PTV_150_250_GE1J = 504,
+ GG2HLL_PTV_GT250 = 505,
+ // ttH
+ TTH_FWDH = 600,
+ TTH_PTH_0_60 = 601,
+ TTH_PTH_60_120 = 602,
+ TTH_PTH_120_200 = 603,
+ TTH_PTH_200_300 = 604,
+ TTH_PTH_GT300 = 605,
+ // bbH
+ BBH_FWDH = 700, BBH = 701,
+ // tH
+ TH_FWDH = 800, TH = 801
+ };
+ } // namespace Stage1_2
+
+ namespace Stage1_2_Fine {
+ enum Category {
+ UNKNOWN = 0,
+ // Gluon fusion
+ GG2H_FWDH = 100,
+ GG2H_PTH_200_300_PTHJoverPTH_0_15 = 101,
+ GG2H_PTH_300_450_PTHJoverPTH_0_15 = 102,
+ GG2H_PTH_450_650_PTHJoverPTH_0_15 = 103,
+ GG2H_PTH_GT650_PTHJoverPTH_0_15 = 104,
+ GG2H_PTH_200_300_PTHJoverPTH_GT15 = 105,
+ GG2H_PTH_300_450_PTHJoverPTH_GT15 = 106,
+ GG2H_PTH_450_650_PTHJoverPTH_GT15 = 107,
+ GG2H_PTH_GT650_PTHJoverPTH_GT15 = 108,
+ GG2H_0J_PTH_0_10 = 109,
+ GG2H_0J_PTH_GT10 = 110,
+ GG2H_1J_PTH_0_60 = 111,
+ GG2H_1J_PTH_60_120 = 112,
+ GG2H_1J_PTH_120_200 = 113,
+ GG2H_GE2J_MJJ_0_350_PTH_0_60_PTHJJ_0_25 = 114,
+ GG2H_GE2J_MJJ_0_350_PTH_60_120_PTHJJ_0_25 = 115,
+ GG2H_GE2J_MJJ_0_350_PTH_120_200_PTHJJ_0_25 = 116,
+ GG2H_GE2J_MJJ_0_350_PTH_0_60_PTHJJ_GT25 = 117,
+ GG2H_GE2J_MJJ_0_350_PTH_60_120_PTHJJ_GT25 = 118,
+ GG2H_GE2J_MJJ_0_350_PTH_120_200_PTHJJ_GT25 = 119,
+ GG2H_GE2J_MJJ_350_700_PTH_0_200_PTHJJ_0_25 = 120,
+ GG2H_GE2J_MJJ_350_700_PTH_0_200_PTHJJ_GT25 = 121,
+ GG2H_GE2J_MJJ_700_1000_PTH_0_200_PTHJJ_0_25 = 122,
+ GG2H_GE2J_MJJ_700_1000_PTH_0_200_PTHJJ_GT25 = 123,
+ GG2H_GE2J_MJJ_1000_1500_PTH_0_200_PTHJJ_0_25 = 124,
+ GG2H_GE2J_MJJ_1000_1500_PTH_0_200_PTHJJ_GT25 = 125,
+ GG2H_GE2J_MJJ_GT1500_PTH_0_200_PTHJJ_0_25 = 126,
+ GG2H_GE2J_MJJ_GT1500_PTH_0_200_PTHJJ_GT25 = 127,
+ // "VBF"
+ QQ2HQQ_FWDH = 200,
+ QQ2HQQ_0J = 201,
+ QQ2HQQ_1J = 202,
+ QQ2HQQ_GE2J_MJJ_0_60_PTHJJ_0_25 = 203,
+ QQ2HQQ_GE2J_MJJ_60_120_PTHJJ_0_25 = 204,
+ QQ2HQQ_GE2J_MJJ_120_350_PTHJJ_0_25 = 205,
+ QQ2HQQ_GE2J_MJJ_0_60_PTHJJ_GT25 = 206,
+ QQ2HQQ_GE2J_MJJ_60_120_PTHJJ_GT25 = 207,
+ QQ2HQQ_GE2J_MJJ_120_350_PTHJJ_GT25 = 208,
+ QQ2HQQ_GE2J_MJJ_350_700_PTH_0_200_PTHJJ_0_25 = 209,
+ QQ2HQQ_GE2J_MJJ_350_700_PTH_0_200_PTHJJ_GT25 = 210,
+ QQ2HQQ_GE2J_MJJ_700_1000_PTH_0_200_PTHJJ_0_25 = 211,
+ QQ2HQQ_GE2J_MJJ_700_1000_PTH_0_200_PTHJJ_GT25 = 212,
+ QQ2HQQ_GE2J_MJJ_1000_1500_PTH_0_200_PTHJJ_0_25 = 213,
+ QQ2HQQ_GE2J_MJJ_1000_1500_PTH_0_200_PTHJJ_GT25 = 214,
+ QQ2HQQ_GE2J_MJJ_GT1500_PTH_0_200_PTHJJ_0_25 = 215,
+ QQ2HQQ_GE2J_MJJ_GT1500_PTH_0_200_PTHJJ_GT25 = 216,
+ QQ2HQQ_GE2J_MJJ_350_700_PTH_GT200_PTHJJ_0_25 = 217,
+ QQ2HQQ_GE2J_MJJ_350_700_PTH_GT200_PTHJJ_GT25 = 218,
+ QQ2HQQ_GE2J_MJJ_700_1000_PTH_GT200_PTHJJ_0_25 = 219,
+ QQ2HQQ_GE2J_MJJ_700_1000_PTH_GT200_PTHJJ_GT25 = 220,
+ QQ2HQQ_GE2J_MJJ_1000_1500_PTH_GT200_PTHJJ_0_25 = 221,
+ QQ2HQQ_GE2J_MJJ_1000_1500_PTH_GT200_PTHJJ_GT25 = 222,
+ QQ2HQQ_GE2J_MJJ_GT1500_PTH_GT200_PTHJJ_0_25 = 223,
+ QQ2HQQ_GE2J_MJJ_GT1500_PTH_GT200_PTHJJ_GT25 = 224,
+ // qq -> WH
+ QQ2HLNU_FWDH = 300,
+ QQ2HLNU_PTV_0_75_0J = 301,
+ QQ2HLNU_PTV_75_150_0J = 302,
+ QQ2HLNU_PTV_150_250_0J = 303,
+ QQ2HLNU_PTV_250_400_0J = 304,
+ QQ2HLNU_PTV_GT400_0J = 305,
+ QQ2HLNU_PTV_0_75_1J = 306,
+ QQ2HLNU_PTV_75_150_1J = 307,
+ QQ2HLNU_PTV_150_250_1J = 308,
+ QQ2HLNU_PTV_250_400_1J = 309,
+ QQ2HLNU_PTV_GT400_1J = 310,
+ QQ2HLNU_PTV_0_75_GE2J = 311,
+ QQ2HLNU_PTV_75_150_GE2J = 312,
+ QQ2HLNU_PTV_150_250_GE2J = 313,
+ QQ2HLNU_PTV_250_400_GE2J = 314,
+ QQ2HLNU_PTV_GT400_GE2J = 315,
+ // qq -> ZH
+ QQ2HLL_FWDH = 400,
+ QQ2HLL_PTV_0_75_0J = 401,
+ QQ2HLL_PTV_75_150_0J = 402,
+ QQ2HLL_PTV_150_250_0J = 403,
+ QQ2HLL_PTV_250_400_0J = 404,
+ QQ2HLL_PTV_GT400_0J = 405,
+ QQ2HLL_PTV_0_75_1J = 406,
+ QQ2HLL_PTV_75_150_1J = 407,
+ QQ2HLL_PTV_150_250_1J = 408,
+ QQ2HLL_PTV_250_400_1J = 409,
+ QQ2HLL_PTV_GT400_1J = 410,
+ QQ2HLL_PTV_0_75_GE2J = 411,
+ QQ2HLL_PTV_75_150_GE2J = 412,
+ QQ2HLL_PTV_150_250_GE2J = 413,
+ QQ2HLL_PTV_250_400_GE2J = 414,
+ QQ2HLL_PTV_GT400_GE2J = 415,
+ // gg -> ZH
+ GG2HLL_FWDH = 500,
+ GG2HLL_PTV_0_75_0J = 501,
+ GG2HLL_PTV_75_150_0J = 502,
+ GG2HLL_PTV_150_250_0J = 503,
+ GG2HLL_PTV_250_400_0J = 504,
+ GG2HLL_PTV_GT400_0J = 505,
+ GG2HLL_PTV_0_75_1J = 506,
+ GG2HLL_PTV_75_150_1J = 507,
+ GG2HLL_PTV_150_250_1J = 508,
+ GG2HLL_PTV_250_400_1J = 509,
+ GG2HLL_PTV_GT400_1J = 510,
+ GG2HLL_PTV_0_75_GE2J = 511,
+ GG2HLL_PTV_75_150_GE2J = 512,
+ GG2HLL_PTV_150_250_GE2J = 513,
+ GG2HLL_PTV_250_400_GE2J = 514,
+ GG2HLL_PTV_GT400_GE2J = 515,
+ // ttH
+ TTH_FWDH = 600,
+ TTH_PTH_0_60 = 601,
+ TTH_PTH_60_120 = 602,
+ TTH_PTH_120_200 = 603,
+ TTH_PTH_200_300 = 604,
+ TTH_PTH_300_450 = 605,
+ TTH_PTH_GT450 = 606,
+ // bbH
+ BBH_FWDH = 700, BBH = 701,
+ // tH
+ TH_FWDH = 800, TH = 801
+ };
+ } // namespace Stage1_2_Fine
+
+
+#ifdef ROOT_TLorentzVector
+
+ typedef TLorentzVector TLV;
+ typedef std::vector<TLV> TLVs;
+
+ template <class vec4>
+ TLV MakeTLV(vec4 const p) { return TLV(p.px(),p.py(),p.pz(),p.E()); }
+
+ template <class Vvec4>
+ inline TLVs MakeTLVs(Vvec4 const &rivet_jets){
+ TLVs jets; for ( auto jet:rivet_jets ) jets.push_back(MakeTLV(jet));
+ return jets;
+ }
+
+ // Structure holding information about the current event:
+ // Four-momenta and event classification according to the
+ // Higgs Template Cross Section
+ struct HiggsClassification {
+ // Higgs production mode
+ HTXS::HiggsProdMode prodMode;
+ // The Higgs boson
+ TLV higgs;
+ // The Higgs boson decay products
+ TLV p4decay_higgs;
+ // Associated vector bosons
+ TLV V;
+ // The V-boson decay products
+ TLV p4decay_V;
+ // Jets are built ignoring Higgs decay products and leptons from V decays
+ // jets with pT > 25 GeV and 30 GeV
+ TLVs jets25, jets30;
+ // Event categorization according to YR4 wrtietup
+ // https://cds.cern.ch/record/2138079
+ HTXS::Stage0::Category stage0_cat;
+ HTXS::Stage1::Category stage1_cat_pTjet25GeV;
+ HTXS::Stage1::Category stage1_cat_pTjet30GeV;
+ HTXS::Stage1_2::Category stage1_2_cat_pTjet25GeV;
+ HTXS::Stage1_2::Category stage1_2_cat_pTjet30GeV;
+ HTXS::Stage1_2_Fine::Category stage1_2_fine_cat_pTjet25GeV;
+ HTXS::Stage1_2_Fine::Category stage1_2_fine_cat_pTjet30GeV;
+ // Flag for Z->vv decay mode (needed to split QQ2ZH and GG2ZH)
+ bool isZ2vvDecay;
+ // Error code :: classification was succesful or some error occured
+ HTXS::ErrorCode errorCode;
+ };
+
+ template <class category>
+ inline HiggsClassification* Rivet2Root(category const &htxs_cat_rivet){
+ HTXS::HiggsClassification* cat = new HTXS::HiggsClassification;
+ cat->prodMode = htxs_cat_rivet.prodMode;
+ cat->errorCode = htxs_cat_rivet.errorCode;
+ cat->higgs = MakeTLV(htxs_cat_rivet.higgs);
+ cat->V = MakeTLV(htxs_cat_rivet.V);
+ cat->p4decay_higgs = MakeTLV(htxs_cat_rivet.p4decay_higgs);
+ cat->p4decay_V = MakeTLV(htxs_cat_rivet.p4decay_V);
+ cat->jets25 = MakeTLVs(htxs_cat_rivet.jets25);
+ cat->jets30 = MakeTLVs(htxs_cat_rivet.jets30);
+ cat->stage0_cat = htxs_cat_rivet.stage0_cat;
+ cat->stage1_cat_pTjet25GeV = htxs_cat_rivet.stage1_cat_pTjet25GeV;
+ cat->stage1_cat_pTjet30GeV = htxs_cat_rivet.stage1_cat_pTjet30GeV;
+ cat->stage1_2_cat_pTjet25GeV = htxs_cat_rivet.stage1_2_cat_pTjet25GeV;
+ cat->stage1_2_cat_pTjet30GeV = htxs_cat_rivet.stage1_2_cat_pTjet30GeV;
+ cat->stage1_2_fine_cat_pTjet25GeV = htxs_cat_rivet.stage1_2_fine_cat_pTjet25GeV;
+ cat->stage1_2_fine_cat_pTjet30GeV = htxs_cat_rivet.stage1_2_fine_cat_pTjet30GeV;
+ cat->isZ2vvDecay = htxs_cat_rivet.isZ2vvDecay;
+ return cat;
+ }
+
+
+
+ inline int HTXSstage1_to_HTXSstage1FineIndex(HTXS::Stage1::Category stage1,
+ HiggsProdMode prodMode, tH_type tH) {
+
+ if(stage1==HTXS::Stage1::Category::UNKNOWN) return 0;
+ int P = (int)(stage1 / 100);
+ int F = (int)(stage1 % 100);
+ // 1.a spit tH categories
+ if (prodMode==HiggsProdMode::TH) {
+ // check that tH splitting is valid for Stage-1 FineIndex
+ // else return unknown category
+ if(tH==tH_type::noTH) return 0;
+ // check if forward tH
+ int fwdH = F==0?0:1;
+ return (49 + 2*(tH-1) +fwdH);
+ }
+ // 1.b QQ2HQQ --> split into VBF, WH, ZH -> HQQ
+ // offset vector 1: input is the Higgs prodMode
+ // first two indicies are dummies, given that this is only called for prodMode=2,3,4
+ std::vector<int> pMode_offset = {0,0,13,19,25};
+ if (P==2) return (F + pMode_offset[prodMode]);
+ // 1.c remaining categories
+ // offset vector 2: input is the Stage-1 category P
+ // third index is dummy, given that this is called for category P=0,1,3,4,5,6,7
+ std::vector<int> catP_offset = {0,1,0,31,36,41,45,47};
+ return (F + catP_offset[P]);
+ }
+
+ inline int HTXSstage1_to_HTXSstage1FineIndex(const HiggsClassification &stxs,
+ tH_type tH=noTH, bool jets_pT25 = false) {
+ HTXS::Stage1::Category stage1 =
+ jets_pT25==false?stxs.stage1_cat_pTjet30GeV:
+ stxs.stage1_cat_pTjet25GeV;
+ return HTXSstage1_to_HTXSstage1FineIndex(stage1,stxs.prodMode,tH);
+ }
+
+ inline int HTXSstage1_to_index(HTXS::Stage1::Category stage1) {
+ // the Stage-1 categories
+ int P = (int)(stage1 / 100);
+ int F = (int)(stage1 % 100);
+ std::vector<int> offset{0,1,13,19,24,29,33,35,37,39};
+ // convert to linear values
+ return ( F + offset[P] );
+ }
+
+ //Same for Stage1_2 categories
+ inline int HTXSstage1_2_to_HTXSstage1_2_FineIndex(HTXS::Stage1_2::Category stage1_2,
+ HiggsProdMode prodMode, tH_type tH) {
+
+ if(stage1_2==HTXS::Stage1_2::Category::UNKNOWN) return 0;
+ int P = (int)(stage1_2 / 100);
+ int F = (int)(stage1_2 % 100);
+ // 1.a spit tH categories
+ if (prodMode==HiggsProdMode::TH) {
+ // check that tH splitting is valid for Stage-1 FineIndex
+ // else return unknown category
+ if(tH==tH_type::noTH) return 0;
+ // check if forward tH
+ int fwdH = F==0?0:1;
+ return (94 + 2*(tH-1) +fwdH);
+ }
+ // 1.b QQ2HQQ --> split into VBF, WH, ZH -> HQQ
+ // offset vector 1: input is the Higgs prodMode
+ // first two indicies are dummies, given that this is only called for prodMode=2,3,4
+ std::vector<int> pMode_offset = {0,0,35,46,57};
+ if (P==2) return (F + pMode_offset[prodMode]);
+ // 1.c GG2ZH split into gg->ZH-had and gg->ZH-lep
+ if (prodMode==HiggsProdMode::GG2ZH && P==1) return F + 18;
+ // 1.d remaining categories
+ // offset vector 2: input is the Stage-1 category P
+ // third index is dummy, given that this is called for category P=0,1,3,4,5,6,7
+ std::vector<int> catP_offset = {0,1,0,68,74,80,86,92};
+ return (F + catP_offset[P]);
+ }
+
+ inline int HTXSstage1_2_to_HTXSstage1_2_FineIndex(const HiggsClassification &stxs,
+ tH_type tH=noTH, bool jets_pT25 = false) {
+ HTXS::Stage1_2::Category stage1_2 =
+ jets_pT25==false?stxs.stage1_2_cat_pTjet30GeV:
+ stxs.stage1_2_cat_pTjet25GeV;
+ return HTXSstage1_2_to_HTXSstage1_2_FineIndex(stage1_2,stxs.prodMode,tH);
+ }
+
+ inline int HTXSstage1_2_to_index(HTXS::Stage1_2::Category stage1_2) {
+ // the Stage-1 categories
+ int P = (int)(stage1_2 / 100);
+ int F = (int)(stage1_2 % 100);
+ std::vector<int> offset{0,1,18,29,35,41,47,53,55,57};
+ // convert to linear values
+ return ( F + offset[P] );
+ }
+
+ //Same for Stage1_2_Fine categories
+ inline int HTXSstage1_2_Fine_to_HTXSstage1_2_Fine_FineIndex(HTXS::Stage1_2_Fine::Category Stage1_2_Fine,
+ HiggsProdMode prodMode, tH_type tH) {
+
+ if(Stage1_2_Fine==HTXS::Stage1_2_Fine::Category::UNKNOWN) return 0;
+ int P = (int)(Stage1_2_Fine / 100);
+ int F = (int)(Stage1_2_Fine % 100);
+ // 1.a spit tH categories
+ if (prodMode==HiggsProdMode::TH) {
+ // check that tH splitting is valid for Stage-1 FineIndex
+ // else return unknown category
+ if(tH==tH_type::noTH) return 0;
+ // check if forward tH
+ int fwdH = F==0?0:1;
+ return (189 + 2*(tH-1) +fwdH);
+ }
+ // 1.b QQ2HQQ --> split into VBF, WH, ZH -> HQQ
+ // offset vector 1: input is the Higgs prodMode
+ // first two indicies are dummies, given that this is only called for prodMode=2,3,4
+ std::vector<int> pMode_offset = {0,0,57,82,107};
+ if (P==2) return (F + pMode_offset[prodMode]);
+ // 1.c GG2ZH split into gg->ZH-had and gg->ZH-lep
+ if (prodMode==HiggsProdMode::GG2ZH && P==1) return F + 29;
+ // 1.d remaining categories
+ // offset vector 2: input is the Stage-1 category P
+ // third index is dummy, given that this is called for category P=0,1,3,4,5,6,7
+ std::vector<int> catP_offset = {0,1,0,132,148,164,180,187};
+ return (F + catP_offset[P]);
+ }
+
+ inline int HTXSstage1_2_Fine_to_HTXSstage1_2_Fine_FineIndex(const HiggsClassification &stxs,
+ tH_type tH=noTH, bool jets_pT25 = false) {
+ HTXS::Stage1_2_Fine::Category Stage1_2_Fine =
+ jets_pT25==false?stxs.stage1_2_fine_cat_pTjet30GeV:
+ stxs.stage1_2_fine_cat_pTjet25GeV;
+ return HTXSstage1_2_Fine_to_HTXSstage1_2_Fine_FineIndex(Stage1_2_Fine,stxs.prodMode,tH);
+ }
+
+ inline int HTXSstage1_2_Fine_to_index(HTXS::Stage1_2_Fine::Category Stage1_2_Fine) {
+ // the Stage-1_2_Fine categories
+ int P = (int)(Stage1_2_Fine / 100);
+ int F = (int)(Stage1_2_Fine % 100);
+ std::vector<int> offset{0,1,29,54,70,86,102,109,111,113};
+ // convert to linear values
+ return ( F + offset[P] );
+ }
+
+
+#endif // ROOT_TLorentzVector
+
+} // namespace HTXS
+
+
+#ifdef RIVET_Particle_HH
+
+namespace Rivet {
+
+ /// @struct HiggsClassification
+ /// @brief Structure holding information about the current event:
+ /// Four-momenta and event classification according to the
+ /// Higgs Template Cross Section
+ struct HiggsClassification {
+ /// Higgs production mode
+ HTXS::HiggsProdMode prodMode;
+ /// The Higgs boson
+ Rivet::Particle higgs;
+ /// Vector boson produced in asscoiation with the Higgs
+ Rivet::Particle V;
+ /// The four momentum sum of all stable decay products orignating from the Higgs boson
+ Rivet::FourMomentum p4decay_higgs;
+ /// The four momentum sum of all stable decay products orignating from the vector boson in associated production
+ Rivet::FourMomentum p4decay_V;
+ /// Jets built ignoring Higgs decay products and leptons from V decays, pT thresholds at 25 GeV and 30 GeV
+ Rivet::Jets jets25, jets30;
+ /// Stage-0 HTXS event classifcation, see: https://cds.cern.ch/record/2138079
+ HTXS::Stage0::Category stage0_cat;
+ /// Stage-1 HTXS event classifcation, see: https://cds.cern.ch/record/2138079
+ HTXS::Stage1::Category stage1_cat_pTjet25GeV;
+ /// Stage-1 HTXS event classifcation, see: https://cds.cern.ch/record/2138079
+ HTXS::Stage1::Category stage1_cat_pTjet30GeV;
+ /// Stage-1_2 STXS event classifcation, see: https://twiki.cern.ch/twiki/bin/view/LHCPhysics/LHCHXSWGFiducialAndSTXS#Stage_1_2
+ HTXS::Stage1_2::Category stage1_2_cat_pTjet25GeV;
+ /// Stage-1_2 STXS event classifcation, see: https://twiki.cern.ch/twiki/bin/view/LHCPhysics/LHCHXSWGFiducialAndSTXS#Stage_1_2
+ HTXS::Stage1_2::Category stage1_2_cat_pTjet30GeV;
+ /// Stage-1_2 STXS event classifcation, see: https://twiki.cern.ch/twiki/bin/view/LHCPhysics/LHCHXSWGFiducialAndSTXS#Stage_1_2
+ HTXS::Stage1_2_Fine::Category stage1_2_fine_cat_pTjet25GeV;
+ /// Stage-1_2 STXS event classifcation, see: https://twiki.cern.ch/twiki/bin/view/LHCPhysics/LHCHXSWGFiducialAndSTXS#Stage_1_2
+ HTXS::Stage1_2_Fine::Category stage1_2_fine_cat_pTjet30GeV;
+ /// Flag to distiguish the Z->vv and Z->l+l- decay modes
+ bool isZ2vvDecay;
+ /// Error code: Whether classification was succesful or some error occured
+ HTXS::ErrorCode errorCode;
+ };
+} // namespace Rivet
+#endif // RIVET_Particle_HH
+
+
+
+#endif // TRUTHRIVETTOOLS_HIGGSTEMPLATECROSSSECTIONSDEFS_H
diff --git a/Generators/TruthRivetTools/TruthRivetTools/HiggsTruthCategoryTool.h b/Generators/TruthRivetTools/TruthRivetTools/HiggsTruthCategoryTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..9052309f23b0bc6ee753d9fbd83815f4017ca2dd
--- /dev/null
+++ b/Generators/TruthRivetTools/TruthRivetTools/HiggsTruthCategoryTool.h
@@ -0,0 +1,56 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+/*
+ * Dual-use tool interface for Rivet routine for classifying MC events according to the Higgs template cross section categories
+ * Authors: Jim Lacey (Carleton University)
+ * <james.lacey@cern.ch,jlacey@physics.carleton.ca>
+ */
+
+#ifndef TRUTHRIVETTOOLS_HIGGSTRUTHCATEGORYTOOL_H
+#define TRUTHRIVETTOOLS_HIGGSTRUTHCATEGORYTOOL_H 1
+
+#include "Rivet/AnalysisHandler.hh"
+#include "TruthRivetTools/HiggsTemplateCrossSections.h"
+
+// To avoid coflict of UNUSED macro of
+// Control/CxxUtils/CxxUtils/unused.h and Rivet/Tools/Utils.hh
+#ifdef UNUSED
+#undef UNUSED
+#endif // UNUSED
+
+// Base classes
+#include "AsgTools/AsgTool.h"
+#include "GenInterfaces/IHiggsTruthCategoryTool.h"
+
+// Return type (non-pointer)
+// Note: the Template XSec Defs *depends* on having included
+// the TLorentzVector header *before* it is included -- it
+// uses the include guard from TLorentzVector to decide
+// what is available
+#include "TLorentzVector.h"
+#include "TruthRivetTools/HiggsTemplateCrossSectionsDefs.h"
+
+#include "AtlasHepMC/GenEvent.h"
+
+class HiggsTruthCategoryTool
+: public asg::AsgTool,
+ public virtual IHiggsTruthCategoryTool
+{
+ public:
+ ASG_TOOL_CLASS( HiggsTruthCategoryTool , IHiggsTruthCategoryTool )
+ HiggsTruthCategoryTool( const std::string& name );
+ ~HiggsTruthCategoryTool() { };
+ public:
+ Rivet::AnalysisHandler *rivetAnaHandler; //!
+ Rivet::HiggsTemplateCrossSections *higgsTemplateCrossSections; //!
+ virtual StatusCode initialize() override;
+ StatusCode finalize () override;
+ HTXS::HiggsClassification* getHiggsTruthCategoryObject(const HepMC::GenEvent& HepMCEvent, const HTXS::HiggsProdMode prodMode) override;
+ private:
+ bool m_isInitialized;
+ bool m_outHistos;
+};
+
+#endif //> !HIGGSTRUTHCLASSIFIER_HIGGSTRUTHCATEGORYTOOL_H
diff --git a/Generators/TruthRivetTools/TruthRivetTools/TruthRivetToolsDict.h b/Generators/TruthRivetTools/TruthRivetTools/TruthRivetToolsDict.h
new file mode 100644
index 0000000000000000000000000000000000000000..462805a379e7c2fcc56ee84d414d7780b17d5dab
--- /dev/null
+++ b/Generators/TruthRivetTools/TruthRivetTools/TruthRivetToolsDict.h
@@ -0,0 +1,11 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// $Id: TruthRivetToolsDict.h 769074 2016-08-22 10:04:26Z krasznaa $
+#ifndef TRUTHRIVETTOOLS_TRUTHRIVETTOOLSDICT_H
+#define TRUTHRIVETTOOLS_TRUTHRIVETTOOLSDICT_H
+
+#include "TruthRivetTools/HiggsTruthCategoryTool.h"
+
+#endif // TRUTHRIVETTOOLS_TRUTHRIVETTOOLSDICT_H
diff --git a/Generators/TruthRivetTools/TruthRivetTools/selection.xml b/Generators/TruthRivetTools/TruthRivetTools/selection.xml
new file mode 100644
index 0000000000000000000000000000000000000000..1068f79214da02683378b29eb62fa74ae58795d3
--- /dev/null
+++ b/Generators/TruthRivetTools/TruthRivetTools/selection.xml
@@ -0,0 +1,4 @@
+<lcgdict>
+ <class name="HiggsTruthCategoryTool" />
+ <class name="IHiggsTruthCategoryTool" />
+</lcgdict>
diff --git a/Generators/TruthRivetTools/share/packages b/Generators/TruthRivetTools/share/packages
new file mode 100644
index 0000000000000000000000000000000000000000..7d04d93db20b23116cc406dce0d1ca491ecba53d
--- /dev/null
+++ b/Generators/TruthRivetTools/share/packages
@@ -0,0 +1,6 @@
+svn+ssh://svn.cern.ch/reps/atlasinst/Institutes/Carleton/HiggsTruthClassifier/BoostExternal/trunk
+svn+ssh://svn.cern.ch/reps/atlasinst/Institutes/Carleton/HiggsTruthClassifier/GSLExternal/trunk
+svn+ssh://svn.cern.ch/reps/atlasinst/Institutes/Carleton/HiggsTruthClassifier/HepMCExternal/trunk
+svn+ssh://svn.cern.ch/reps/atlasinst/Institutes/Carleton/HiggsTruthClassifier/FastJetExternal/trunk
+svn+ssh://svn.cern.ch/reps/atlasinst/Institutes/Carleton/HiggsTruthClassifier/YODAExternal/trunk
+svn+ssh://svn.cern.ch/reps/atlasinst/Institutes/Carleton/HiggsTruthClassifier/RivetExternal/trunk
diff --git a/Generators/TruthRivetTools/src/components/TruthRivetTools_entries.cxx b/Generators/TruthRivetTools/src/components/TruthRivetTools_entries.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..669dd2a7bb8aba977cac94a4b9ad86c627bf771e
--- /dev/null
+++ b/Generators/TruthRivetTools/src/components/TruthRivetTools_entries.cxx
@@ -0,0 +1,3 @@
+#include "TruthRivetTools/HiggsTruthCategoryTool.h"
+
+DECLARE_COMPONENT( HiggsTruthCategoryTool )
\ No newline at end of file