diff --git a/Generators/TruthRivetTools/CMakeLists.txt b/Generators/TruthRivetTools/CMakeLists.txt index 0c7b36f7ec470d9190b1c077802b85b615d2bbc8..d72cc1c682b2397fef146282f1bb0644b7e96478 100644 --- a/Generators/TruthRivetTools/CMakeLists.txt +++ b/Generators/TruthRivetTools/CMakeLists.txt @@ -1,6 +1,4 @@ -################################################################################ -# Package: TruthRivetTools -################################################################################ +# Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration # Declare the name of the package: atlas_subdir( TruthRivetTools ) @@ -24,7 +22,7 @@ atlas_add_library( TruthRivetToolsLib ${FASTJET_INCLUDE_DIRS} ${ROOT_INCLUDE_DIRS} ${GSL_INCLUDE_DIRS} LINK_LIBRARIES ${RIVET_LIBRARIES} ${YODA_LIBRARIES} ${FASTJET_LIBRARIES} ${ROOT_LIBRARIES} ${GSL_LIBRARIES} AsgTools AtlasHepMCLib AtlasHepMCsearchLib - GenInterfacesLib ) + CxxUtils GenInterfacesLib ) atlas_add_component( TruthRivetTools src/components/*.cxx diff --git a/Generators/TruthRivetTools/Root/HiggsTruthCategoryTool.cxx b/Generators/TruthRivetTools/Root/HiggsTruthCategoryTool.cxx index 9c54893190cb53341fe9113757f029a5068fc650..f7c446c9d80f2e045a05f34c3f0ea523a7313f6c 100644 --- a/Generators/TruthRivetTools/Root/HiggsTruthCategoryTool.cxx +++ b/Generators/TruthRivetTools/Root/HiggsTruthCategoryTool.cxx @@ -47,14 +47,15 @@ StatusCode HiggsTruthCategoryTool :: finalize () { } HTXS::HiggsClassification* HiggsTruthCategoryTool :: getHiggsTruthCategoryObject (const HepMC::GenEvent& HepMCEvent, const HTXS::HiggsProdMode prodMode) const { - - static std::once_flag flag; - std::call_once(flag, [&]() { - higgsTemplateCrossSections->setHiggsProdMode(prodMode); - rivetAnaHandler->init(HepMCEvent); - }); + if ( !m_isInitialized.test_and_set() ) { + [&]() { + higgsTemplateCrossSections->setHiggsProdMode(prodMode); + rivetAnaHandler->init(HepMCEvent); + }(); + } // fill histos if flag is specified - if ( m_outHistos ) rivetAnaHandler->analyze(HepMCEvent); + // if ( m_outHistos ) rivetAnaHandler->analyze(HepMCEvent); + // get the category output object containing the template cross section category, // and Higgs, V-boson, jets 4-vectors const Rivet::HiggsClassification htxs_cat_rivet = higgsTemplateCrossSections->classifyEvent(HepMCEvent,prodMode); diff --git a/Generators/TruthRivetTools/TruthRivetTools/ATLAS_CHECK_THREAD_SAFETY b/Generators/TruthRivetTools/TruthRivetTools/ATLAS_CHECK_THREAD_SAFETY new file mode 100644 index 0000000000000000000000000000000000000000..fa6670d43142b3d90b5296313ed11edb46526712 --- /dev/null +++ b/Generators/TruthRivetTools/TruthRivetTools/ATLAS_CHECK_THREAD_SAFETY @@ -0,0 +1 @@ +Generators/TruthRivetTools diff --git a/Generators/TruthRivetTools/TruthRivetTools/HiggsTemplateCrossSections.h b/Generators/TruthRivetTools/TruthRivetTools/HiggsTemplateCrossSections.h index f7c8482a3e2c8496886f295d46d38334722fb2a8..fe077d860c9c9793405315161a53d02215b3db9b 100644 --- a/Generators/TruthRivetTools/TruthRivetTools/HiggsTemplateCrossSections.h +++ b/Generators/TruthRivetTools/TruthRivetTools/HiggsTemplateCrossSections.h @@ -22,6 +22,7 @@ #include "AtlasHepMC/GenVertex.h" #include "AtlasHepMC/GenParticle.h" #include "AtlasHepMC/GenRanges.h" +#include "CxxUtils/checker_macros.h" namespace Rivet { @@ -44,7 +45,7 @@ namespace Rivet { /// @{ /// follow a "propagating" particle and return its last instance - Particle getLastInstance(Particle ptcl) { + Particle getLastInstance(Particle ptcl) const { if ( ptcl.genParticle()->end_vertex() ) { if ( !hasChild(ptcl.genParticle(),ptcl.pid()) ) return ptcl; else return getLastInstance(ptcl.children()[0]); @@ -53,7 +54,7 @@ namespace Rivet { } /// @brief Whether particle p originate from any of the ptcls - bool originateFrom(const Particle& p, const Particles& ptcls ) { + bool originateFrom(const Particle& p, const Particles& ptcls ) const { auto prodVtx = p.genParticle()->production_vertex(); if (prodVtx == nullptr) return false; // for each ancestor, check if it matches any of the input particles @@ -66,33 +67,33 @@ namespace Rivet { } /// @brief Whether particle p originates from p2 - bool originateFrom(const Particle& p, const Particle& p2 ) { + bool originateFrom(const Particle& p, const Particle& p2 ) const { Particles ptcls = {p2}; return originateFrom(p,ptcls); } /// @brief Checks whether the input particle has a child with a given PDGID - bool hasChild(HepMC::ConstGenParticlePtr ptcl, int pdgID) { + bool hasChild(HepMC::ConstGenParticlePtr ptcl, int pdgID) const { for (const Particle& child:Particle(*ptcl).children()) if (child.pid()==pdgID) return true; return false; } /// @brief Checks whether the input particle has a parent with a given PDGID - bool hasParent(HepMC::ConstGenParticlePtr ptcl, int pdgID) { + bool hasParent(HepMC::ConstGenParticlePtr ptcl, int pdgID) const { for (auto parent:Rivet::HepMCUtils::particles(ptcl->production_vertex(),Relatives::PARENTS)) if (parent->pdg_id()==pdgID) return true; return false; } /// @brief Return true is particle decays to quarks - bool quarkDecay(const Particle &p) { + bool quarkDecay(const Particle &p) const { for (const Particle& child:p.children()) if (PID::isQuark(child.pid())) return true; return false; } /// @brief Return true if particle decays to charged leptons. - bool ChLeptonDecay(const Particle &p) { + bool ChLeptonDecay(const Particle &p) const { for (const Particle& child:p.children()) if (PID::isChLepton(child.pid())) return true; return false; @@ -101,7 +102,7 @@ namespace Rivet { /// @brief Returns the classification object with the error code set. /// Prints an warning message, and keeps track of number of errors HiggsClassification error(HiggsClassification &cat, HTXS::ErrorCode err, - std::string msg="", int NmaxWarnings=20) { + std::string msg="", int NmaxWarnings=20) const { // Set the error, and keep statistics cat.errorCode = err; ++m_errorCount[err]; @@ -116,10 +117,8 @@ namespace Rivet { /// @} /// @brief Main classificaion method. - HiggsClassification classifyEvent(const Event& event, const HTXS::HiggsProdMode prodMode ) { + HiggsClassification classifyEvent(const Event& event, const HTXS::HiggsProdMode prodMode ) const { - if (m_HiggsProdMode==HTXS::UNKNOWN) m_HiggsProdMode = prodMode; - // the classification object HiggsClassification cat; cat.prodMode = prodMode; @@ -310,7 +309,7 @@ namespace Rivet { /// @{ /// @brief Return bin index of x given the provided bin edges. 0=first bin, -1=underflow bin. - int getBin(double x, std::vector<double> bins) { + int getBin(double x, std::vector<double> bins) const { if (bins.size()==0||x<bins[0]) return -1; // should not happen! for (size_t i=1;i<bins.size();++i) if (x<bins[i]) return i-1; @@ -320,7 +319,7 @@ namespace Rivet { /// @brief VBF topolog selection /// 0 = fail loose selction: m_jj > 400 GeV and Dy_jj > 2.8 /// 1 pass loose, but fail additional cut pT(Hjj)<25. 2 pass tight selection - int vbfTopology(const Jets &jets, const Particle &higgs) { + int vbfTopology(const Jets &jets, const Particle &higgs) const { if (jets.size()<2) return 0; const FourMomentum &j1=jets[0].momentum(), &j2=jets[1].momentum(); bool VBFtopo = (j1+j2).mass() > 400.0 && std::abs(j1.rapidity()-j2.rapidity()) > 2.8; @@ -330,7 +329,7 @@ namespace Rivet { /// 0 = fail loose selection: m_jj > 350 GeV /// 1 pass loose, but fail additional cut pT(Hjj)<25. 2 pass pT(Hjj)>25 selection /// 3 pass tight (m_jj>700 GeV), but fail additional cut pT(Hjj)<25. 4 pass pT(Hjj)>25 selection - int vbfTopology_Stage1_2(const Jets &jets, const Particle &higgs) { + int vbfTopology_Stage1_2(const Jets &jets, const Particle &higgs) const { if (jets.size()<2) return 0; const FourMomentum &j1=jets[0].momentum(), &j2=jets[1].momentum(); double mjj = (j1+j2).mass(); @@ -344,7 +343,7 @@ namespace Rivet { /// 3 pass 700<m_jj<1000 GeV, but fail additional cut pT(Hjj)<25. 4 pass pT(Hjj)>25 selection /// 5 pass 1000<m_jj<1500 GeV, but fail additional cut pT(Hjj)<25. 6 pass pT(Hjj)>25 selection /// 7 pass m_jj>1500 GeV, but fail additional cut pT(Hjj)<25. 8 pass pT(Hjj)>25 selection - int vbfTopology_Stage1_2_Fine(const Jets &jets, const Particle &higgs) { + int vbfTopology_Stage1_2_Fine(const Jets &jets, const Particle &higgs) const { if (jets.size()<2) return 0; const FourMomentum &j1=jets[0].momentum(), &j2=jets[1].momentum(); double mjj = (j1+j2).mass(); @@ -356,12 +355,12 @@ namespace Rivet { } /// @brief Whether the Higgs is produced in association with a vector boson (VH) - bool isVH(HTXS::HiggsProdMode p) { return p==HTXS::WH || p==HTXS::QQ2ZH || p==HTXS::GG2ZH; } + bool isVH(HTXS::HiggsProdMode p) const { return p==HTXS::WH || p==HTXS::QQ2ZH || p==HTXS::GG2ZH; } /// @brief Stage-0 HTXS categorization HTXS::Stage0::Category getStage0Category(const HTXS::HiggsProdMode prodMode, const Particle &higgs, - const Particle &V) { + const Particle &V) const { using namespace HTXS::Stage0; int ctrlHiggs = std::abs(higgs.rapidity())<2.5; // Special cases first, qq→Hqq @@ -378,7 +377,7 @@ namespace Rivet { HTXS::Stage1::Category getStage1Category(const HTXS::HiggsProdMode prodMode, const Particle &higgs, const Jets &jets, - const Particle &V) { + const Particle &V) const { using namespace HTXS::Stage1; int Njets=jets.size(), ctrlHiggs = std::abs(higgs.rapidity())<2.5, fwdHiggs = !ctrlHiggs; double pTj1 = jets.size() ? jets[0].momentum().pt() : 0; @@ -444,7 +443,7 @@ namespace Rivet { HTXS::Stage1_2::Category getStage1_2_Category(const HTXS::HiggsProdMode prodMode, const Particle &higgs, const Jets &jets, - const Particle &V) { + const Particle &V) const { using namespace HTXS::Stage1_2; int Njets=jets.size(), ctrlHiggs = std::abs(higgs.rapidity())<2.5, fwdHiggs = !ctrlHiggs; int vbfTopo = vbfTopology_Stage1_2(jets,higgs); @@ -521,7 +520,7 @@ namespace Rivet { HTXS::Stage1_2_Fine::Category getStage1_2_Fine_Category(const HTXS::HiggsProdMode prodMode, const Particle &higgs, const Jets &jets, - const Particle &V) { + const Particle &V) const { using namespace HTXS::Stage1_2_Fine; int Njets=jets.size(), ctrlHiggs = std::abs(higgs.rapidity())<2.5, fwdHiggs = !ctrlHiggs; int vbfTopo = vbfTopology_Stage1_2_Fine(jets,higgs); @@ -766,7 +765,7 @@ namespace Rivet { private: double m_sumw; HTXS::HiggsProdMode m_HiggsProdMode; - std::map<HTXS::ErrorCode,size_t> m_errorCount; + mutable std::array<std::atomic<size_t>, HTXS::NUM_ERRORCODES> m_errorCount ATLAS_THREAD_SAFE {}; Histo1DPtr m_hist_stage0; Histo1DPtr m_hist_stage1_pTjet25, m_hist_stage1_pTjet30; Histo1DPtr m_hist_stage1_2_pTjet25, m_hist_stage1_2_pTjet30; diff --git a/Generators/TruthRivetTools/TruthRivetTools/HiggsTemplateCrossSectionsDefs.h b/Generators/TruthRivetTools/TruthRivetTools/HiggsTemplateCrossSectionsDefs.h index 10492c7aab971b9887d91e8a2f74a4982efec655..44453bf2c66e8355fbd6936c34b96781c3e14f44 100644 --- a/Generators/TruthRivetTools/TruthRivetTools/HiggsTemplateCrossSectionsDefs.h +++ b/Generators/TruthRivetTools/TruthRivetTools/HiggsTemplateCrossSectionsDefs.h @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration */ #ifndef TRUTHRIVETTOOLS_HIGGSTEMPLATECROSSSECTIONSDEFS_H @@ -19,7 +19,8 @@ namespace HTXS { HS_VTX_IDENTIFICATION = 5, ///< failed to identify hard scatter vertex VH_IDENTIFICATION = 6, ///< failed to identify associated vector boson VH_DECAY_IDENTIFICATION = 7, ///< failed to identify associated vector boson decay products - TOP_W_IDENTIFICATION = 8 ///< failed to identify top decay + TOP_W_IDENTIFICATION = 8, ///< failed to identify top decay + NUM_ERRORCODES ///< number of error codes (keep this unnumbered and last) }; /// Higgs production modes, corresponding to input sample diff --git a/Generators/TruthRivetTools/TruthRivetTools/HiggsTruthCategoryTool.h b/Generators/TruthRivetTools/TruthRivetTools/HiggsTruthCategoryTool.h index a8e4716e2b2d393e143e04eea79cc01e31d2ea0d..705a7aabdc3d2a83b7d3e437fc3d38f6b12cac77 100644 --- a/Generators/TruthRivetTools/TruthRivetTools/HiggsTruthCategoryTool.h +++ b/Generators/TruthRivetTools/TruthRivetTools/HiggsTruthCategoryTool.h @@ -11,6 +11,7 @@ #ifndef TRUTHRIVETTOOLS_HIGGSTRUTHCATEGORYTOOL_H #define TRUTHRIVETTOOLS_HIGGSTRUTHCATEGORYTOOL_H 1 +#include "CxxUtils/checker_macros.h" #include "Rivet/AnalysisHandler.hh" #include "TruthRivetTools/HiggsTemplateCrossSections.h" @@ -49,6 +50,7 @@ class HiggsTruthCategoryTool StatusCode finalize () override; HTXS::HiggsClassification* getHiggsTruthCategoryObject(const HepMC::GenEvent& HepMCEvent, const HTXS::HiggsProdMode prodMode) const override; private: + mutable std::atomic_flag m_isInitialized ATLAS_THREAD_SAFE = ATOMIC_FLAG_INIT; bool m_outHistos; };