From 972d51a21df035c7971451f8021161c9e15e3cc1 Mon Sep 17 00:00:00 2001 From: Johannes Josef Junggeburth Date: Mon, 29 Aug 2022 08:53:05 +0200 Subject: [PATCH 1/8] Revive the TruthCategories decorator as AthReentrant algorithm --- .../TruthCategoriesDecorator.h | 86 ++-- .../src/TruthCategoriesDecorator.cxx | 382 ++++++++++-------- 2 files changed, 256 insertions(+), 212 deletions(-) diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/DerivationFrameworkHiggs/TruthCategoriesDecorator.h b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/DerivationFrameworkHiggs/TruthCategoriesDecorator.h index 6fabb2dfe8a..3ab9324974f 100644 --- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/DerivationFrameworkHiggs/TruthCategoriesDecorator.h +++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/DerivationFrameworkHiggs/TruthCategoriesDecorator.h @@ -1,5 +1,5 @@ /* - 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 @@ -8,67 +8,81 @@ #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 "xAODTruth/TruthEventContainer.h" #include "xAODEventInfo/EventInfo.h" + // Note: must include TLorentzVector before the next one -#include "TLorentzVector.h" -#include "TruthRivetTools/HiggsTemplateCrossSectionsDefs.h" -class IHiggsTruthCategoryTool; -class IxAODtoHepMCTool; +#include "TruthRivetTools/HiggsTemplateCrossSectionsDefs.h" +#include "GenInterfaces/IHiggsTruthCategoryTool.h" +#include "GenInterfaces/IxAODtoHepMCTool.h" namespace DerivationFramework { - class TruthCategoriesDecorator : public AthAlgTool, public IAugmentationTool { + class TruthCategoriesDecorator : public AthReentrantAlgorithm { public: - TruthCategoriesDecorator(const std::string& t, const std::string& n, const IInterface* p); - ~TruthCategoriesDecorator(); + TruthCategoriesDecorator(const std::string& n, ISvcLocator* p); + virtual ~TruthCategoriesDecorator() =default; StatusCode initialize(); - StatusCode finalize(); - virtual StatusCode addBranches() const; + StatusCode execute(const EventContext& ctx) const; private: - ToolHandle m_xAODtoHepMCTool; - ToolHandle m_higgsTruthCatTool; + ToolHandle m_xAODtoHepMCTool{this, "HepMCTool", "xAODtoHepMCTool"}; + ToolHandle m_higgsTruthCatTool{this, "CategoryTool", "HiggsTruthCategoryTool" }; - // Path to TEnv file containing MC-channel-numbre <-> HiggsProdMode map - TEnv *m_config; - std::string m_configPath; + // 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 - 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=" ") ; + Gaudi::Property m_detailLevel{this, "DetailLevel", 3}; - // 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 ) ; + 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};} + }; + using P4DecorMap = std::map; + P4DecorMap m_p4_decors{}; }; /// class } /// namespace diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx index 8d5a042fa49..09cb4a7de52 100644 --- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx +++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx @@ -4,239 +4,269 @@ #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() {} - + TruthCategoriesDecorator::TruthCategoriesDecorator(const std::string& n, ISvcLocator* p) : AthReentrantAlgorithm(n, p) {} + StatusCode TruthCategoriesDecorator::initialize() { - ATH_MSG_INFO("Initialize " ); + ATH_MSG_DEBUG("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_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"} )); } - ATH_MSG_INFO("Retrieved tool: " << m_higgsTruthCatTool); + for (const std::string& key : decor_keys) m_STXSDecors.emplace_back(m_evtInfoKey.key()+"."+key); + for (const auto& p4_pair : m_p4_decors){ + std::vector keys = p4_pair.second.vect(); + for (EvtInfoDecorKey& key: keys) m_STXSDecors.push_back(std::move(key)); + } + /// Declare everything that's written by the algorithm + ATH_CHECK(m_STXSDecors.initialize()); // Open the TEnv configuration file - m_config = new TEnv(); - int status = m_config->ReadFile(PathResolverFindCalibFile(m_configPath).c_str(),EEnvLevel(0)); - if ( status != 0 ) { + TEnv config{}; + if (config.ReadFile(PathResolverFindCalibFile(m_configPath).c_str(),EEnvLevel(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; - } - - - - 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; - } - - 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 "< 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)); } - // This is perfectly fine if we aren't running on a Higgs sample - ATH_MSG_INFO("Did not manage to extract Higgs production mode for MC channel number " << - mc_channel_number << ". HTXS categorization will hence not be derived."); - return HTXS::HiggsProdMode::UNKNOWN; + return StatusCode::SUCCESS; } // Save a TLV as 4 floats - void TruthCategoriesDecorator::decorateFourVec(const xAOD::EventInfo *eventInfo, const TString& prefix, const TLorentzVector& p4) { - eventInfo->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; + 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 "<; + 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; + return StatusCode::SUCCESS; } // 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::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 "<>; + 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); + } + return StatusCode::SUCCESS; } - StatusCode TruthCategoriesDecorator::addBranches() const{ - + StatusCode TruthCategoriesDecorator::execute(const EventContext& ctx ) const{ + using IntDecorator = SG::AuxElement::Decorator; + using FloatDecorator = SG::AuxElement::Decorator; // 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; + SG::ReadHandle eventInfo{m_evtInfoKey, ctx}; + if (!eventInfo.isValid()) { + ATH_MSG_FATAL("Failed to retrieve "<auxdecor("HTXS_prodMode") = (int)prodMode; + + + 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 "<prod; + const HTXS::tH_type th_type = smp_itr->th_type; + - // Retrieve the xAOD truth - const xAOD::TruthEventContainer* xTruthEventContainer = nullptr; - ATH_CHECK( evtStore()->retrieve(xTruthEventContainer,"TruthEvents") ); - + // Retrieve the xAOD truth + SG::ReadHandle xTruthEventContainer{m_truthEvtKey, ctx}; + if (!xTruthEventContainer.isValid()) { + ATH_MSG_FATAL( "Failed to retrieve "< HepMC - std::vector hepmc_evts = m_xAODtoHepMCTool->getHepMCEvents( xTruthEventContainer, eventInfo ); + 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; } - + // classify event according to simplified template cross section - HTXS::HiggsClassification *htxs = m_higgsTruthCatTool->getHiggsTruthCategoryObject(hepmc_evts[0],prodMode); + std::unique_ptr htxs{ m_higgsTruthCatTool->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; - + // 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 - 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); - + 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_pTjet30"}; + static const IntDecorator dec_stage1IdxPt30{"HTXS_Stage1_FineIndex_pTjet25"}; + + dec_stage1CatPt25(*eventInfo) = htxs->stage1_cat_pTjet25GeV; + dec_stage1CatPt30(*eventInfo) = htxs->stage1_cat_pTjet25GeV; + dec_stage1IdxPt25(*eventInfo) = HTXSstage1_to_HTXSstage1FineIndex(*htxs,th_type); + dec_stage1IdxPt30(*eventInfo) = 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); + 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_pTjet30"}; + static const IntDecorator dec_stage1p2_IdxPt30{"HTXS_Stage1_2_FineIndex_pTjet25"}; + + dec_stage1p2_CatPt25(*eventInfo) = htxs->stage1_2_cat_pTjet25GeV; + dec_stage1p2_CatPt30(*eventInfo) = htxs->stage1_2_cat_pTjet30GeV; + dec_stage1p2_IdxPt25(*eventInfo) = HTXSstage1_2_to_HTXSstage1_2_FineIndex(*htxs,th_type); + dec_stage1p2_IdxPt30(*eventInfo) = 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; - - 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); - - eventInfo->auxdecor("HTXS_Njets_pTjet25") = (int)htxs->jets25.size(); - eventInfo->auxdecor("HTXS_Njets_pTjet30") = (int)htxs->jets30.size(); + 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_pTjet30"}; + static const IntDecorator dec_stage1p2_Fine_IdxPt30{"HTXS_Stage1_2_Fine_FineIndex_pTjet25"}; + + 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); + dec_stage1p2_Fine_IdxPt30(*eventInfo) = HTXSstage1_2_Fine_to_HTXSstage1_2_Fine_FineIndex(*htxs,th_type,true); + + static const IntDecorator dec_NJets25{"HTXS_Njets_pTjet25"}; + static const IntDecorator dec_NJets30{"HTXS_Njets_pTjet30"}; + dec_NJets25(*eventInfo) = htxs->jets25.size(); + dec_NJets25(*eventInfo) = htxs->jets30.size(); - eventInfo->auxdecor("HTXS_isZ2vvDecay") = (bool)htxs->isZ2vvDecay; + static const IntDecorator dec_isZnunu{"HTXS_isZ2vvDecay"}; + dec_isZnunu(*eventInfo) = htxs->isZ2vvDecay; // At the very least, save the Higgs boson pT - if (m_detailLevel==0) eventInfo->auxdecor("HTXS_Higgs_pt") = htxs->higgs.Pt()*CLHEP::GeV; - + 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) - decorateFourVec(eventInfo,"HTXS_Higgs",htxs->higgs); - decorateFourVec(eventInfo,"HTXS_V",htxs->V); + 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); + 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); + ATH_CHECK(decorateFourVec(ctx,"HTXS_Higgs_decay",htxs->p4decay_higgs)); + ATH_CHECK(decorateFourVec(ctx,"HTXS_V_decay",htxs->p4decay_V)); } - delete htxs; - return StatusCode::SUCCESS; - - } // addBranches + } } // namespace -- GitLab From 543b8318cda03259626ae7bc15a45320a6784426 Mon Sep 17 00:00:00 2001 From: Johannes Josef Junggeburth Date: Mon, 29 Aug 2022 09:22:12 +0200 Subject: [PATCH 2/8] Clang-format --- .../TruthCategoriesDecorator.h | 130 +++-- .../src/TruthCategoriesDecorator.cxx | 513 +++++++++--------- 2 files changed, 329 insertions(+), 314 deletions(-) diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/DerivationFrameworkHiggs/TruthCategoriesDecorator.h b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/DerivationFrameworkHiggs/TruthCategoriesDecorator.h index 3ab9324974f..4965e9756b8 100644 --- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/DerivationFrameworkHiggs/TruthCategoriesDecorator.h +++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/DerivationFrameworkHiggs/TruthCategoriesDecorator.h @@ -2,9 +2,9 @@ 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 @@ -14,77 +14,75 @@ #include "StoreGate/WriteDecorHandleKeyArray.h" // EDM: typedefs, so must be included and not forward-referenced -#include "xAODTruth/TruthEventContainer.h" #include "xAODEventInfo/EventInfo.h" - +#include "xAODTruth/TruthEventContainer.h" // Note: must include TLorentzVector before the next one -#include "TruthRivetTools/HiggsTemplateCrossSectionsDefs.h" #include "GenInterfaces/IHiggsTruthCategoryTool.h" #include "GenInterfaces/IxAODtoHepMCTool.h" +#include "TruthRivetTools/HiggsTemplateCrossSectionsDefs.h" namespace DerivationFramework { - - 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};} - }; - using P4DecorMap = std::map; - P4DecorMap m_p4_decors{}; - }; /// 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}; } + }; + using P4DecorMap = std::map; + P4DecorMap m_p4_decors{}; + }; /// class + +} // namespace DerivationFramework #endif diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx index 09cb4a7de52..5b80ebbc60a 100644 --- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx +++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx @@ -2,271 +2,288 @@ Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration */ - #include "DerivationFrameworkHiggs/TruthCategoriesDecorator.h" -#include -#include -#include -#include -#include #include +#include +#include #include #include -#include -#include #include +#include #include +#include +#include +#include +#include namespace DerivationFramework { - 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"} )); - } - - for (const std::string& key : decor_keys) m_STXSDecors.emplace_back(m_evtInfoKey.key()+"."+key); - for (const auto& p4_pair : m_p4_decors){ - std::vector keys = p4_pair.second.vect(); - for (EvtInfoDecorKey& key: keys) m_STXSDecors.push_back(std::move(key)); - } - /// 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 "< 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())); + 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"})); } - m_htxs_samples.push_back(std::move(smp)); - } - 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 "<; - 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; - return StatusCode::SUCCESS; - } - - // 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 "<>; - 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); - } - return StatusCode::SUCCESS; - } - - 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 "<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 "<prod; - const HTXS::tH_type th_type = smp_itr->th_type; - - - // Retrieve the xAOD truth - SG::ReadHandle xTruthEventContainer{m_truthEvtKey, ctx}; - if (!xTruthEventContainer.isValid()) { - ATH_MSG_FATAL( "Failed to retrieve "< 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; - } - - // classify event according to simplified template cross section - std::unique_ptr htxs{ m_higgsTruthCatTool->getHiggsTruthCategoryObject(hepmc_evts[0],prodMode)}; - - // 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_pTjet30"}; - static const IntDecorator dec_stage1IdxPt30{"HTXS_Stage1_FineIndex_pTjet25"}; - - dec_stage1CatPt25(*eventInfo) = htxs->stage1_cat_pTjet25GeV; - dec_stage1CatPt30(*eventInfo) = htxs->stage1_cat_pTjet25GeV; - dec_stage1IdxPt25(*eventInfo) = HTXSstage1_to_HTXSstage1FineIndex(*htxs,th_type); - dec_stage1IdxPt30(*eventInfo) = HTXSstage1_to_HTXSstage1FineIndex(*htxs,th_type, true); - - - // 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_pTjet30"}; - static const IntDecorator dec_stage1p2_IdxPt30{"HTXS_Stage1_2_FineIndex_pTjet25"}; - - dec_stage1p2_CatPt25(*eventInfo) = htxs->stage1_2_cat_pTjet25GeV; - dec_stage1p2_CatPt30(*eventInfo) = htxs->stage1_2_cat_pTjet30GeV; - dec_stage1p2_IdxPt25(*eventInfo) = HTXSstage1_2_to_HTXSstage1_2_FineIndex(*htxs,th_type); - dec_stage1p2_IdxPt30(*eventInfo) = HTXSstage1_2_to_HTXSstage1_2_FineIndex(*htxs,th_type,true); - - // 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_pTjet30"}; - static const IntDecorator dec_stage1p2_Fine_IdxPt30{"HTXS_Stage1_2_Fine_FineIndex_pTjet25"}; - - 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); - dec_stage1p2_Fine_IdxPt30(*eventInfo) = HTXSstage1_2_Fine_to_HTXSstage1_2_Fine_FineIndex(*htxs,th_type,true); - - static const IntDecorator dec_NJets25{"HTXS_Njets_pTjet25"}; - static const IntDecorator dec_NJets30{"HTXS_Njets_pTjet30"}; - dec_NJets25(*eventInfo) = htxs->jets25.size(); - dec_NJets25(*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; + for (const std::string& key : decor_keys) m_STXSDecors.emplace_back(m_evtInfoKey.key() + "." + key); + for (const auto& p4_pair : m_p4_decors) { + std::vector keys = p4_pair.second.vect(); + for (EvtInfoDecorKey& key : keys) m_STXSDecors.push_back(std::move(key)); + } + /// 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; + } + + 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 (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)); + + // 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; + return StatusCode::SUCCESS; } - 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)); + // 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); + } + return StatusCode::SUCCESS; } - 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)); + 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; + } + + 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; + } + + const HTXS::HiggsProdMode prodMode = smp_itr->prod; + const HTXS::tH_type th_type = smp_itr->th_type; + + // 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; + } + + // classify event according to simplified template cross section + std::unique_ptr htxs{m_higgsTruthCatTool->getHiggsTruthCategoryObject(hepmc_evts[0], prodMode)}; + + // 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_pTjet30"}; + static const IntDecorator dec_stage1IdxPt30{"HTXS_Stage1_FineIndex_pTjet25"}; + + dec_stage1CatPt25(*eventInfo) = htxs->stage1_cat_pTjet25GeV; + dec_stage1CatPt30(*eventInfo) = htxs->stage1_cat_pTjet25GeV; + dec_stage1IdxPt25(*eventInfo) = HTXSstage1_to_HTXSstage1FineIndex(*htxs, th_type); + dec_stage1IdxPt30(*eventInfo) = HTXSstage1_to_HTXSstage1FineIndex(*htxs, th_type, true); + + // 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_pTjet30"}; + static const IntDecorator dec_stage1p2_IdxPt30{"HTXS_Stage1_2_FineIndex_pTjet25"}; + + dec_stage1p2_CatPt25(*eventInfo) = htxs->stage1_2_cat_pTjet25GeV; + dec_stage1p2_CatPt30(*eventInfo) = htxs->stage1_2_cat_pTjet30GeV; + dec_stage1p2_IdxPt25(*eventInfo) = HTXSstage1_2_to_HTXSstage1_2_FineIndex(*htxs, th_type); + dec_stage1p2_IdxPt30(*eventInfo) = HTXSstage1_2_to_HTXSstage1_2_FineIndex(*htxs, th_type, true); + + // 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_pTjet30"}; + static const IntDecorator dec_stage1p2_Fine_IdxPt30{"HTXS_Stage1_2_Fine_FineIndex_pTjet25"}; + + 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); + dec_stage1p2_Fine_IdxPt30(*eventInfo) = HTXSstage1_2_Fine_to_HTXSstage1_2_Fine_FineIndex(*htxs, th_type, true); + + static const IntDecorator dec_NJets25{"HTXS_Njets_pTjet25"}; + static const IntDecorator dec_NJets30{"HTXS_Njets_pTjet30"}; + dec_NJets25(*eventInfo) = htxs->jets25.size(); + dec_NJets25(*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 + 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 + ATH_CHECK(decorateFourVec(ctx, "HTXS_Higgs_decay", htxs->p4decay_higgs)); + ATH_CHECK(decorateFourVec(ctx, "HTXS_V_decay", htxs->p4decay_V)); + } + + return StatusCode::SUCCESS; } - return StatusCode::SUCCESS; - } - -} // namespace +} // namespace DerivationFramework -- GitLab From dade532aab9bbae62cdcf026931d6c2d19e4b224 Mon Sep 17 00:00:00 2001 From: Johannes Josef Junggeburth Date: Mon, 29 Aug 2022 09:55:41 +0200 Subject: [PATCH 3/8] Write configuration --- .../GenInterfaces/IHiggsTruthCategoryTool.h | 11 +++--- Generators/TruthConverters/CMakeLists.txt | 2 + .../python/TruthConvertersCfg.py | 8 ++++ Generators/TruthRivetTools/CMakeLists.txt | 2 + .../python/TruthRivetToolsCfg.py | 8 ++++ .../DerivationFrameworkHiggs/CMakeLists.txt | 2 +- .../python/TruthCategories.py | 13 ++----- .../python/TruthCategoriesConfig.py | 38 +++++++++++++++++++ .../src/TruthCategoriesDecorator.cxx | 2 +- 9 files changed, 68 insertions(+), 18 deletions(-) create mode 100644 Generators/TruthConverters/python/TruthConvertersCfg.py create mode 100644 Generators/TruthRivetTools/python/TruthRivetToolsCfg.py create mode 100644 PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/TruthCategoriesConfig.py diff --git a/Generators/GenInterfaces/GenInterfaces/IHiggsTruthCategoryTool.h b/Generators/GenInterfaces/GenInterfaces/IHiggsTruthCategoryTool.h index dd8b1e0de73..52dce88eea0 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 90155d36f1e..14e0bedec85 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/python/TruthConvertersCfg.py b/Generators/TruthConverters/python/TruthConvertersCfg.py new file mode 100644 index 00000000000..fe209e8c5f2 --- /dev/null +++ b/Generators/TruthConverters/python/TruthConvertersCfg.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 xAODtoHEPToolCfg(flags,name="xAODToHepMCTool", **kwargs): + result = ComponentAccumulator() + result.setPrivateTools(CompFactory.xAODtoHepMCTool(name,**kwargs)) + return result \ No newline at end of file diff --git a/Generators/TruthRivetTools/CMakeLists.txt b/Generators/TruthRivetTools/CMakeLists.txt index bce762fd485..0c7b36f7ec4 100644 --- a/Generators/TruthRivetTools/CMakeLists.txt +++ b/Generators/TruthRivetTools/CMakeLists.txt @@ -33,3 +33,5 @@ atlas_add_component( TruthRivetTools atlas_add_dictionary( TruthRivetToolsDict TruthRivetTools/TruthRivetToolsDict.h TruthRivetTools/selection.xml LINK_LIBRARIES TruthRivetToolsLib ) + +atlas_install_python_modules( python/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} ) diff --git a/Generators/TruthRivetTools/python/TruthRivetToolsCfg.py b/Generators/TruthRivetTools/python/TruthRivetToolsCfg.py new file mode 100644 index 00000000000..8bb2ad6a0a6 --- /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 42bbbe30c58..e7affbd8f35 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/python/TruthCategories.py b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/TruthCategories.py index 4d01b47e11f..1e449333c35 100644 --- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/TruthCategories.py +++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/TruthCategories.py @@ -1,17 +1,10 @@ -# 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 +if globalflags.DataSource()=='geant4': import AthenaCommon.CfgMgr as CfgMgr - 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 00000000000..fc4b3d7fa8b --- /dev/null +++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/TruthCategoriesConfig.py @@ -0,0 +1,38 @@ +# 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 + + 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)) + + 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 5b80ebbc60a..8740609400d 100644 --- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx +++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx @@ -1,5 +1,5 @@ /* - 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" -- GitLab From e719695c825b4c8ca49314acafb5add61ed69d34 Mon Sep 17 00:00:00 2001 From: Johannes Josef Junggeburth Date: Mon, 29 Aug 2022 13:21:52 +0200 Subject: [PATCH 4/8] Add debug lines for testing macro --- .../python/TruthCategoriesConfig.py | 2 ++ .../src/TruthCategoriesDecorator.cxx | 11 ++++++++++- 2 files changed, 12 insertions(+), 1 deletion(-) diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/TruthCategoriesConfig.py b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/TruthCategoriesConfig.py index fc4b3d7fa8b..f9bb2e9c0a5 100644 --- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/TruthCategoriesConfig.py +++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/TruthCategoriesConfig.py @@ -20,6 +20,8 @@ if __name__ == "__main__": ConfigFlags.Input.Files = args.input ConfigFlags.Concurrency.NumThreads=args.threads ConfigFlags.Concurrency.NumConcurrentEvents=args.threads + from AthenaCommon.Constants import DEBUG + ConfigFlags.Exec.OutputLevel = DEBUG #Global Output Level ConfigFlags.lock() ConfigFlags.dump() diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx index 8740609400d..8f0a8ea9a0e 100644 --- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx +++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx @@ -138,6 +138,7 @@ namespace DerivationFramework { dec_eta(*eventInfo) = p4.Eta(); dec_phi(*eventInfo) = p4.Phi(); dec_m(*eventInfo) = p4.M() * Gaudi::Units::GeV; + ATH_MSG_DEBUG("Decorate "<jets25.size(); - dec_NJets25(*eventInfo) = htxs->jets30.size(); + dec_NJets30(*eventInfo) = htxs->jets30.size(); static const IntDecorator dec_isZnunu{"HTXS_isZ2vvDecay"}; dec_isZnunu(*eventInfo) = htxs->isZ2vvDecay; @@ -282,6 +284,13 @@ namespace DerivationFramework { 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; } -- GitLab From 5ba9d0769a853f8f2af356c9bdb828b0d273f6a7 Mon Sep 17 00:00:00 2001 From: Johannes Josef Junggeburth Date: Mon, 29 Aug 2022 13:27:04 +0200 Subject: [PATCH 5/8] Whitespace fixes and it's not very wise to set the avalanche scheduler to DEBUG --- .../python/TruthCategoriesConfig.py | 5 ++--- .../src/TruthCategoriesDecorator.cxx | 8 ++++---- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/TruthCategoriesConfig.py b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/TruthCategoriesConfig.py index f9bb2e9c0a5..d8043d52ddf 100644 --- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/TruthCategoriesConfig.py +++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/TruthCategoriesConfig.py @@ -21,8 +21,7 @@ if __name__ == "__main__": ConfigFlags.Concurrency.NumThreads=args.threads ConfigFlags.Concurrency.NumConcurrentEvents=args.threads from AthenaCommon.Constants import DEBUG - ConfigFlags.Exec.OutputLevel = DEBUG #Global Output Level - + ConfigFlags.lock() ConfigFlags.dump() @@ -33,7 +32,7 @@ if __name__ == "__main__": from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg cfg.merge(PoolReadCfg(ConfigFlags)) - cfg.merge(TruthCategoriesDecoratorCfg(ConfigFlags)) + cfg.merge(TruthCategoriesDecoratorCfg(ConfigFlags, OutputLevel = DEBUG)) sc = cfg.run(ConfigFlags.Exec.MaxEvents) if not sc.isSuccess(): diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx index 8f0a8ea9a0e..cfc8308539a 100644 --- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx +++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx @@ -33,12 +33,12 @@ namespace DerivationFramework { std::set decor_keys{ "HTXS_prodMode", - " HTXS_errorCode", - " HTXS_Stage0_Category", + "HTXS_errorCode", + "HTXS_Stage0_Category", // Stage 1 binning "HTXS_Stage1_Category_pTjet25", - " HTXS_Stage1_Category_pTjet30", - " HTXS_Stage1_FineIndex_pTjet30", + "HTXS_Stage1_Category_pTjet30", + "HTXS_Stage1_FineIndex_pTjet30", "HTXS_Stage1_FineIndex_pTjet25", // Stage 1.2 binning "HTXS_Stage1_2_Category_pTjet25", -- GitLab From d4f413d181e79aa5062aa5f46d962cf53e3ae463 Mon Sep 17 00:00:00 2001 From: Johannes Josef Junggeburth Date: Mon, 29 Aug 2022 14:43:14 +0200 Subject: [PATCH 6/8] Few runtime fixes.. Double deletion... uninitialized keys... etc --- .../TruthConverters/Root/xAODtoHepMCTool.cxx | 385 ++++++++++-------- .../Root/HiggsTruthCategoryTool.cxx | 4 +- .../TruthCategoriesDecorator.h | 6 + .../python/TruthCategories.py | 1 - .../src/TruthCategoriesDecorator.cxx | 10 +- 5 files changed, 233 insertions(+), 173 deletions(-) diff --git a/Generators/TruthConverters/Root/xAODtoHepMCTool.cxx b/Generators/TruthConverters/Root/xAODtoHepMCTool.cxx index 2caacc61a22..3099309437d 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/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/DerivationFrameworkHiggs/TruthCategoriesDecorator.h b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/DerivationFrameworkHiggs/TruthCategoriesDecorator.h index 4965e9756b8..13881ddd0e6 100644 --- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/DerivationFrameworkHiggs/TruthCategoriesDecorator.h +++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/DerivationFrameworkHiggs/TruthCategoriesDecorator.h @@ -78,6 +78,12 @@ namespace DerivationFramework { 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{}; diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/TruthCategories.py b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/TruthCategories.py index 1e449333c35..2bf27e20fa9 100644 --- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/TruthCategories.py +++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/TruthCategories.py @@ -4,7 +4,6 @@ from DerivationFrameworkCore.DerivationFrameworkMaster import DerivationFramewor from AthenaCommon.GlobalFlags import globalflags if globalflags.DataSource()=='geant4': - import AthenaCommon.CfgMgr as CfgMgr from DerivationFrameworkHiggs.DerivationFrameworkHiggsConf import DerivationFramework__TruthCategoriesDecorator DFHTXSdecorator = DerivationFramework__TruthCategoriesDecorator(name = "DFHTXSdecorator") DerivationFrameworkJob += DFHTXSdecorator \ No newline at end of file diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx index cfc8308539a..ce663c1cdfc 100644 --- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx +++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx @@ -74,9 +74,11 @@ namespace DerivationFramework { } for (const std::string& key : decor_keys) m_STXSDecors.emplace_back(m_evtInfoKey.key() + "." + key); - for (const auto& p4_pair : m_p4_decors) { + for (auto& p4_pair : m_p4_decors) { std::vector keys = p4_pair.second.vect(); - for (EvtInfoDecorKey& key : keys) m_STXSDecors.push_back(std::move(key)); + 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()); @@ -208,11 +210,13 @@ namespace DerivationFramework { // ANGRY MESSAGE HERE ATH_MSG_FATAL("The HEP MC GenEvent conversion failed"); return StatusCode::FAILURE; + } else { + ATH_MSG_DEBUG("Found "< 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"}; -- GitLab From 457b16bb883cf9d23e5792a8359c8a53a747a058 Mon Sep 17 00:00:00 2001 From: Johannes Junggeburth Date: Tue, 30 Aug 2022 12:39:43 +0200 Subject: [PATCH 7/8] Address typo comments --- .../src/TruthCategoriesDecorator.cxx | 23 +++++++++++-------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx index ce663c1cdfc..defdaf93364 100644 --- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx +++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx @@ -232,31 +232,34 @@ namespace DerivationFramework { static const IntDecorator dec_stage1IdxPt30{"HTXS_Stage1_FineIndex_pTjet25"}; dec_stage1CatPt25(*eventInfo) = htxs->stage1_cat_pTjet25GeV; - dec_stage1CatPt30(*eventInfo) = htxs->stage1_cat_pTjet25GeV; - dec_stage1IdxPt25(*eventInfo) = HTXSstage1_to_HTXSstage1FineIndex(*htxs, th_type); - dec_stage1IdxPt30(*eventInfo) = HTXSstage1_to_HTXSstage1FineIndex(*htxs, th_type, true); + 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_pTjet30"}; - static const IntDecorator dec_stage1p2_IdxPt30{"HTXS_Stage1_2_FineIndex_pTjet25"}; + 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; - dec_stage1p2_IdxPt25(*eventInfo) = HTXSstage1_2_to_HTXSstage1_2_FineIndex(*htxs, th_type); - dec_stage1p2_IdxPt30(*eventInfo) = HTXSstage1_2_to_HTXSstage1_2_FineIndex(*htxs, th_type, true); + /// 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_pTjet30"}; - static const IntDecorator dec_stage1p2_Fine_IdxPt30{"HTXS_Stage1_2_Fine_FineIndex_pTjet25"}; + 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); - dec_stage1p2_Fine_IdxPt30(*eventInfo) = HTXSstage1_2_Fine_to_HTXSstage1_2_Fine_FineIndex(*htxs, th_type, true); + 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"}; -- GitLab From b19dab2cbb63dcea4103356c4df1a7b10c85fdde Mon Sep 17 00:00:00 2001 From: Johannes Josef Junggeburth Date: Tue, 30 Aug 2022 13:03:03 +0200 Subject: [PATCH 8/8] Alright sometimes is the local copy better than the browser editign --- .../src/TruthCategoriesDecorator.cxx | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx index defdaf93364..0913ff31f65 100644 --- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx +++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx @@ -228,8 +228,8 @@ namespace DerivationFramework { // 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_pTjet30"}; - static const IntDecorator dec_stage1IdxPt30{"HTXS_Stage1_FineIndex_pTjet25"}; + 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; @@ -257,7 +257,6 @@ namespace DerivationFramework { 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); 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); -- GitLab