diff --git a/Generators/GenInterfaces/GenInterfaces/IHiggsTruthCategoryTool.h b/Generators/GenInterfaces/GenInterfaces/IHiggsTruthCategoryTool.h index dd8b1e0de73c0940884f400ac840c8c78f177623..52dce88eea096773f0ee2834f692466f7f45c293 100644 --- a/Generators/GenInterfaces/GenInterfaces/IHiggsTruthCategoryTool.h +++ b/Generators/GenInterfaces/GenInterfaces/IHiggsTruthCategoryTool.h @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration */ /* @@ -13,18 +13,17 @@ #include "AsgTools/IAsgTool.h" #include "AtlasHepMC/GenEvent.h" - +#include "TruthRivetTools/HiggsTemplateCrossSectionsDefs.h" namespace HTXS { struct 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 ~IHiggsTruthCategoryTool() = default; + public: virtual HTXS::HiggsClassification* getHiggsTruthCategoryObject(const HepMC::GenEvent& HepMCEvent, const HTXS::HiggsProdMode prodMode) const =0; }; diff --git a/Generators/TruthConverters/CMakeLists.txt b/Generators/TruthConverters/CMakeLists.txt index 90155d36f1ec0aab37afed9f9df164d805e66bb8..14e0bedec85a5d6f17ee438aa5f4061031295741 100644 --- a/Generators/TruthConverters/CMakeLists.txt +++ b/Generators/TruthConverters/CMakeLists.txt @@ -21,3 +21,5 @@ atlas_add_component( TruthConverters atlas_add_dictionary( TruthConvertersDict TruthConverters/xAODtoHepMCDict.h TruthConverters/selection.xml LINK_LIBRARIES TruthConvertersLib ) + +atlas_install_python_modules( python/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} ) \ No newline at end of file diff --git a/Generators/TruthConverters/Root/xAODtoHepMCTool.cxx b/Generators/TruthConverters/Root/xAODtoHepMCTool.cxx index 2caacc61a223d0ab11ef8ed57f375c2f52876fef..3099309437d79590f4a0b7f0e34223384b280eba 100644 --- a/Generators/TruthConverters/Root/xAODtoHepMCTool.cxx +++ b/Generators/TruthConverters/Root/xAODtoHepMCTool.cxx @@ -2,61 +2,67 @@ Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration */ - #include "TruthConverters/xAODtoHepMCTool.h" -xAODtoHepMCTool::xAODtoHepMCTool( const std::string& name ) - : asg::AsgTool( name ), - m_momFac(1.), //<-- MeV/GeV conversion factor - m_lenFac(1.), //<-- length conversion factor - m_signalOnly(true), - m_maxCount(0) +xAODtoHepMCTool::xAODtoHepMCTool(const std::string &name) + : asg::AsgTool(name), + m_momFac(1.), //<-- MeV/GeV conversion factor + m_lenFac(1.), //<-- length conversion factor + m_signalOnly(true), + m_maxCount(0) { } - -StatusCode xAODtoHepMCTool::initialize() { - ATH_MSG_INFO ("Initializing xAODtoHepMCTool "<< name() << "..."); - ATH_MSG_INFO ("SignalOnly = " << m_signalOnly); - m_evtCount = 0; +StatusCode xAODtoHepMCTool::initialize() +{ + ATH_MSG_INFO("Initializing xAODtoHepMCTool " << name() << "..."); + ATH_MSG_INFO("SignalOnly = " << m_signalOnly); + m_evtCount = 0; m_badSuggest = 0; - m_noProdVtx = 0; - m_badBeams = 0; + m_noProdVtx = 0; + m_badBeams = 0; return StatusCode::SUCCESS; } - - -StatusCode xAODtoHepMCTool :: finalize () +StatusCode xAODtoHepMCTool ::finalize() { - ATH_MSG_INFO ("=============================================================="); - ATH_MSG_INFO ("========== xAOD -> HepMC Tool :: Run Summary =========="); - ATH_MSG_INFO ("=============================================================="); - if( m_badSuggest ){ - ATH_MSG_INFO ("Number of suggest_barcode failures = " << m_badSuggest); - ATH_MSG_INFO (" ... check input particle/vertex barcodes."); - } else { - ATH_MSG_INFO ("No suggest_barcode failures"); + ATH_MSG_INFO("=============================================================="); + ATH_MSG_INFO("========== xAOD -> HepMC Tool :: Run Summary =========="); + ATH_MSG_INFO("=============================================================="); + if (m_badSuggest) + { + ATH_MSG_INFO("Number of suggest_barcode failures = " << m_badSuggest); + ATH_MSG_INFO(" ... check input particle/vertex barcodes."); + } + else + { + ATH_MSG_INFO("No suggest_barcode failures"); + } + if (m_noProdVtx) + { + ATH_MSG_INFO("Number of events that have a missing production vertex = " << m_noProdVtx); + } + else + { + ATH_MSG_INFO("No missing production vertices"); } - if ( m_noProdVtx ) { - ATH_MSG_INFO ("Number of events that have a missing production vertex = " << m_noProdVtx); - }else{ - ATH_MSG_INFO ("No missing production vertices"); + if (m_badBeams) + { + ATH_MSG_INFO("Number events with improperly defined beamparticles = " << m_badBeams); + ATH_MSG_INFO("Inconsistencies with the beam particles storage in the event record!"); + ATH_MSG_INFO("... check xAOD::TruthEvent record ..."); } - if ( m_badBeams ){ - ATH_MSG_INFO ("Number events with improperly defined beamparticles = " << m_badBeams); - ATH_MSG_INFO ("Inconsistencies with the beam particles storage in the event record!"); - ATH_MSG_INFO ("... check xAOD::TruthEvent record ..."); - }else{ - ATH_MSG_INFO ("No events with undefined beams."); + else + { + ATH_MSG_INFO("No events with undefined beams."); } - ATH_MSG_INFO ("=============================================================="); - ATH_MSG_INFO ("=================== End Run Summary ==================="); - ATH_MSG_INFO ("=============================================================="); + ATH_MSG_INFO("=============================================================="); + ATH_MSG_INFO("=================== End Run Summary ==================="); + ATH_MSG_INFO("=============================================================="); return StatusCode::SUCCESS; } -std::vector xAODtoHepMCTool :: getHepMCEvents(const xAOD::TruthEventContainer* xTruthEventContainer, const xAOD::EventInfo* eventInfo) const +std::vector xAODtoHepMCTool ::getHepMCEvents(const xAOD::TruthEventContainer *xTruthEventContainer, const xAOD::EventInfo *eventInfo) const { ++m_evtCount; bool doPrint = m_evtCount < m_maxCount; @@ -65,75 +71,89 @@ std::vector xAODtoHepMCTool :: getHepMCEvents(const xAOD::Truth // Loop over xAOD truth container // Signal event is always first, followed by pileup events ATH_MSG_DEBUG("Start event loop"); - for (const xAOD::TruthEvent* xAODEvent : *xTruthEventContainer) { - if( doPrint ) ATH_MSG_DEBUG("XXX Printing xAOD Event"); - if( doPrint ) printxAODEvent( xAODEvent,eventInfo ); + for (const xAOD::TruthEvent *xAODEvent : *xTruthEventContainer) + { + if (doPrint) + ATH_MSG_DEBUG("XXX Printing xAOD Event"); + if (doPrint) + printxAODEvent(xAODEvent, eventInfo); // Create GenEvent for each xAOD truth event ATH_MSG_DEBUG("Create new GenEvent"); - HepMC::GenEvent hepmcEvent = createHepMCEvent(xAODEvent,eventInfo); + HepMC::GenEvent hepmcEvent = createHepMCEvent(xAODEvent, eventInfo); // Insert into McEventCollection - mcEventCollection.push_back(hepmcEvent); - if( doPrint ) ATH_MSG_DEBUG("XXX Printing HepMC Event"); - if( doPrint ) HepMC::Print::line(std::cout,hepmcEvent); + mcEventCollection.push_back(std::move(hepmcEvent)); + if (doPrint) + ATH_MSG_DEBUG("XXX Printing HepMC Event"); + if (doPrint) + HepMC::Print::line(std::cout, hepmcEvent); // Quit if signal only - if( m_signalOnly ) break; - } + if (m_signalOnly) + break; + } return mcEventCollection; } +HepMC::GenEvent xAODtoHepMCTool::createHepMCEvent(const xAOD::TruthEvent *xEvt, const xAOD::EventInfo *eventInfo) const +{ -HepMC::GenEvent xAODtoHepMCTool::createHepMCEvent(const xAOD::TruthEvent* xEvt, const xAOD::EventInfo* eventInfo) const { - - ///EVENT LEVEL + /// EVENT LEVEL HepMC::GenEvent genEvt; - long long int evtNum = eventInfo->eventNumber(); + long long int evtNum = eventInfo->eventNumber(); genEvt.set_event_number(evtNum); - ATH_MSG_DEBUG("Start createHepMCEvent for event " < vertexMap; + std::map vertexMap; // Loop over all of the particles in the event, call particle builder - // Call suggest_barcode only after insertion! - for (auto tlink:xEvt->truthParticleLinks()) { - if (!tlink.isValid()) continue; - const xAOD::TruthParticle* xPart = *tlink; - + // Call suggest_barcode only after insertion! + for (auto tlink : xEvt->truthParticleLinks()) + { + if (!tlink.isValid()) + continue; + const xAOD::TruthParticle *xPart = *tlink; + // sanity check - if (xPart == nullptr) { + if (xPart == nullptr) + { ATH_MSG_WARNING("xAOD TruthParticle is equal to NULL. This should not happen!"); continue; } - if( !xPart->hasProdVtx() && !xPart->hasDecayVtx() ){ - ATH_MSG_WARNING("xAOD particle with no vertices, bc = "<barcode()); + if (!xPart->hasProdVtx() && !xPart->hasDecayVtx()) + { + ATH_MSG_WARNING("xAOD particle with no vertices, bc = " << xPart->barcode()); continue; } - // skip particles with barcode > 200000 --> Geant4 secondaries - if ( xPart->barcode() > 200000 ) continue; + // skip particles with barcode > 200000 --> Geant4 secondaries + if (xPart->barcode() > 200000) + continue; - // Create GenParticle - //presumably the GenEvent takes ownership of this, but creating a unique_ptr here as that will only happen if there's an associated vertex + // Create GenParticle + // presumably the GenEvent takes ownership of this, but creating a unique_ptr here as that will only happen if there's an associated vertex #ifdef HEPMC3 - auto hepmcParticle=createHepMCParticle(xPart) ; + auto hepmcParticle = createHepMCParticle(xPart); #else - std::unique_ptr hepmcParticle( createHepMCParticle(xPart) ); + std::unique_ptr hepmcParticle(createHepMCParticle(xPart)); #endif int bcpart = xPart->barcode(); - + // status 10902 should be treated just as status 2 - if ( hepmcParticle->status() == 10902 ) hepmcParticle->set_status(2); + if (hepmcParticle->status() == 10902) + hepmcParticle->set_status(2); // Get the production and decay vertices - if( xPart->hasProdVtx() ) { - const xAOD::TruthVertex* xAODProdVtx = xPart->prodVtx(); - // skip production vertices with barcode > 200000 --> Geant4 secondaries - if ( std::abs(xAODProdVtx->barcode()) > 200000 ) continue; + if (xPart->hasProdVtx()) + { + const xAOD::TruthVertex *xAODProdVtx = xPart->prodVtx(); + // skip production vertices with barcode > 200000 --> Geant4 secondaries + if (std::abs(xAODProdVtx->barcode()) > 200000) + continue; bool prodVtxSeenBefore(false); // is this new? - auto hepmcProdVtx = vertexHelper(xAODProdVtx,vertexMap,prodVtxSeenBefore); + auto hepmcProdVtx = vertexHelper(xAODProdVtx, vertexMap, prodVtxSeenBefore); // Set the decay/production links #ifdef HEPMC3 hepmcProdVtx->add_particle_out(hepmcParticle); @@ -141,28 +161,42 @@ HepMC::GenEvent xAODtoHepMCTool::createHepMCEvent(const xAOD::TruthEvent* xEvt, hepmcProdVtx->add_particle_out(hepmcParticle.get()); #endif // Insert into Event - if (!prodVtxSeenBefore){ - genEvt.add_vertex(hepmcProdVtx); - if( !HepMC::suggest_barcode(hepmcProdVtx,xAODProdVtx->barcode()) ){ - ATH_MSG_WARNING("suggest_barcode failed for vertex "<barcode()); - ++m_badSuggest; - } - } - if( !HepMC::suggest_barcode(hepmcParticle,bcpart) ){ - ATH_MSG_DEBUG("suggest_barcode failed for particle " <barcode())) + { + ATH_MSG_WARNING("suggest_barcode failed for vertex " << xAODProdVtx->barcode()); + ++m_badSuggest; + } + } + if (!HepMC::suggest_barcode(hepmcParticle, bcpart)) + { + ATH_MSG_VERBOSE("suggest_barcode failed for particle " << bcpart); + ++m_badSuggest; } bcpart = 0; - } else { - ATH_MSG_DEBUG("No production vertex found for particle "<barcode()); } - - if( xPart->hasDecayVtx() ){ - const xAOD::TruthVertex* xAODDecayVtx = xPart->decayVtx(); - // skip decay vertices with barcode > 200000 --> Geant4 secondaries - if ( std::abs(xAODDecayVtx->barcode()) > 200000 ) continue; + else + { + ATH_MSG_VERBOSE("No production vertex found for particle " << xPart->barcode()); + } + + if (xPart->hasDecayVtx()) + { + const xAOD::TruthVertex *xAODDecayVtx = xPart->decayVtx(); + // skip decay vertices with barcode > 200000 --> Geant4 secondaries + if (std::abs(xAODDecayVtx->barcode()) > 200000) + { +/// Avoid double deletion +#ifndef HEPMC3 + if (xPart->hasProdVtx()) + hepmcParticle.release(); +#endif + continue; + } bool decayVtxSeenBefore(false); // is this new? - auto hepmcDecayVtx = vertexHelper(xAODDecayVtx,vertexMap,decayVtxSeenBefore); + auto hepmcDecayVtx = vertexHelper(xAODDecayVtx, vertexMap, decayVtxSeenBefore); // Set the decay/production links #ifdef HEPMC3 hepmcDecayVtx->add_particle_in(hepmcParticle); @@ -170,20 +204,24 @@ HepMC::GenEvent xAODtoHepMCTool::createHepMCEvent(const xAOD::TruthEvent* xEvt, hepmcDecayVtx->add_particle_in(hepmcParticle.get()); #endif // Insert into Event - if (!decayVtxSeenBefore){ - genEvt.add_vertex(hepmcDecayVtx); - if( !HepMC::suggest_barcode(hepmcDecayVtx,xAODDecayVtx->barcode()) ){ - ATH_MSG_WARNING("suggest_barcode failed for vertex " - <barcode()); - ++m_badSuggest; - } + if (!decayVtxSeenBefore) + { + genEvt.add_vertex(hepmcDecayVtx); + if (!HepMC::suggest_barcode(hepmcDecayVtx, xAODDecayVtx->barcode())) + { + ATH_MSG_WARNING("suggest_barcode failed for vertex " + << xAODDecayVtx->barcode()); + ++m_badSuggest; + } } - if( bcpart != 0 ){ - if( !HepMC::suggest_barcode(hepmcParticle,bcpart) ){ - ATH_MSG_DEBUG("suggest_barcode failed for particle " < &vertexMap, - bool &seenBefore) const { - +HepMC::GenVertexPtr xAODtoHepMCTool::vertexHelper(const xAOD::TruthVertex *xaodVertex, + std::map &vertexMap, + bool &seenBefore) const +{ + HepMC::GenVertexPtr hepmcVertex; - std::map::iterator vMapItr; - vMapItr=vertexMap.find(xaodVertex); + std::map::iterator vMapItr; + vMapItr = vertexMap.find(xaodVertex); // Vertex seen before? - if (vMapItr!=vertexMap.end()) { - // YES: use the HepMC::Vertex already in the map + if (vMapItr != vertexMap.end()) + { + // YES: use the HepMC::Vertex already in the map hepmcVertex = (*vMapItr).second; seenBefore = true; - } else { + } + else + { // NO: create a new HepMC::Vertex vertexMap[xaodVertex] = createHepMCVertex(xaodVertex); - hepmcVertex=vertexMap[xaodVertex]; + hepmcVertex = vertexMap[xaodVertex]; seenBefore = false; } return hepmcVertex; - } // Create the HepMC GenParticle // Call suggest_barcode after insertion! -HepMC::GenParticlePtr xAODtoHepMCTool::createHepMCParticle(const xAOD::TruthParticle* particle) const { - ATH_MSG_VERBOSE("Creating GenParticle for barcode " <barcode()); - const HepMC::FourVector fourVec( m_momFac * particle->px(), m_momFac * particle->py(), m_momFac * particle->pz(), m_momFac * particle->e() ); - auto hepmcParticle=HepMC::newGenParticlePtr(fourVec, particle->pdgId(), particle->status()); - hepmcParticle->set_generated_mass( m_momFac * particle->m()); +HepMC::GenParticlePtr xAODtoHepMCTool::createHepMCParticle(const xAOD::TruthParticle *particle) const +{ + ATH_MSG_VERBOSE("Creating GenParticle for barcode " << particle->barcode()); + const HepMC::FourVector fourVec(m_momFac * particle->px(), m_momFac * particle->py(), m_momFac * particle->pz(), m_momFac * particle->e()); + auto hepmcParticle = HepMC::newGenParticlePtr(fourVec, particle->pdgId(), particle->status()); + hepmcParticle->set_generated_mass(m_momFac * particle->m()); return hepmcParticle; } // Create the HepMC GenVertex // Call suggest_barcode after insertion! -HepMC::GenVertexPtr xAODtoHepMCTool::createHepMCVertex(const xAOD::TruthVertex* vertex) const { - ATH_MSG_VERBOSE("Creating GenVertex for barcode " <barcode()); - HepMC::FourVector prod_pos( m_lenFac * vertex->x(), m_lenFac * vertex->y(),m_lenFac * vertex->z(), m_lenFac * vertex->t() ); - auto genVertex=HepMC::newGenVertexPtr(prod_pos); +HepMC::GenVertexPtr xAODtoHepMCTool::createHepMCVertex(const xAOD::TruthVertex *vertex) const +{ + ATH_MSG_VERBOSE("Creating GenVertex for barcode " << vertex->barcode()); + HepMC::FourVector prod_pos(m_lenFac * vertex->x(), m_lenFac * vertex->y(), m_lenFac * vertex->z(), m_lenFac * vertex->t()); + auto genVertex = HepMC::newGenVertexPtr(prod_pos); return genVertex; } // Print xAODTruth Event. The printout is particle oriented, unlike the // HepMC particle/vertex printout. Geant and pileup particles with // barcode>100000 are omitted. -void xAODtoHepMCTool::printxAODEvent(const xAOD::TruthEvent* event, const xAOD::EventInfo* eventInfo) const{ - +void xAODtoHepMCTool::printxAODEvent(const xAOD::TruthEvent *event, const xAOD::EventInfo *eventInfo) const +{ + std::vector bcPars; std::vector bcKids; long long int evtNum = eventInfo->eventNumber(); - std::cout <<"======================================================================================" <nTruthParticles(); - for(int i=0; itruthParticle(i); - if (part==nullptr) continue; + for (int i = 0; i < nPart; ++i) + { + const xAOD::TruthParticle *part = event->truthParticle(i); + if (part == nullptr) + continue; int bc = part->barcode(); - if( bc > 100000 ) continue; + if (bc > 100000) + continue; int id = part->pdgId(); - if ( id != 25 ) continue; + if (id != 25) + continue; int stat = part->status(); - float px = part->px()/1000.; - float py = part->py()/1000.; - float pz = part->pz()/1000.; - float e = part->e()/1000.; + float px = part->px() / 1000.; + float py = part->py() / 1000.; + float pz = part->pz() / 1000.; + float e = part->e() / 1000.; bcPars.clear(); bcKids.clear(); - if( part->hasProdVtx() ){ - const xAOD::TruthVertex* pvtx = part->prodVtx(); - if( pvtx ) bcPars.push_back(pvtx->barcode()); + if (part->hasProdVtx()) + { + const xAOD::TruthVertex *pvtx = part->prodVtx(); + if (pvtx) + bcPars.push_back(pvtx->barcode()); } - if( part->hasDecayVtx() ){ - const xAOD::TruthVertex* dvtx = part->decayVtx(); - if( dvtx ) bcKids.push_back(dvtx->barcode()); + if (part->hasDecayVtx()) + { + const xAOD::TruthVertex *dvtx = part->decayVtx(); + if (dvtx) + bcKids.push_back(dvtx->barcode()); } - std::cout <addAnalysis(&(*higgsTemplateCrossSections)); + rivetAnaHandler->addAnalysis(higgsTemplateCrossSections); return StatusCode::SUCCESS; } diff --git a/Generators/TruthRivetTools/python/TruthRivetToolsCfg.py b/Generators/TruthRivetTools/python/TruthRivetToolsCfg.py new file mode 100644 index 0000000000000000000000000000000000000000..8bb2ad6a0a667fb4a3ed70a36651f857c277ca0e --- /dev/null +++ b/Generators/TruthRivetTools/python/TruthRivetToolsCfg.py @@ -0,0 +1,8 @@ +# Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration +from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator +from AthenaConfiguration.ComponentFactory import CompFactory + +def HiggsTruthCategoryToolCfg(flags, name="HiggsTruthCategoryTool", **kwargs): + result = ComponentAccumulator() + result.setPrivateTools(CompFactory.HiggsTruthCategoryTool(name,**kwargs)) + return result \ No newline at end of file diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/CMakeLists.txt b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/CMakeLists.txt index 42bbbe30c58fb73e0527b8c5ce16acaefff66edd..e7affbd8f35c8f218b24343d0c987d5a758ab8d7 100644 --- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/CMakeLists.txt +++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/CMakeLists.txt @@ -17,7 +17,7 @@ atlas_add_component( DerivationFrameworkHiggs xAODEventInfo xAODEgamma xAODJet xAODMuon xAODTracking xAODTruth xAODTrigger xAODPFlow FourMomUtils TrkVKalVrtFitterLib TrkVertexFitterInterfaces TrkVertexAnalysisUtilsLib TrigDecisionToolLib VxVertex GammaORToolsLib egammaInterfacesLib PFlowUtilsLib PhotonVertexSelectionLib VrtSecInclusiveLib PFlowUtilsLib - EgammaAnalysisInterfacesLib MuonAnalysisInterfacesLib DerivationFrameworkInterfaces CaloDetDescrLib + EgammaAnalysisInterfacesLib MuonAnalysisInterfacesLib DerivationFrameworkInterfaces CaloDetDescrLib MuonCondSvcLib ExpressionEvaluationLib AthContainers PathResolver TruthConvertersLib TruthRivetToolsLib TruthUtils ElectronPhotonSelectorToolsLib PFlowUtilsLib JpsiUpsilonToolsLib) diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/DerivationFrameworkHiggs/TruthCategoriesDecorator.h b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/DerivationFrameworkHiggs/TruthCategoriesDecorator.h index 6fabb2dfe8ac7f431cc9948db9df0cf420e19b02..13881ddd0e60fe238f010060cec393eecd7a3437 100644 --- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/DerivationFrameworkHiggs/TruthCategoriesDecorator.h +++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/DerivationFrameworkHiggs/TruthCategoriesDecorator.h @@ -1,76 +1,94 @@ /* - Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration */ -// Tool to decorate the EventInfo object with truth categories informations +// Tool to decorate the EventInfo object with truth categories informations // Authors: T. Guillemin, J. Lacey, D. Gillberg - + #ifndef DerivationFrameworkHiggs_TruthCategoriesDecorator_H #define DerivationFrameworkHiggs_TruthCategoriesDecorator_H -#include -#include -#include -#include -#include - -#include "AthenaBaseComps/AthAlgTool.h" +#include "AthenaBaseComps/AthReentrantAlgorithm.h" #include "GaudiKernel/ToolHandle.h" -#include "DerivationFrameworkInterfaces/IAugmentationTool.h" +#include "StoreGate/ReadHandleKey.h" +#include "StoreGate/WriteDecorHandleKeyArray.h" // EDM: typedefs, so must be included and not forward-referenced -#include "xAODTruth/TruthParticleContainer.h" #include "xAODEventInfo/EventInfo.h" +#include "xAODTruth/TruthEventContainer.h" // Note: must include TLorentzVector before the next one -#include "TLorentzVector.h" -#include "TruthRivetTools/HiggsTemplateCrossSectionsDefs.h" -class IHiggsTruthCategoryTool; -class IxAODtoHepMCTool; +#include "GenInterfaces/IHiggsTruthCategoryTool.h" +#include "GenInterfaces/IxAODtoHepMCTool.h" +#include "TruthRivetTools/HiggsTemplateCrossSectionsDefs.h" namespace DerivationFramework { - - class TruthCategoriesDecorator : public AthAlgTool, public IAugmentationTool { - - public: - TruthCategoriesDecorator(const std::string& t, const std::string& n, const IInterface* p); - ~TruthCategoriesDecorator(); - StatusCode initialize(); - StatusCode finalize(); - virtual StatusCode addBranches() const; - - private: - ToolHandle m_xAODtoHepMCTool; - ToolHandle m_higgsTruthCatTool; - - // Path to TEnv file containing MC-channel-numbre <-> HiggsProdMode map - TEnv *m_config; - std::string m_configPath; - - // Detail level. Steers amount of decoration. - // 0: basic information. Categoization ints, Higgs prod mode, Njets, Higgs pT - // 1: the above + Higgs boson 4-vec + associated V 4-vec - // 2: the above + truth jets built excluding the Higgs boson decay - // 3: the above + 4-vector sum of all decay products from Higgs boson and V-boson - int m_detailLevel; - - // methods to locate the TEnv input file - bool fileExists(TString fileName) { return gSystem->AccessPathName(fileName) == false; } - - // Converts a string to a vector of integers - static std::vector vectorize(const TString& string, const TString& sep=" ") ; - - // Method to access the production mode for a given MC channel number - HTXS::HiggsProdMode getHiggsProductionMode(uint32_t mc_channel_number,HTXS::tH_type &th_type) const; - - // Methods for decoration of four vectors - static void decorateFourVec ( const xAOD::EventInfo *eventInfo, const TString& prefix, const TLorentzVector& p4 ) ; - static void decorateFourVecs ( const xAOD::EventInfo *eventInfo, const TString& prefix, const std::vector& p4s ) ; - - - }; /// class - -} /// namespace + + class TruthCategoriesDecorator : public AthReentrantAlgorithm { + public: + TruthCategoriesDecorator(const std::string& n, ISvcLocator* p); + virtual ~TruthCategoriesDecorator() = default; + StatusCode initialize(); + StatusCode execute(const EventContext& ctx) const; + + private: + ToolHandle m_xAODtoHepMCTool{this, "HepMCTool", "xAODtoHepMCTool"}; + ToolHandle m_higgsTruthCatTool{this, "CategoryTool", "HiggsTruthCategoryTool"}; + + // Path to TEnv file containing MC-channel-numbre <-> HiggsProdMode map + Gaudi::Property m_configPath{this, "ConfigPath", "DerivationFrameworkHiggs/HiggsMCsamples.cfg"}; + struct HTXSSample { + /// Higgs production modes, corresponding to input sample + HTXS::HiggsProdMode prod{HTXS::HiggsProdMode::UNKNOWN}; + /// Additional identifier flag for TH production modes + HTXS::tH_type th_type{HTXS::tH_type::noTH}; + // DSIDs satisfying this production mode + std::set dsids{}; + }; + + std::vector m_htxs_samples{}; + + // Detail level. Steers amount of decoration. + // 0: basic information. Categoization ints, Higgs prod mode, Njets, Higgs pT + // 1: the above + Higgs boson 4-vec + associated V 4-vec + // 2: the above + truth jets built excluding the Higgs boson decay + // 3: the above + 4-vector sum of all decay products from Higgs boson and V-boson + Gaudi::Property m_detailLevel{this, "DetailLevel", 3}; + + // Methods for decoration of four vectors + StatusCode decorateFourVec(const EventContext& ctx, const std::string& prefix, const TLorentzVector& p4) const; + StatusCode decorateFourVecs(const EventContext& ctx, const std::string& prefix, const std::vector& p4s) const; + + SG::ReadHandleKey m_truthEvtKey{this, "TruthEvtKey", "TruthEvents"}; + SG::ReadHandleKey m_evtInfoKey{this, "EvtKey", "EventInfo"}; + + using EvtInfoDecorKey = SG::WriteDecorHandleKey; + SG::WriteDecorHandleKeyArray m_STXSDecors{this, "ToolDecorations", {}, "List of all keys decorated by the alg"}; + + /// Set of DecorHandleKeys to write the four momenta needed for the HTXS categorization. + struct FourMomDecoration { + FourMomDecoration(const SG::ReadHandleKey& ev_key, const std::string& prefix) : + pt{ev_key.key() + "." + prefix + "_pt"}, + eta{ev_key.key() + "." + prefix + "_eta"}, + phi{ev_key.key() + "." + prefix + "_phi"}, + m{ev_key.key() + "." + prefix + "_m"} {} + EvtInfoDecorKey pt; + EvtInfoDecorKey eta; + EvtInfoDecorKey phi; + EvtInfoDecorKey m; + std::vector vect() const { return {pt, eta, phi, m}; } + StatusCode initialize(){ + if (!pt.initialize().isSuccess() || !eta.initialize().isSuccess() || !phi.initialize().isSuccess() || !m.initialize().isSuccess()){ + return StatusCode::FAILURE; + } + return StatusCode::SUCCESS; + } + }; + using P4DecorMap = std::map; + P4DecorMap m_p4_decors{}; + }; /// class + +} // namespace DerivationFramework #endif diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/TruthCategories.py b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/TruthCategories.py index 4d01b47e11f9ee98e04a0270efab450857df2266..2bf27e20fa9e222a6248202077a8963355c459de 100644 --- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/TruthCategories.py +++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/TruthCategories.py @@ -1,17 +1,9 @@ -# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration +# Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration from DerivationFrameworkCore.DerivationFrameworkMaster import DerivationFrameworkJob from AthenaCommon.GlobalFlags import globalflags -if globalflags.DataSource()=='geant4': - from AthenaCommon.AppMgr import ToolSvc - import AthenaCommon.CfgMgr as CfgMgr - +if globalflags.DataSource()=='geant4': from DerivationFrameworkHiggs.DerivationFrameworkHiggsConf import DerivationFramework__TruthCategoriesDecorator DFHTXSdecorator = DerivationFramework__TruthCategoriesDecorator(name = "DFHTXSdecorator") - ToolSvc += DFHTXSdecorator - - ToolSvc += DFHTXSdecorator - DerivationFrameworkJob += CfgMgr.DerivationFramework__CommonAugmentation("TruthCategoriesCommonKernel", - AugmentationTools = [DFHTXSdecorator] - ) + DerivationFrameworkJob += DFHTXSdecorator \ No newline at end of file diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/TruthCategoriesConfig.py b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/TruthCategoriesConfig.py new file mode 100644 index 0000000000000000000000000000000000000000..d8043d52ddfd829e2ef28767e472d21e7ba9f2ba --- /dev/null +++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/TruthCategoriesConfig.py @@ -0,0 +1,39 @@ +# Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration +from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator +from AthenaConfiguration.ComponentFactory import CompFactory + +def TruthCategoriesDecoratorCfg(flags, name="TruthCategoriesDecorator", **kwargs): + result = ComponentAccumulator() + if not flags.Input.isMC: return result + from TruthConverters.TruthConvertersCfg import xAODtoHEPToolCfg + from TruthRivetTools.TruthRivetToolsCfg import HiggsTruthCategoryToolCfg + kwargs.setdefault("HepMCTool", result.popToolsAndMerge(xAODtoHEPToolCfg(flags))) + kwargs.setdefault("CategoryTool", result.popToolsAndMerge(HiggsTruthCategoryToolCfg(flags))) + the_alg = CompFactory.DerivationFramework.TruthCategoriesDecorator(name, **kwargs) + result.addEventAlgo(the_alg, primary = True) + return result + +if __name__ == "__main__": + from MuonConfig.MuonConfigUtils import SetupMuonStandaloneArguments + from AthenaConfiguration.AllConfigFlags import ConfigFlags + args = SetupMuonStandaloneArguments() + ConfigFlags.Input.Files = args.input + ConfigFlags.Concurrency.NumThreads=args.threads + ConfigFlags.Concurrency.NumConcurrentEvents=args.threads + from AthenaCommon.Constants import DEBUG + + ConfigFlags.lock() + ConfigFlags.dump() + + from AthenaConfiguration.MainServicesConfig import MainServicesCfg + cfg = MainServicesCfg(ConfigFlags) + msgService = cfg.getService('MessageSvc') + msgService.Format = "S:%s E:%e % F%128W%S%7W%R%T %0W%M" + + from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg + cfg.merge(PoolReadCfg(ConfigFlags)) + cfg.merge(TruthCategoriesDecoratorCfg(ConfigFlags, OutputLevel = DEBUG)) + + sc = cfg.run(ConfigFlags.Exec.MaxEvents) + if not sc.isSuccess(): + exit(1) diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx index 8d5a042fa49f8dac332acb893a550c91b18b79b9..0913ff31f655201d200ebbebf4897382f6aa1e2b 100644 --- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx +++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx @@ -1,242 +1,304 @@ /* - Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration */ - #include "DerivationFrameworkHiggs/TruthCategoriesDecorator.h" -#include "xAODEventInfo/EventInfo.h" -#include "xAODTruth/TruthParticleContainer.h" -#include "xAODTruth/TruthVertex.h" -#include "xAODJet/JetContainer.h" -#include "TruthUtils/PIDHelpers.h" -#include "PathResolver/PathResolver.h" - -#include "GenInterfaces/IHiggsTruthCategoryTool.h" -#include "GenInterfaces/IxAODtoHepMCTool.h" - -#include "CLHEP/Units/SystemOfUnits.h" -// Note: must include TLorentzVector before the next one -#include "TLorentzVector.h" -#include "TruthRivetTools/HiggsTemplateCrossSectionsDefs.h" - -#include -#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include namespace DerivationFramework { - - TruthCategoriesDecorator::TruthCategoriesDecorator(const std::string& t, const std::string& n, const IInterface* p): - AthAlgTool(t,n,p), - m_xAODtoHepMCTool("xAODtoHepMCTool"), - m_higgsTruthCatTool("HiggsTruthCategoryTool"), - m_config(nullptr), - m_configPath("") - { - declareInterface(this); - declareProperty("ConfigPath",m_configPath="DerivationFrameworkHiggs/HiggsMCsamples.cfg"); - declareProperty("DetailLevel",m_detailLevel=3); - } - - TruthCategoriesDecorator::~TruthCategoriesDecorator() {} - - - StatusCode TruthCategoriesDecorator::initialize() { - - ATH_MSG_INFO("Initialize " ); - - // FOR xAOD->HEPMC :: xAODtoHepMC tool - if (m_xAODtoHepMCTool.retrieve().isFailure()) { - ATH_MSG_FATAL("Failed to retrieve tool: " << m_xAODtoHepMCTool); - return StatusCode::FAILURE; - } - ATH_MSG_INFO("Retrieved tool: " << m_xAODtoHepMCTool); - - - // Higgs truth category tool - if (m_higgsTruthCatTool.retrieve().isFailure()) { - ATH_MSG_FATAL("Failed to retrieve tool: " << m_higgsTruthCatTool); - return StatusCode::FAILURE; - } - ATH_MSG_INFO("Retrieved tool: " << m_higgsTruthCatTool); - - // Open the TEnv configuration file - m_config = new TEnv(); - int status = m_config->ReadFile(PathResolverFindCalibFile(m_configPath).c_str(),EEnvLevel(0)); - if ( status != 0 ) { - ATH_MSG_FATAL("Failed to open TEnv file "< TruthCategoriesDecorator::vectorize(const TString& str, const TString& sep) { - std::vector result; - TObjArray *strings = str.Tokenize(sep.Data()); - if (strings->GetEntries()==0) { delete strings; return result; } - TIter istr(strings); - while (TObjString* os=(TObjString*)istr()) { - result.push_back(atol(os->GetString())); - } - delete strings; - return result; - } - + TruthCategoriesDecorator::TruthCategoriesDecorator(const std::string& n, ISvcLocator* p) : AthReentrantAlgorithm(n, p) {} + + StatusCode TruthCategoriesDecorator::initialize() { + ATH_MSG_DEBUG("Initialize "); + + // FOR xAOD->HEPMC :: xAODtoHepMC tool + ATH_CHECK(m_xAODtoHepMCTool.retrieve()); + ATH_CHECK(m_higgsTruthCatTool.retrieve()); + + ATH_CHECK(m_truthEvtKey.initialize()); + ATH_CHECK(m_evtInfoKey.initialize()); + + std::set decor_keys{ + "HTXS_prodMode", + "HTXS_errorCode", + "HTXS_Stage0_Category", + // Stage 1 binning + "HTXS_Stage1_Category_pTjet25", + "HTXS_Stage1_Category_pTjet30", + "HTXS_Stage1_FineIndex_pTjet30", + "HTXS_Stage1_FineIndex_pTjet25", + // Stage 1.2 binning + "HTXS_Stage1_2_Category_pTjet25", + "HTXS_Stage1_2_Category_pTjet30", + "HTXS_Stage1_2_FineIndex_pTjet30", + "HTXS_Stage1_2_FineIndex_pTjet25", + // stage 1.2. finer binning + "HTXS_Stage1_2_Fine_Category_pTjet25", + "HTXS_Stage1_2_Fine_Category_pTjet30", + "HTXS_Stage1_2_Fine_FineIndex_pTjet30", + "HTXS_Stage1_2_Fine_FineIndex_pTjet25", + + "HTXS_Njets_pTjet25", + "HTXS_Njets_pTjet30", + "HTXS_isZ2vvDecay", + }; + + if (!m_detailLevel) decor_keys.insert("HTXS_Higgs_pt"); + // The Higgs and the associated V (last instances prior to decay) + if (m_detailLevel > 0) { + m_p4_decors.insert(std::make_pair("HTXS_Higgs", FourMomDecoration{m_evtInfoKey, "HTXS_Higgs"})); + m_p4_decors.insert(std::make_pair("HTXS_V", FourMomDecoration{m_evtInfoKey, "HTXS_V"})); + } + // Jets built excluding Higgs decay products + if (m_detailLevel > 1) { + m_p4_decors.insert(std::make_pair("HTXS_V_jets25", FourMomDecoration{m_evtInfoKey, "HTXS_V_jets25"})); + m_p4_decors.insert(std::make_pair("HTXS_V_jets30", FourMomDecoration{m_evtInfoKey, "HTXS_V_jets30"})); + } + // Everybody might not want this ... but good for validation + if (m_detailLevel > 2) { + m_p4_decors.insert(std::make_pair("HTXS_Higgs_decay", FourMomDecoration{m_evtInfoKey, "HTXS_Higgs_decay"})); + m_p4_decors.insert(std::make_pair("HTXS_V_decay", FourMomDecoration{m_evtInfoKey, "HTXS_V_decay"})); + } - HTXS::HiggsProdMode TruthCategoriesDecorator::getHiggsProductionMode(uint32_t mc_channel_number, HTXS::tH_type &th_type) const { - if (m_config==nullptr) { - ATH_MSG_ERROR("TEnv file pointer is NULL? Bad configuration."); - return HTXS::HiggsProdMode::UNKNOWN; - } + for (const std::string& key : decor_keys) m_STXSDecors.emplace_back(m_evtInfoKey.key() + "." + key); + for (auto& p4_pair : m_p4_decors) { + std::vector keys = p4_pair.second.vect(); + for (EvtInfoDecorKey key : keys) m_STXSDecors.push_back(key); + ATH_CHECK(p4_pair.second.initialize()); + + } + /// Declare everything that's written by the algorithm + ATH_CHECK(m_STXSDecors.initialize()); + // Open the TEnv configuration file + TEnv config{}; + if (config.ReadFile(PathResolverFindCalibFile(m_configPath).c_str(), EEnvLevel(0))) { + ATH_MSG_FATAL("Failed to open TEnv file " << m_configPath); + return StatusCode::FAILURE; + } - static const std::vector prodModes = { - "GGF", "VBF", "WH", "QQ2ZH", "GG2ZH", "TTH", "BBH", "TH", "THQB", "WHT" - }; - for (const TString& prodMode : prodModes) { - - // loop over each mcID belonging to the production mode - for ( int mcID : vectorize(m_config->GetValue("HTXS.MCsamples."+prodMode,"")) ){ - if (mcID==(int)mc_channel_number) { - ATH_MSG_INFO("Higgs production for MC channel number "<auxdecor((prefix+"_pt").Data()) = p4.Pt()*CLHEP::GeV; - eventInfo->auxdecor((prefix+"_eta").Data()) = p4.Eta(); - eventInfo->auxdecor((prefix+"_phi").Data()) = p4.Phi(); - eventInfo->auxdecor((prefix+"_m").Data()) = p4.M()*CLHEP::GeV; - } - - // Save a vector of TLVs as vectors of float - void TruthCategoriesDecorator::decorateFourVecs(const xAOD::EventInfo *eventInfo, const TString& prefix, const std::vector& p4s) { - std::vector pt, eta, phi, mass; - for (const auto& p4:p4s) { pt.push_back(p4.Pt()*CLHEP::GeV); eta.push_back(p4.Eta()); phi.push_back(p4.Phi()); mass.push_back(p4.M()*CLHEP::GeV); } - eventInfo->auxdecor >((prefix+"_pt").Data()) = pt; - eventInfo->auxdecor >((prefix+"_eta").Data()) = eta; - eventInfo->auxdecor >((prefix+"_phi").Data()) = phi; - eventInfo->auxdecor >((prefix+"_m").Data()) = mass; - } - - StatusCode TruthCategoriesDecorator::addBranches() const{ - - // Retrieve the xAOD event info - const xAOD::EventInfo *eventInfo = nullptr; - ATH_CHECK( evtStore()->retrieve(eventInfo,"EventInfo") ); - - // Extract the prodocution mode the first time - static bool first = true; - static HTXS::HiggsProdMode prodMode = HTXS::HiggsProdMode::UNKNOWN; - static HTXS::tH_type th_type = HTXS::tH_type::noTH; - if (first) { - uint32_t mcChannelNumber = eventInfo->mcChannelNumber(); - if(mcChannelNumber==0) mcChannelNumber = eventInfo->runNumber(); // EVNT input - prodMode = getHiggsProductionMode(mcChannelNumber,th_type); - first = false; + for (const std::string prod_mode : {"GGF", "VBF", "WH", "QQ2ZH", "GG2ZH", "TTH", "BBH", "TH", "THQB", "WHT"}) { + HTXSSample smp{}; + if (prod_mode == "GGF") + smp.prod = HTXS::HiggsProdMode::GGF; + else if (prod_mode == "VBF") + smp.prod = HTXS::HiggsProdMode::VBF; + else if (prod_mode == "WH") + smp.prod = HTXS::HiggsProdMode::WH; + else if (prod_mode == "QQ2ZH") + smp.prod = HTXS::HiggsProdMode::QQ2ZH; + else if (prod_mode == "GG2ZH") + smp.prod = HTXS::HiggsProdMode::GG2ZH; + else if (prod_mode == "TTH") + smp.prod = HTXS::HiggsProdMode::TTH; + else if (prod_mode == "BBH") + smp.prod = HTXS::HiggsProdMode::BBH; + else if (prod_mode == "TH") + smp.prod = HTXS::HiggsProdMode::TH; + else if (prod_mode == "THQB") { + smp.th_type = HTXS::tH_type::THQB; + smp.prod = HTXS::HiggsProdMode::TH; + } else if (prod_mode == "WHT") { + smp.th_type = HTXS::tH_type::TWH; + smp.prod = HTXS::HiggsProdMode::TH; + } + std::vector dsid_str{}; + MuonCalib::MdtStringUtils::tokenize(config.GetValue(Form("HTXS.MCsamples.%s", prod_mode.c_str()), ""), dsid_str, " "); + for (const std::string& dsid : dsid_str) { smp.dsids.insert(std::atoi(dsid.c_str())); } + m_htxs_samples.push_back(std::move(smp)); + } + return StatusCode::SUCCESS; } - // If the production mode is unkown, the categorization will return -99 - // Set HTXS_prodMode decoration to indicate that HTXS categorization was indeed run - if ( prodMode == HTXS::HiggsProdMode::UNKNOWN) { - eventInfo->auxdecor("HTXS_prodMode") = (int)prodMode; - return StatusCode::SUCCESS; + // Save a TLV as 4 floats + StatusCode TruthCategoriesDecorator::decorateFourVec(const EventContext& ctx, const std::string& prefix, + const TLorentzVector& p4) const { + P4DecorMap::const_iterator itr = m_p4_decors.find(prefix); + if (itr == m_p4_decors.end()) { + ATH_MSG_FATAL("decorateFourVec() -- Key " << prefix << " is unknown"); + return StatusCode::FAILURE; + } + using FloatDecorator = SG::WriteDecorHandle; + FloatDecorator dec_pt{itr->second.pt, ctx}; + FloatDecorator dec_eta{itr->second.eta, ctx}; + FloatDecorator dec_phi{itr->second.phi, ctx}; + FloatDecorator dec_m{itr->second.m, ctx}; + const xAOD::EventInfo* eventInfo = dec_pt.cptr(); + dec_pt(*eventInfo) = p4.Pt() * Gaudi::Units::GeV; + dec_eta(*eventInfo) = p4.Eta(); + dec_phi(*eventInfo) = p4.Phi(); + dec_m(*eventInfo) = p4.M() * Gaudi::Units::GeV; + ATH_MSG_DEBUG("Decorate "<retrieve(xTruthEventContainer,"TruthEvents") ); - - // convert xAOD -> HepMC - std::vector hepmc_evts = m_xAODtoHepMCTool->getHepMCEvents( xTruthEventContainer, eventInfo ); - - if (hepmc_evts.empty()) { - // ANGRY MESSAGE HERE - return StatusCode::FAILURE; + // Save a vector of TLVs as vectors of float + StatusCode TruthCategoriesDecorator::decorateFourVecs(const EventContext& ctx, const std::string& prefix, + const std::vector& p4s) const { + P4DecorMap::const_iterator itr = m_p4_decors.find(prefix); + if (itr == m_p4_decors.end()) { + ATH_MSG_FATAL("decorateFourVecs() -- Key " << prefix << " is unknown"); + return StatusCode::FAILURE; + } + using FloatDecorator = SG::WriteDecorHandle>; + FloatDecorator dec_pt{itr->second.pt, ctx}; + FloatDecorator dec_eta{itr->second.eta, ctx}; + FloatDecorator dec_phi{itr->second.phi, ctx}; + FloatDecorator dec_m{itr->second.m, ctx}; + const xAOD::EventInfo* eventInfo = dec_pt.cptr(); + + std::vector&pt{dec_pt(*eventInfo)}, eta{dec_eta(*eventInfo)}, phi{dec_phi(*eventInfo)}, mass{dec_m(*eventInfo)}; + for (const TLorentzVector& p4 : p4s) { + pt.push_back(p4.Pt() * Gaudi::Units::GeV); + eta.push_back(p4.Eta()); + phi.push_back(p4.Phi()); + mass.push_back(p4.M() * Gaudi::Units::GeV); + ATH_MSG_DEBUG("Decorate "<getHiggsTruthCategoryObject(hepmc_evts[0],prodMode); - - // Decorate the enums - eventInfo->auxdecor("HTXS_prodMode") = (int)htxs->prodMode; - eventInfo->auxdecor("HTXS_errorCode") = (int)htxs->errorCode; - eventInfo->auxdecor("HTXS_Stage0_Category") = (int)htxs->stage0_cat; - - // Stage-1 binning - eventInfo->auxdecor("HTXS_Stage1_Category_pTjet25") = (int)htxs->stage1_cat_pTjet25GeV; - eventInfo->auxdecor("HTXS_Stage1_Category_pTjet30") = (int)htxs->stage1_cat_pTjet30GeV; - - eventInfo->auxdecor("HTXS_Stage1_FineIndex_pTjet30") = HTXSstage1_to_HTXSstage1FineIndex(*htxs,th_type); - eventInfo->auxdecor("HTXS_Stage1_FineIndex_pTjet25") = HTXSstage1_to_HTXSstage1FineIndex(*htxs,th_type,true); - - // Stage-1.2 binning - eventInfo->auxdecor("HTXS_Stage1_2_Category_pTjet25") = (int)htxs->stage1_2_cat_pTjet25GeV; - eventInfo->auxdecor("HTXS_Stage1_2_Category_pTjet30") = (int)htxs->stage1_2_cat_pTjet30GeV; - - eventInfo->auxdecor("HTXS_Stage1_2_FineIndex_pTjet30") = HTXSstage1_2_to_HTXSstage1_2_FineIndex(*htxs,th_type); - eventInfo->auxdecor("HTXS_Stage1_2_FineIndex_pTjet25") = HTXSstage1_2_to_HTXSstage1_2_FineIndex(*htxs,th_type,true); - - // Stage-1.2 finer binning - eventInfo->auxdecor("HTXS_Stage1_2_Fine_Category_pTjet25") = (int)htxs->stage1_2_fine_cat_pTjet25GeV; - eventInfo->auxdecor("HTXS_Stage1_2_Fine_Category_pTjet30") = (int)htxs->stage1_2_fine_cat_pTjet30GeV; + StatusCode TruthCategoriesDecorator::execute(const EventContext& ctx) const { + using IntDecorator = SG::AuxElement::Decorator; + using FloatDecorator = SG::AuxElement::Decorator; + // Retrieve the xAOD event info + SG::ReadHandle eventInfo{m_evtInfoKey, ctx}; + if (!eventInfo.isValid()) { + ATH_MSG_FATAL("Failed to retrieve " << m_evtInfoKey.fullKey()); + return StatusCode::FAILURE; + } - eventInfo->auxdecor("HTXS_Stage1_2_Fine_FineIndex_pTjet30") = HTXSstage1_2_Fine_to_HTXSstage1_2_Fine_FineIndex(*htxs,th_type); - eventInfo->auxdecor("HTXS_Stage1_2_Fine_FineIndex_pTjet25") = HTXSstage1_2_Fine_to_HTXSstage1_2_Fine_FineIndex(*htxs,th_type,true); + int mcChannelNumber = eventInfo->mcChannelNumber(); + if (!mcChannelNumber) mcChannelNumber = eventInfo->runNumber(); // EVNT input + + static const IntDecorator dec_prodMode{"HTXS_prodMode"}; + std::vector::const_iterator smp_itr = + std::find_if(m_htxs_samples.begin(), m_htxs_samples.end(), + [mcChannelNumber](const HTXSSample& smp) { return smp.dsids.count(mcChannelNumber); }); + if (smp_itr == m_htxs_samples.end()) { + dec_prodMode(*eventInfo) = HTXS::HiggsProdMode::UNKNOWN; + ATH_MSG_DEBUG("The sample " << mcChannelNumber + << " is not a Higgs sample and not of further interest. Decorate the prod mode information only"); + return StatusCode::SUCCESS; + } - eventInfo->auxdecor("HTXS_Njets_pTjet25") = (int)htxs->jets25.size(); - eventInfo->auxdecor("HTXS_Njets_pTjet30") = (int)htxs->jets30.size(); - - eventInfo->auxdecor("HTXS_isZ2vvDecay") = (bool)htxs->isZ2vvDecay; + const HTXS::HiggsProdMode prodMode = smp_itr->prod; + const HTXS::tH_type th_type = smp_itr->th_type; - // At the very least, save the Higgs boson pT - if (m_detailLevel==0) eventInfo->auxdecor("HTXS_Higgs_pt") = htxs->higgs.Pt()*CLHEP::GeV; + // Retrieve the xAOD truth + SG::ReadHandle xTruthEventContainer{m_truthEvtKey, ctx}; + if (!xTruthEventContainer.isValid()) { + ATH_MSG_FATAL("Failed to retrieve " << m_truthEvtKey.fullKey()); + return StatusCode::FAILURE; + } + // convert xAOD -> HepMC + std::vector hepmc_evts = m_xAODtoHepMCTool->getHepMCEvents(xTruthEventContainer.cptr(), eventInfo.cptr()); + + if (hepmc_evts.empty()) { + // ANGRY MESSAGE HERE + ATH_MSG_FATAL("The HEP MC GenEvent conversion failed"); + return StatusCode::FAILURE; + } else { + ATH_MSG_DEBUG("Found "<0) { - // The Higgs and the associated V (last instances prior to decay) - decorateFourVec(eventInfo,"HTXS_Higgs",htxs->higgs); - decorateFourVec(eventInfo,"HTXS_V",htxs->V); - } + // classify event according to simplified template cross section + std::unique_ptr htxs{m_higgsTruthCatTool->getHiggsTruthCategoryObject(hepmc_evts[0], prodMode)}; + ATH_MSG_DEBUG("Truth categorization done "); + // Decorate the enums + static const IntDecorator dec_errorCode{"HTXS_errorCode"}; + static const IntDecorator dec_stage0Cat{"HTXS_Stage0_Category"}; + + dec_prodMode(*eventInfo) = htxs->prodMode; + dec_errorCode(*eventInfo) = htxs->errorCode; + dec_stage0Cat(*eventInfo) = htxs->stage0_cat; + + // Stage-1 binning + static const IntDecorator dec_stage1CatPt25{"HTXS_Stage1_Category_pTjet25"}; + static const IntDecorator dec_stage1CatPt30{"HTXS_Stage1_Category_pTjet30"}; + static const IntDecorator dec_stage1IdxPt25{"HTXS_Stage1_FineIndex_pTjet25"}; + static const IntDecorator dec_stage1IdxPt30{"HTXS_Stage1_FineIndex_pTjet30"}; + + dec_stage1CatPt25(*eventInfo) = htxs->stage1_cat_pTjet25GeV; + dec_stage1CatPt30(*eventInfo) = htxs->stage1_cat_pTjet30GeV; + /// Last argument switches between 25 GeV (true) and 30 GeV jets (false) + dec_stage1IdxPt25(*eventInfo) = HTXSstage1_to_HTXSstage1FineIndex(*htxs, th_type, true); + dec_stage1IdxPt30(*eventInfo) = HTXSstage1_to_HTXSstage1FineIndex(*htxs, th_type, false); + + // Stage-1.2 binning + static const IntDecorator dec_stage1p2_CatPt25{"HTXS_Stage1_2_Category_pTjet25"}; + static const IntDecorator dec_stage1p2_CatPt30{"HTXS_Stage1_2_Category_pTjet30"}; + static const IntDecorator dec_stage1p2_IdxPt25{"HTXS_Stage1_2_FineIndex_pTjet25"}; + static const IntDecorator dec_stage1p2_IdxPt30{"HTXS_Stage1_2_FineIndex_pTjet30"}; + + dec_stage1p2_CatPt25(*eventInfo) = htxs->stage1_2_cat_pTjet25GeV; + dec_stage1p2_CatPt30(*eventInfo) = htxs->stage1_2_cat_pTjet30GeV; + /// Last argument switches between 25 (true) / 30 (false) GeV jets. + dec_stage1p2_IdxPt25(*eventInfo) = HTXSstage1_2_to_HTXSstage1_2_FineIndex(*htxs, th_type, true); + dec_stage1p2_IdxPt30(*eventInfo) = HTXSstage1_2_to_HTXSstage1_2_FineIndex(*htxs, th_type, false); + + // Stage-1.2 finer binning + static const IntDecorator dec_stage1p2_Fine_CatPt25{"HTXS_Stage1_2_Fine_Category_pTjet25"}; + static const IntDecorator dec_stage1p2_Fine_CatPt30{"HTXS_Stage1_2_Fine_Category_pTjet30"}; + static const IntDecorator dec_stage1p2_Fine_IdxPt25{"HTXS_Stage1_2_Fine_FineIndex_pTjet25"}; + static const IntDecorator dec_stage1p2_Fine_IdxPt30{"HTXS_Stage1_2_Fine_FineIndex_pTjet30"}; + + dec_stage1p2_Fine_CatPt25(*eventInfo) = htxs->stage1_2_fine_cat_pTjet25GeV; + dec_stage1p2_Fine_CatPt30(*eventInfo) = htxs->stage1_2_fine_cat_pTjet30GeV; + dec_stage1p2_Fine_IdxPt25(*eventInfo) = HTXSstage1_2_Fine_to_HTXSstage1_2_Fine_FineIndex(*htxs, th_type, true); + dec_stage1p2_Fine_IdxPt30(*eventInfo) = HTXSstage1_2_Fine_to_HTXSstage1_2_Fine_FineIndex(*htxs, th_type, false); + + static const IntDecorator dec_NJets25{"HTXS_Njets_pTjet25"}; + static const IntDecorator dec_NJets30{"HTXS_Njets_pTjet30"}; + dec_NJets25(*eventInfo) = htxs->jets25.size(); + dec_NJets30(*eventInfo) = htxs->jets30.size(); + + static const IntDecorator dec_isZnunu{"HTXS_isZ2vvDecay"}; + dec_isZnunu(*eventInfo) = htxs->isZ2vvDecay; + + // At the very least, save the Higgs boson pT + if (!m_detailLevel) { + static const FloatDecorator dec_higgsPt{"HTXS_Higgs_pt"}; + dec_higgsPt(*eventInfo) = htxs->higgs.Pt() * Gaudi::Units::GeV; + } + if (m_detailLevel > 0) { + // The Higgs and the associated V (last instances prior to decay) + ATH_CHECK(decorateFourVec(ctx, "HTXS_Higgs", htxs->higgs)); + ATH_CHECK(decorateFourVec(ctx, "HTXS_V", htxs->V)); + } - if (m_detailLevel>1) { - // Jets built excluding Higgs decay products - decorateFourVecs(eventInfo,"HTXS_V_jets25",htxs->jets25); - decorateFourVecs(eventInfo,"HTXS_V_jets30",htxs->jets30); - } + if (m_detailLevel > 1) { + // Jets built excluding Higgs decay products + ATH_CHECK(decorateFourVecs(ctx, "HTXS_V_jets25", htxs->jets25)); + ATH_CHECK(decorateFourVecs(ctx, "HTXS_V_jets30", htxs->jets30)); + } - if (m_detailLevel>2) { - // Everybody might not want this ... but good for validation - decorateFourVec(eventInfo,"HTXS_Higgs_decay",htxs->p4decay_higgs); - decorateFourVec(eventInfo,"HTXS_V_decay",htxs->p4decay_V); + if (m_detailLevel > 2) { + // Everybody might not want this ... but good for validation + ATH_CHECK(decorateFourVec(ctx, "HTXS_Higgs_decay", htxs->p4decay_higgs)); + ATH_CHECK(decorateFourVec(ctx, "HTXS_V_decay", htxs->p4decay_V)); + } + /// Summary of the HTXS categorization. Used mainly in the testing algorithm + ATH_MSG_DEBUG("production mode: " << dec_prodMode(*eventInfo) << ", errorCode: " << dec_errorCode(*eventInfo) << ", Stage0: " << dec_stage0Cat(*eventInfo) + << ", " << dec_stage1CatPt25(*eventInfo) << ", Stage1 -- Jet25 " << dec_stage1CatPt25(*eventInfo) << " Idx " << dec_stage1IdxPt25(*eventInfo) + << ", Jet30: " << dec_stage1CatPt30(*eventInfo) << " Idx "<< dec_stage1IdxPt30(*eventInfo) << ", Stage 1.2 -- Jet25: " + << dec_stage1p2_CatPt25(*eventInfo) << " Idx: " << dec_stage1p2_IdxPt25(*eventInfo) << " Jet30: " << dec_stage1p2_CatPt30(*eventInfo) + << " Idx: " << dec_stage1p2_IdxPt30(*eventInfo) << ", HTXS NJets 25" << dec_NJets25(*eventInfo) + << " HTXS NJets 30 " << dec_NJets30(*eventInfo) << " Z->nunu" << dec_isZnunu(*eventInfo)); + + return StatusCode::SUCCESS; } - delete htxs; - - return StatusCode::SUCCESS; - - } // addBranches - -} // namespace +} // namespace DerivationFramework