diff --git a/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/AsgAnalysisAlgorithms/AsgAnalysisAlgorithmsDict.h b/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/AsgAnalysisAlgorithms/AsgAnalysisAlgorithmsDict.h index 9514cfab34b003978aff022e0da74047c152aa03..5a07aec0b3a51ede12c83e920784befa275713c7 100644 --- a/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/AsgAnalysisAlgorithms/AsgAnalysisAlgorithmsDict.h +++ b/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/AsgAnalysisAlgorithms/AsgAnalysisAlgorithmsDict.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 */ /// @author Nils Krumnack @@ -27,6 +27,7 @@ #include #include #include +#include #include #include #include diff --git a/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/AsgAnalysisAlgorithms/SysTruthWeightAlg.h b/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/AsgAnalysisAlgorithms/SysTruthWeightAlg.h new file mode 100644 index 0000000000000000000000000000000000000000..e6428ebe185a68f32101a84a84ad89c9c58ba926 --- /dev/null +++ b/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/AsgAnalysisAlgorithms/SysTruthWeightAlg.h @@ -0,0 +1,57 @@ +/* + Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration +*/ + +/// @author Miha Muskinja + +#ifndef ASG_ANALYSIS_ALGORITHMS__PMG_HF_PRODUCTION_FRACTION_ALG_H +#define ASG_ANALYSIS_ALGORITHMS__PMG_HF_PRODUCTION_FRACTION_ALG_H + +#include +#include +#include +#include +#include +#include + +namespace CP +{ +/// \brief an algorithm for calling \ref ISysTruthWeightTool + +class SysTruthWeightAlg final : public EL::AnaAlgorithm +{ + /// \brief the standard constructor +public: + SysTruthWeightAlg(const std::string& name, ISvcLocator* pSvcLocator); + +public: + StatusCode initialize() override; + +public: + StatusCode execute() override; + + /// \brief the truth particle container to use for the calculation +private: + std::string m_truthParticleContainer; + + /// \brief the tool +private: + ToolHandle m_sysTruthWeightTool; + + /// \brief the systematics list we run +private: + SysListHandle m_systematicsList{this}; + + /// \brief the event collection we run on +private: + SysReadHandle m_eventInfoHandle{ + this, "eventInfo", "EventInfo", "the event info object to run on"}; + + /// \brief the decoration for the truth weights +private: + SysWriteDecorHandle m_decoration{ + this, "decoration", "", "the decoration for the truth weights"}; +}; +} // namespace CP + +#endif diff --git a/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/AsgAnalysisAlgorithms/selection.xml b/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/AsgAnalysisAlgorithms/selection.xml index 1d8e58ad2376e48ebb4498bed41d2c1503956a1b..f50c3eacb3577567ba96b066bae1484eba746e64 100644 --- a/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/AsgAnalysisAlgorithms/selection.xml +++ b/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/AsgAnalysisAlgorithms/selection.xml @@ -20,6 +20,7 @@ + diff --git a/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/Root/SysTruthWeightAlg.cxx b/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/Root/SysTruthWeightAlg.cxx new file mode 100644 index 0000000000000000000000000000000000000000..c36011dfa0b258d246b8f9c2595632b7a3178d70 --- /dev/null +++ b/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/Root/SysTruthWeightAlg.cxx @@ -0,0 +1,69 @@ +/* + Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration +*/ + +/// @author Miha Muskinja + +// +// includes +// + +// EDM include(s): +#include "xAODTruth/TruthParticleContainer.h" + +// Local include(s): +#include + +// +// method implementations +// + +namespace CP +{ +SysTruthWeightAlg::SysTruthWeightAlg(const std::string &name, ISvcLocator *pSvcLocator) + : AnaAlgorithm(name, pSvcLocator), m_sysTruthWeightTool("PMGTools::PMGHFProductionFractionTool", this) +{ + // The truth particle container + declareProperty("TruthParticleContainer", m_truthParticleContainer = "TruthParticles"); + + // The production fraction rw tool + declareProperty("sysTruthWeightTool", m_sysTruthWeightTool); +} + +StatusCode SysTruthWeightAlg::initialize() +{ + if (m_decoration.empty()) + { + ANA_MSG_ERROR("no decoration name set"); + return StatusCode::FAILURE; + } + + ANA_CHECK(m_sysTruthWeightTool.retrieve()); + ANA_CHECK(m_eventInfoHandle.initialize(m_systematicsList, SG::AllowEmpty)); + ANA_CHECK(m_decoration.initialize(m_systematicsList, m_eventInfoHandle, SG::AllowEmpty)); + ANA_CHECK(m_systematicsList.addSystematics(*m_sysTruthWeightTool)); + ANA_CHECK(m_systematicsList.initialize()); + + return StatusCode::SUCCESS; +} + +StatusCode SysTruthWeightAlg::execute() +{ + + // Retreive the truth particles container + const xAOD::TruthParticleContainer *truthParticles = nullptr; + ANA_CHECK(evtStore()->retrieve(truthParticles, m_truthParticleContainer)); + + for (const auto &sys : m_systematicsList.systematicsVector()) + { + ANA_CHECK(m_sysTruthWeightTool->applySystematicVariation(sys)); + const xAOD::EventInfo *eventInfo = nullptr; + ANA_CHECK(m_eventInfoHandle.retrieve(eventInfo, sys)); + + m_decoration.set(*eventInfo, m_sysTruthWeightTool->getWeight(truthParticles), sys); + } + + // return + return StatusCode::SUCCESS; +} +} // namespace CP diff --git a/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/python/AsgAnalysisAlgorithmsTest.py b/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/python/AsgAnalysisAlgorithmsTest.py index 2f9a7ead3f5d3ed97eb1e16d8fe0cb638135e31d..002af543b71264aaf6b696f5a3f844e307527c66 100644 --- a/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/python/AsgAnalysisAlgorithmsTest.py +++ b/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/python/AsgAnalysisAlgorithmsTest.py @@ -1,4 +1,4 @@ -# Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration +# Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration # # @author Nils Krumnack # @author Tadej Novak @@ -205,7 +205,7 @@ def makeGeneratorAlgorithmsSequence (dataType) : # Include, and then set up the generator analysis sequence: from AsgAnalysisAlgorithms.GeneratorAnalysisSequence import \ makeGeneratorAnalysisSequence - generatorSequence = makeGeneratorAnalysisSequence( dataType, saveCutBookkeepers=True, runNumber=284500, cutBookkeepersSystematics=True ) + generatorSequence = makeGeneratorAnalysisSequence( dataType, saveCutBookkeepers=True, runNumber=284500, cutBookkeepersSystematics=True, prodFractionWeight=True ) algSeq += generatorSequence # Set up an ntuple to check the job with: @@ -218,6 +218,7 @@ def makeGeneratorAlgorithmsSequence (dataType) : 'EventInfo.runNumber -> runNumber', 'EventInfo.eventNumber -> eventNumber', 'EventInfo.generatorWeight_%SYS% -> generatorWeight_%SYS%', + 'EventInfo.prodFracWeight_%SYS% -> prodFracWeight_%SYS%', ] algSeq += ntupleMaker treeFiller = createAlgorithm( 'CP::TreeFillerAlg', 'TreeFiller' ) diff --git a/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/python/GeneratorAnalysisSequence.py b/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/python/GeneratorAnalysisSequence.py index 244c9cb25e4ddf3555cc04b4078367f0a385dfc3..36c152d9f7c45d5b528140e56cafbf8a039a8081 100644 --- a/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/python/GeneratorAnalysisSequence.py +++ b/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/python/GeneratorAnalysisSequence.py @@ -1,4 +1,4 @@ -# Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration +# Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration # AnaAlgorithm import(s): from AnaAlgorithm.AnaAlgSequence import AnaAlgSequence @@ -7,7 +7,9 @@ from AnaAlgorithm.DualUseConfig import createAlgorithm, addPrivateTool def makeGeneratorAnalysisSequence( dataType, saveCutBookkeepers=False, runNumber=0, - cutBookkeepersSystematics=False ): + cutBookkeepersSystematics=False, + generator="default", + prodFractionWeight=False ): """Create a generator analysis algorithm sequence Keyword arguments: @@ -15,6 +17,8 @@ def makeGeneratorAnalysisSequence( dataType, saveCutBookkeepers -- save cut bokkeepers information into output file runNumber -- MC run number cutBookkeepersSystematics -- store CutBookkeepers systematics + generator -- Generator for HF production fraction weights + prodFractionWeight -- Calculate the HF production fraction weights """ if dataType not in ["mc", "afii"] : @@ -23,6 +27,27 @@ def makeGeneratorAnalysisSequence( dataType, if saveCutBookkeepers and not runNumber: raise ValueError ("invalid run number: " + 0) + if generator not in ["default", "Pythia8", "Sherpa221", "Sherpa228", "Sherpa2210", "Herwig7", "Herwig713", "Herwig721", "amc@NLO"]: + raise ValueError ("invalid generator type: " + generator) + + # MC/MC scale factors configuration + if generator == "default": + DSID = "410470" + elif generator == "Sherpa221": + DSID = "410250" + elif generator == "Sherpa228": + DSID = "421152" + elif generator == "Sherpa2210": + DSID = "700122" + elif generator == "Herwig7": + DSID = "410558" + elif generator == "Herwig713": + DSID = "411233" + elif generator == "Herwig721": + DSID = "600666" + elif generator == "amc@NLO": + DSID = "410464" + # Create the analysis algorithm sequence object: seq = AnaAlgSequence( "GeneratorAnalysisSequence" ) @@ -38,8 +63,15 @@ def makeGeneratorAnalysisSequence( dataType, alg = createAlgorithm( 'CP::PMGTruthWeightAlg', 'PMGTruthWeightAlg' ) addPrivateTool( alg, 'truthWeightTool', 'PMGTools::PMGTruthWeightTool' ) alg.decoration = 'generatorWeight_%SYS%' - seq.append( alg, inputPropName = None ) + # Production fraction weights + if prodFractionWeight: + alg = createAlgorithm('CP::SysTruthWeightAlg', 'SysTruthWeightAlg') + addPrivateTool(alg, 'sysTruthWeightTool', 'PMGTools::PMGHFProductionFractionTool') + alg.decoration = 'prodFracWeight_%SYS%' + alg.sysTruthWeightTool.ShowerGenerator = DSID + seq.append( alg, inputPropName = 'eventInfo') + # Return the sequence: return seq diff --git a/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/share/GeneratorAlgorithmsTest_eljob.py b/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/share/GeneratorAlgorithmsTest_eljob.py index 9406cda8fdfe2eb90d169e4452e0eb3f7adb482e..996bf791931bc0100d9722ca8f1c81243f6ef384 100755 --- a/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/share/GeneratorAlgorithmsTest_eljob.py +++ b/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/share/GeneratorAlgorithmsTest_eljob.py @@ -1,6 +1,6 @@ #!/usr/bin/env python # -# Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration +# Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration # # @author Tadej Novak diff --git a/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/src/components/AsgAnalysisAlgorithms_entries.cxx b/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/src/components/AsgAnalysisAlgorithms_entries.cxx index cdeed836fb66c39a6d94595985ca7a4ad8fac09d..95d84cf630e7f78a2107fb57d122767344e1a092 100644 --- a/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/src/components/AsgAnalysisAlgorithms_entries.cxx +++ b/PhysicsAnalysis/Algorithms/AsgAnalysisAlgorithms/src/components/AsgAnalysisAlgorithms_entries.cxx @@ -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 */ /// @author Nils Krumnack @@ -26,6 +26,7 @@ #include #include #include +#include #include #include #include @@ -49,6 +50,7 @@ DECLARE_NAMESPACE_ALGORITHM_FACTORY (CP, ObjectCutFlowHistAlg) DECLARE_NAMESPACE_ALGORITHM_FACTORY (CP, OverlapRemovalAlg) DECLARE_NAMESPACE_ALGORITHM_FACTORY (CP, PileupReweightingAlg) DECLARE_NAMESPACE_ALGORITHM_FACTORY (CP, PMGTruthWeightAlg) +DECLARE_NAMESPACE_ALGORITHM_FACTORY (CP, SysTruthWeightAlg) DECLARE_NAMESPACE_ALGORITHM_FACTORY (CP, SysListDumperAlg) DECLARE_NAMESPACE_ALGORITHM_FACTORY (CP, TreeFillerAlg) DECLARE_NAMESPACE_ALGORITHM_FACTORY (CP, TreeMakerAlg) @@ -72,6 +74,7 @@ DECLARE_FACTORY_ENTRIES(AsgAnalysisAlgorithms) { DECLARE_NAMESPACE_ALGORITHM (CP, OverlapRemovalAlg) DECLARE_NAMESPACE_ALGORITHM (CP, PileupReweightingAlg) DECLARE_NAMESPACE_ALGORITHM (CP, PMGTruthWeightAlg) + DECLARE_NAMESPACE_ALGORITHM (CP, SysTruthWeightAlg) DECLARE_NAMESPACE_ALGORITHM (CP, SysListDumperAlg) DECLARE_NAMESPACE_ALGORITHM (CP, TreeFillerAlg) DECLARE_NAMESPACE_ALGORITHM (CP, TreeMakerAlg) diff --git a/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/PMGHFProductionFractionTool.h b/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/PMGHFProductionFractionTool.h new file mode 100644 index 0000000000000000000000000000000000000000..a9679ebe5e8c672af229469b9e98d7508dbd6a81 --- /dev/null +++ b/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/PMGHFProductionFractionTool.h @@ -0,0 +1,130 @@ +/* + Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef PMGTOOLS_PMGHFPRODUCTIONFRACTIONTOOL_H +#define PMGTOOLS_PMGHFPRODUCTIONFRACTIONTOOL_H + +// EDM include(s): +#include "AsgTools/AsgTool.h" +#include "PATInterfaces/SystematicsCache.h" + +// Interface include(s): +#include "PMGAnalysisInterfaces/ISysTruthWeightTool.h" + +namespace PMGTools +{ +/// Interface for xAOD Truth Weight Tool which retrieves +/// the charm / bottom hadron content of an event and derives +/// an event weight based on the HF production fractions. +/// +/// @author Miha Muskinja +/// +class PMGHFProductionFractionTool : public virtual ISysTruthWeightTool, + public asg::AsgTool +{ + /// Create a proper constructor for Athena + ASG_TOOL_CLASS(PMGHFProductionFractionTool, PMGTools::ISysTruthWeightTool) + +public: + /// Create a constructor for standalone usage + PMGHFProductionFractionTool(const std::string &name); + + /// @name Function(s) implementing the asg::IAsgTool interface + /// @{ + + /// Function initialising the tool + virtual StatusCode initialize() override; + + /// @} + + /// @name Function(s) implementing the ISysTruthWeightTool interface + /// @{ + + /// Implements interface from ISysTruthWeightTool + virtual float getWeight(const xAOD::TruthParticleContainer *truthParticles) const override; + + /// Which systematics have an effect on the tool's behaviour? + virtual CP::SystematicSet affectingSystematics() const override; + + /// Use specific systematic + virtual CP::SystematicCode applySystematicVariation(const CP::SystematicSet &systConfig) override; + + /// Copied from SystematicsTool.cxx + virtual CP::SystematicSet recommendedSystematics() const override { return affectingSystematics(); }; + + /// Copied from SystematicsTool.cxx + virtual bool isAffectedBySystematic(const CP::SystematicVariation &systematic) const override + { + return (affectingSystematics().find(systematic) != affectingSystematics().end()); + }; + + /// @} + +private: + /// Setup weights + StatusCode setupProductionFractions(); + + /// Read production fractions from input file + StatusCode readProductionFractionsFile(std::string, std::map> *); + + /// Checks if a particle originates from a bottom decay + bool fromBdecay(const xAOD::TruthParticle *particle) const; + + /// Loops back the decay chain through particles with the same pdgId (e.g. photon emission) + const xAOD::TruthParticle *getInitialParticle(const xAOD::TruthParticle *tp) const; + + /// Print the current production fractions + void printCurrentProdFractions() const; + + /// Systematics set of the weight systematics + CP::SystematicSet m_systematicsSet; + + /// Path to calibration area + std::string m_calibrationAreaPath; + + /// MC Shower generator map + std::map m_showerGeneratorMap; + + /// MC Shower generator map file name + std::string m_showerGeneratorMapFile; + + /// MC Shower generator software (valid options: Pythia8) + std::string m_showerGenerator; + + /// Input file with charm production fractions + std::string m_charmFilename; + + /// Input file with bottom production fractions + std::string m_bottomFilename; + + /// The fiducial charm/bottom pT cut (in GeV) + float m_fiducialPtCut; + + /// The fiducial charm/bottom eta cut + float m_fiducialEtaCut; + + /// Charm production fraction weights + std::map> m_charmProdFractionWeights; + + /// Bottom production fraction weights + std::map> m_bottomProdFractionWeights; + + /// Struct for the production fractions + struct ParameterSet { + std::map charmWeights; + std::map bottomWeights; + }; + + /// calculate the parameter set for the given systematic + StatusCode setSystematicVariation(const CP::SystematicSet& systConfig, ParameterSet& param) const; + + /// The SystematicsCache object + CP::SystematicsCache m_Parameters{this}; + + /// Current production fractions (depending on applied systematic) + const ParameterSet* m_currentParameters{nullptr}; +}; +} // namespace PMGTools + +#endif // PMGTOOLS_PMGHFProductionFractionTool_H diff --git a/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/PMGToolsDict.h b/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/PMGToolsDict.h index 5066d2443fe33c4049bf2691302086a69439887f..d79b320badac12bee89bd13ed85969582d3e6173 100644 --- a/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/PMGToolsDict.h +++ b/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/PMGToolsDict.h @@ -1,7 +1,7 @@ // Dear emacs, this is -*- c++ -*- /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration */ // $Id: PMGToolsDict.h 762199 2016-07-15 13:47:19Z aknue $ @@ -13,5 +13,6 @@ #include "PMGTools/PMGDecayProductsSelectionTool.h" #include "PMGTools/PMGSherpa22VJetsWeightTool.h" #include "PMGTools/PMGTruthWeightTool.h" +#include "PMGTools/PMGHFProductionFractionTool.h" #endif // PMGTOOLS_PMGTOOLSDICT_H diff --git a/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/selection.xml b/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/selection.xml index 7141c2b611a64d6824a1bae5f8cc9c336eb40daa..7b384ec099cc1347abe3fee20546b34c8b4e8ec8 100644 --- a/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/selection.xml +++ b/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/selection.xml @@ -6,5 +6,6 @@ + diff --git a/PhysicsAnalysis/AnalysisCommon/PMGTools/Root/PMGHFProductionFractionTool.cxx b/PhysicsAnalysis/AnalysisCommon/PMGTools/Root/PMGHFProductionFractionTool.cxx new file mode 100644 index 0000000000000000000000000000000000000000..a22e635e05d84de0c9d95251b4638bdb2f664a9d --- /dev/null +++ b/PhysicsAnalysis/AnalysisCommon/PMGTools/Root/PMGHFProductionFractionTool.cxx @@ -0,0 +1,542 @@ +/* + Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration +*/ + +// $Id: PMGHFProductionFractionTool.cxx + +// Systematics include(s): +#include "PathResolver/PathResolver.h" +#include "PATInterfaces/SystematicRegistry.h" + +// Local include(s): +#include "PMGTools/PMGHFProductionFractionTool.h" + +// C++ +#include +#include +#include + +namespace fs = boost::filesystem; + +namespace PMGTools +{ +PMGHFProductionFractionTool::PMGHFProductionFractionTool(const std::string &name): + asg::AsgTool(name), + m_calibrationAreaPath(""), + m_showerGeneratorMap({}), + m_showerGeneratorMapFile(""), + m_showerGenerator(""), + m_charmFilename(""), + m_bottomFilename(""), + m_charmProdFractionWeights({}), + m_bottomProdFractionWeights({}) +{ + // Calibration area path + declareProperty("CalibrationAreaPath", m_calibrationAreaPath = "xAODBTaggingEfficiency/13TeV/HFReweighting-2022-04-19_v1", + "Path to the calibration area"); + + // MC Shower generator map + declareProperty("ShowerGeneratorMapFile", m_showerGeneratorMapFile = "HFProductionFractionsMap.txt", + "Input file name for the MC shower generator map"); + + // MC Shower generator software + declareProperty("ShowerGenerator", m_showerGenerator = "", + "MC Shower generator software"); + + // The fiducial charm/bottom pT cut (in MeV) + declareProperty("FiducialPtCut", m_fiducialPtCut = 5000, + "The fiducial charm/bottom pT cut. The recommended value is 5000 (MeV)."); + + // The fiducial charm/bottom eta cut + declareProperty("FiducialEtaCut", m_fiducialEtaCut = 2.5, + "The fiducial charm/bottom eta cut. The recommended value is 2.5."); +} + +StatusCode PMGHFProductionFractionTool::initialize() +{ + // Tell the user what's happening: + ATH_MSG_INFO("Initializing " << name() << "..."); + + // Check for correct settings + if (m_showerGenerator.empty()) + { + ATH_MSG_WARNING("The property `ShowerGenerator' was not set!"); + return StatusCode::FAILURE; + } + + // read input files and calculate weights + ANA_CHECK(setupProductionFractions()); + + // setup the SystematicsCache object + m_Parameters.initialize(affectingSystematics(), [this](const CP::SystematicSet &systConfig, ParameterSet ¶m) + { return setSystematicVariation(systConfig, param); }); + + // Set default for running without systematics + ANA_CHECK(applySystematicVariation(CP::SystematicSet())); + + // Register systematics with the registry + CP::SystematicRegistry ®istry = CP::SystematicRegistry::getInstance(); + if (registry.registerSystematics(*this) != CP::SystematicCode::Ok) + { + ATH_MSG_ERROR("Unkown systematic list"); + return StatusCode::FAILURE; + } + + // print + printCurrentProdFractions(); + + // Return gracefully + return StatusCode::SUCCESS; +} + +StatusCode PMGHFProductionFractionTool::readProductionFractionsFile(std::string filename, std::map> *weights) +{ + /* Example file structure: + # Charm production fraction from literature (Eur. Phys. J. C (2016) 76:397) + 411 421 431 4000 + NOSYS 0.2404 0.6086 0.0802 0.0623 + PROD_FRAC_CHARM_EIG_1 0.2405 0.6087 0.0770 0.0650 + PROD_FRAC_CHARM_EIG_2 0.2439 0.6108 0.0779 0.0593 + PROD_FRAC_CHARM_EIG_3 0.2461 0.6013 0.0808 0.0632 + */ + std::ifstream infile(filename); + std::string line; + int pdg1 = -1, pdg2 = -1, pdg3 = -1, pdg4 = -1; + while (std::getline(infile, line)) + { + std::istringstream iss(line); + if (line.rfind("#", 0) == 0 || line.empty()) + { + continue; + } + if (pdg1 < 0) + { + if (!(iss >> pdg1 >> pdg2 >> pdg3 >> pdg4)) + { + ATH_MSG_ERROR("Invalid formatting: " << line); + return StatusCode::FAILURE; + } + } + else + { + std::string sys; + float f1, f2, f3, f4; + if (!(iss >> sys >> f1 >> f2 >> f3 >> f4)) + { + ATH_MSG_ERROR("Invalid formatting: " << line); + return StatusCode::FAILURE; + } + if (sys == "NOSYS") + { + sys = ""; + } + if (affectingSystematics().find(CP::SystematicVariation(sys)) == affectingSystematics().end()) + { + ATH_MSG_ERROR("Systematic uncertaitny not recognized: " << sys); + return StatusCode::FAILURE; + } + weights->insert({CP::SystematicVariation(sys), {{pdg1, f1}, {pdg2, f2}, {pdg3, f3}, {pdg4, f4}}}); + } + } + + // Close file + infile.close(); + + // return + return StatusCode::SUCCESS; +} + +StatusCode PMGHFProductionFractionTool::setupProductionFractions() +{ + // Step 0: find available files + // ________________________________________________________________ + + // Charm production fractions data file + std::string mapFilename = PathResolverFindCalibFile((fs::path(m_calibrationAreaPath) / fs::path(m_showerGeneratorMapFile)).string()); + if (mapFilename.empty()) + { + ATH_MSG_ERROR("Input file not found " << m_showerGeneratorMapFile); + return StatusCode::FAILURE; + } + ATH_MSG_INFO("Found input file " << mapFilename); + + // Find all available MC generators + std::ifstream infile(mapFilename); + std::string line; + while (std::getline(infile, line)) + { + std::istringstream iss(line); + if (line.rfind("#", 0) == 0 || line.empty()) + { + continue; + } + + std::string id, name; + if (!(iss >> id >> name)) + { + ATH_MSG_ERROR("Invalid formatting: " << line); + return StatusCode::FAILURE; + } + if (id == "4") + { + m_charmFilename = name; + } + else if (id == "5") + { + m_bottomFilename = name; + } + else + { + m_showerGeneratorMap[id] = name; + } + } + + // Close file + infile.close(); + + // Check if the requested MC shower version is available + if (m_showerGeneratorMap.find(m_showerGenerator) == m_showerGeneratorMap.end()) + { + ATH_MSG_ERROR("MC shower generator " << m_showerGenerator << " not found!"); + ATH_MSG_INFO("Available generators are:"); + for (auto &gen : m_showerGeneratorMap) + { + ATH_MSG_INFO(gen.first); + } + return StatusCode::FAILURE; + } + + // Step 1: read HF production fractions from literature + // ________________________________________________________________ + + // Charm production fractions data file + std::string charmFilename = PathResolverFindCalibFile((fs::path(m_calibrationAreaPath) / fs::path(m_charmFilename)).string()); + if (charmFilename.empty()) + { + ATH_MSG_ERROR("Input file not found " << m_charmFilename); + return StatusCode::FAILURE; + } + ATH_MSG_INFO("Found input file " << charmFilename); + + // Bottom production fractions data file + std::string bottomFilename = PathResolverFindCalibFile((fs::path(m_calibrationAreaPath) / fs::path(m_bottomFilename)).string()); + if (bottomFilename.empty()) + { + ATH_MSG_ERROR("Input file not found " << m_bottomFilename); + return StatusCode::FAILURE; + } + ATH_MSG_INFO("Found input file " << bottomFilename); + + // Configure the production fractions from files + ATH_CHECK(readProductionFractionsFile(charmFilename, &m_charmProdFractionWeights)); + ATH_CHECK(readProductionFractionsFile(bottomFilename, &m_bottomProdFractionWeights)); + + // Step 2: read HF production fractions for a given MC shower + // ________________________________________________________________ + + // Read the config files from cvmfs + std::string filename = PathResolverFindCalibFile((fs::path(m_calibrationAreaPath) / fs::path(m_showerGeneratorMap[m_showerGenerator])).string()); + if (filename.empty()) + { + ATH_MSG_ERROR("Input file not found " << m_showerGeneratorMap[m_showerGenerator]); + return StatusCode::FAILURE; + } + ATH_MSG_INFO("Found input file " << filename); + + // Configure the weights from the file + /* Example file structure: + #Sherpa(v2.2.1) (DSID 410250) + #PDG ID | Production Fraction + 511 0.27243 + 521 0.273029 + 531 0.0893853 + 5000 0.365034 + 411 0.142884 + 421 0.379585 + 431 0.113005 + 4000 0.364526 + */ + infile = std::ifstream(filename); + std::map charm; + std::map bottom; + while (std::getline(infile, line)) + { + std::istringstream iss(line); + if (line.rfind("#", 0) == 0 || line.empty()) + { + continue; + } + + int pdgId; + float fraction; + if (!(iss >> pdgId >> fraction)) + { + ATH_MSG_ERROR("Invalid formatting: " << line); + return StatusCode::FAILURE; + } + if (pdgId == 4000 || pdgId < 500) + { + ATH_MSG_DEBUG("Charm meson " << pdgId << " has a fraction of " << fraction); + charm[pdgId] = fraction; + } + else + { + ATH_MSG_DEBUG("Bottom meson " << pdgId << " has a fraction of " << fraction); + bottom[pdgId] = fraction; + } + } + + // Close file + infile.close(); + + // Step 3: calcualte the weight as `w = f(literature) / f(MC)' + // ________________________________________________________________ + + // charm weights + for (auto &sys : m_charmProdFractionWeights) + { + for (auto &kv : charm) + { + m_charmProdFractionWeights[sys.first][kv.first] /= kv.second; + } + } + + // bottom weights + for (auto &sys : m_bottomProdFractionWeights) + { + for (auto &kv : bottom) + { + m_bottomProdFractionWeights[sys.first][kv.first] /= kv.second; + } + } + + // Return gracefully + return StatusCode::SUCCESS; +} + +bool PMGHFProductionFractionTool::fromBdecay(const xAOD::TruthParticle *particle) const +{ + for (unsigned int i = 0; i < particle->nParents(); i++) + { + if (particle->parent(i)->isHeavyHadron()) + { + if (particle->parent(i)->isBottomHadron()) + { + return true; + } + else + { + return fromBdecay(particle->parent(i)); + } + } + } + return false; +} + +const xAOD::TruthParticle *PMGHFProductionFractionTool::getInitialParticle(const xAOD::TruthParticle *tp) const +{ + if (tp->nParents() < 1) + { + return tp; + } + + const xAOD::TruthParticle *out = nullptr; + for (unsigned int i = 0; i < tp->nParents(); i++) + { + if (tp->parent(i) && (tp->parent(i)->absPdgId() == tp->absPdgId())) + { + const xAOD::TruthParticle *initialParticle = getInitialParticle(tp->parent(i)); + if (initialParticle && out && (initialParticle != out)) + { + // this should never happen, but if it does crash the program + throw std::runtime_error("Contradictory information in the truth decay chain!"); + } + out = initialParticle; + } + } + + if (out) + { + return out; + } + return tp; +} + +float PMGHFProductionFractionTool::getWeight(const xAOD::TruthParticleContainer *truthParticles) const +{ + // Clear vectors + std::set bottomHadrons; + std::set charmHadrons; + + // Find truth particles hadrons + for (const xAOD::TruthParticle *tp : *truthParticles) + { + if (!(tp->status() == 2 || tp->status() == 1)) + { + continue; + } + if (tp->isBottomHadron()) + { + if (tp->pt() >= m_fiducialPtCut && std::abs(tp->eta()) < m_fiducialEtaCut) + { + bottomHadrons.insert(getInitialParticle(tp)); + } + } + if (tp->isCharmHadron()) + { + if (!fromBdecay(getInitialParticle(tp)) && + tp->pt() >= m_fiducialPtCut && std::abs(tp->eta()) < m_fiducialEtaCut) + { + charmHadrons.insert(getInitialParticle(tp)); + } + } + } + + // The weight + float weight = 1.0; + + // Calculate the charm weight + for (auto *tp : charmHadrons) + { + int pdgId = tp->absPdgId(); + float w = 1.0; + if (tp->isCharmMeson() && m_currentParameters->charmWeights.find(pdgId) != m_currentParameters->charmWeights.end()) + { + w = m_currentParameters->charmWeights.at(pdgId); + } + else if (pdgId == 4122 || // Lambda_c+ + pdgId == 4232 || // Xi_c+ + pdgId == 4132 || // Xi_c0 + pdgId == 4332) // Omega_c0 + { + w = m_currentParameters->charmWeights.at(4000); + } + ATH_MSG_DEBUG("Charm weight for pdgId " << pdgId << ", pT " << tp->pt() << ", eta " << tp->eta() << ": " << w); + weight *= w; + } + + // Calculate the bottom weight + for (auto *tp : bottomHadrons) + { + int pdgId = tp->absPdgId(); + float w = 1.0; + if (tp->isBottomMeson() && m_currentParameters->bottomWeights.find(pdgId) != m_currentParameters->bottomWeights.end()) + { + w = m_currentParameters->bottomWeights.at(pdgId); + } + else if (pdgId == 5122 || // Lambda_b0 + pdgId == 5132 || // Xi_b- + pdgId == 5232 || // Xi_b0 + pdgId == 5332) // Omega_b- + { + w = m_currentParameters->bottomWeights.at(5000); + } + ATH_MSG_DEBUG("Bottom weight for pdgId " << pdgId << ", pT " << tp->pt() << ", eta " << tp->eta() << ": " << w); + weight *= w; + } + + return weight; +} + +void PMGHFProductionFractionTool::printCurrentProdFractions() const +{ + // MC Shower + ATH_MSG_INFO("The tool is currently configured to give weights for " << m_showerGenerator); + + // Charm production fraction weights + ATH_MSG_INFO("Currently set weights for charm production fractions:"); + for (auto const &generator : m_charmProdFractionWeights) + { + ATH_MSG_INFO("sys: " << generator.first.name()); + for (auto const &kv : generator.second) + { + ATH_MSG_INFO(" - " << kv.first << " " << kv.second); + } + } + + // Bottom production fraction weights + ATH_MSG_INFO("Currently set weights for bottom production fractions:"); + for (auto const &generator : m_bottomProdFractionWeights) + { + ATH_MSG_INFO("sys: " << generator.first.name()); + for (auto const &kv : generator.second) + { + ATH_MSG_INFO(" - " << kv.first << " " << kv.second); + } + } +} + +StatusCode PMGHFProductionFractionTool::setSystematicVariation(const CP::SystematicSet &systConfig, ParameterSet ¶m) const +{ + // retreive the previosuly calcualted charm weights + if (m_charmProdFractionWeights.find(CP::SystematicVariation("")) == m_charmProdFractionWeights.end()) + { + ATH_MSG_ERROR("Nominal charm weights not found"); + return StatusCode::FAILURE; + } + auto charmWeights = m_charmProdFractionWeights.at(CP::SystematicVariation("")); + if (systConfig.name().rfind("PROD_FRAC_CHARM", 0) == 0) + { + const CP::SystematicVariation ¤t_sys = *systConfig.begin(); + if (m_charmProdFractionWeights.find(current_sys) == m_charmProdFractionWeights.end()) + { + ATH_MSG_ERROR("Unknown systematic uncertainty " << systConfig.name()); + return StatusCode::FAILURE; + } + charmWeights = m_charmProdFractionWeights.at(current_sys); + } + + // retreive the previosuly calcualted bottom weights + if (m_bottomProdFractionWeights.find(CP::SystematicVariation("")) == m_bottomProdFractionWeights.end()) + { + ATH_MSG_ERROR("Nominal charm weights not found"); + return StatusCode::FAILURE; + } + auto bottomWeights = m_bottomProdFractionWeights.at(CP::SystematicVariation("")); + if (systConfig.name().rfind("PROD_FRAC_BOTTOM", 0) == 0) + { + const CP::SystematicVariation ¤t_sys = *systConfig.begin(); + if (m_bottomProdFractionWeights.find(current_sys) == m_bottomProdFractionWeights.end()) + { + ATH_MSG_ERROR("Unknown systematic uncertainty " << systConfig.name()); + return StatusCode::FAILURE; + } + bottomWeights = m_bottomProdFractionWeights.at(current_sys); + } + + // set the parameters + param.charmWeights = charmWeights; + param.bottomWeights = bottomWeights; + + return StatusCode::SUCCESS; +} + +CP::SystematicSet PMGHFProductionFractionTool::affectingSystematics() const +{ + CP::SystematicSet result; + + // three eigenvectors for charm + result.insert(CP::SystematicVariation("PROD_FRAC_CHARM_EIG_1")); + result.insert(CP::SystematicVariation("PROD_FRAC_CHARM_EIG_2")); + result.insert(CP::SystematicVariation("PROD_FRAC_CHARM_EIG_3")); + + // two eigenvectors for bottom (f(B0) = f(B+)) + result.insert(CP::SystematicVariation("PROD_FRAC_BOTTOM_EIG_1")); + result.insert(CP::SystematicVariation("PROD_FRAC_BOTTOM_EIG_2")); + + return result; +} + +CP::SystematicCode PMGHFProductionFractionTool::applySystematicVariation(const CP::SystematicSet &systConfig) +{ + // Only a single systematic can be applied at a time + if (systConfig.size() > 1) + { + return CP::SystematicCode::Unsupported; + } + + // Otherwise let's assume that things are correct + return m_Parameters.get(systConfig, m_currentParameters); +} + +} // namespace PMGTools diff --git a/PhysicsAnalysis/AnalysisCommon/PMGTools/src/components/PMGTools_entries.cxx b/PhysicsAnalysis/AnalysisCommon/PMGTools/src/components/PMGTools_entries.cxx index 9ea3ca34d935297496078f3c3a3454404c0c95d4..e216ac53d5cd6a27283b6642425b91b95e0f1646 100644 --- a/PhysicsAnalysis/AnalysisCommon/PMGTools/src/components/PMGTools_entries.cxx +++ b/PhysicsAnalysis/AnalysisCommon/PMGTools/src/components/PMGTools_entries.cxx @@ -1,4 +1,8 @@ +/* + Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration +*/ + #include "GaudiKernel/DeclareFactoryEntries.h" #include "PMGTools/PMGCrossSectionTool.h" @@ -6,6 +10,7 @@ #include "PMGTools/PMGSherpa22VJetsWeightTool.h" #include "PMGTools/PMGSherpaVjetsSysTool.h" #include "PMGTools/PMGTruthWeightTool.h" +#include "PMGTools/PMGHFProductionFractionTool.h" using namespace PMGTools; @@ -14,6 +19,7 @@ DECLARE_TOOL_FACTORY( PMGDecayProductsSelectionTool ) DECLARE_TOOL_FACTORY( PMGSherpa22VJetsWeightTool ) DECLARE_TOOL_FACTORY( PMGSherpaVjetsSysTool ) DECLARE_TOOL_FACTORY( PMGTruthWeightTool ) +DECLARE_TOOL_FACTORY( PMGHFProductionFractionTool ) DECLARE_FACTORY_ENTRIES( PMGTools ) { @@ -22,4 +28,5 @@ DECLARE_FACTORY_ENTRIES( PMGTools ) DECLARE_TOOL( PMGSherpa22VJetsWeightTool ); DECLARE_TOOL( PMGSherpaVjetsSysTool ); DECLARE_TOOL( PMGTruthWeightTool ); + DECLARE_TOOL( PMGHFProductionFractionTool ); } diff --git a/PhysicsAnalysis/Interfaces/PMGAnalysisInterfaces/CMakeLists.txt b/PhysicsAnalysis/Interfaces/PMGAnalysisInterfaces/CMakeLists.txt index 9bca6b5e8fec958a0851efe4fadef89ab6416e0a..a16dcda657be5e137a641349cf11a3c38bd5fbf9 100644 --- a/PhysicsAnalysis/Interfaces/PMGAnalysisInterfaces/CMakeLists.txt +++ b/PhysicsAnalysis/Interfaces/PMGAnalysisInterfaces/CMakeLists.txt @@ -1,3 +1,5 @@ +# Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration + # The name of the package: atlas_subdir( PMGAnalysisInterfaces ) @@ -5,6 +7,7 @@ atlas_subdir( PMGAnalysisInterfaces ) atlas_depends_on_subdirs( PUBLIC Control/AthToolSupport/AsgTools + Event/xAOD/xAODTruth PhysicsAnalysis/AnalysisCommon/PATInterfaces ) # Component(s) in the package: @@ -12,7 +15,7 @@ atlas_add_library( PMGAnalysisInterfacesLib PMGAnalysisInterfaces/*.h INTERFACE PUBLIC_HEADERS PMGAnalysisInterfaces - LINK_LIBRARIES AsgTools PATInterfaces ) + LINK_LIBRARIES AsgTools PATInterfaces xAODTruth ) atlas_add_dictionary( PMGAnalysisInterfacesDict PMGAnalysisInterfaces/PMGAnalysisInterfacesDict.h diff --git a/PhysicsAnalysis/Interfaces/PMGAnalysisInterfaces/PMGAnalysisInterfaces/ISysTruthWeightTool.h b/PhysicsAnalysis/Interfaces/PMGAnalysisInterfaces/PMGAnalysisInterfaces/ISysTruthWeightTool.h new file mode 100644 index 0000000000000000000000000000000000000000..066c852dc32c1ed3e0f5312a2ff7aa895adc5a18 --- /dev/null +++ b/PhysicsAnalysis/Interfaces/PMGAnalysisInterfaces/PMGAnalysisInterfaces/ISysTruthWeightTool.h @@ -0,0 +1,36 @@ +/* + Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef PMGTOOLS_ISYSTRUTHWEIGHTTOOLTOOL_H +#define PMGTOOLS_ISYSTRUTHWEIGHTTOOLTOOL_H + +// STL include(s): +#include +#include + +// EDM include(s): +#include "PATInterfaces/ISystematicsTool.h" +#include "xAODTruth/TruthParticleContainer.h" + +namespace PMGTools +{ +/// Interface for xAOD Truth Weight Tool which returns +/// an event weight based on truth particle container +/// +/// @author Miha Muskinja +/// +class ISysTruthWeightTool : public virtual CP::ISystematicsTool +{ + /// Declare the interface that the class provides + ASG_TOOL_INTERFACE(xAOD::ISysTruthWeightTool) + +public: + /// Return the weight corresponding to this event + virtual float getWeight(const xAOD::TruthParticleContainer *truthParticles) const = 0; + +}; // class ISysTruthWeightTool + +} // namespace PMGTools + +#endif // PMGTOOLS_ISYSTRUTHWEIGHTTOOLTOOL_H diff --git a/PhysicsAnalysis/Interfaces/PMGAnalysisInterfaces/PMGAnalysisInterfaces/PMGAnalysisInterfacesDict.h b/PhysicsAnalysis/Interfaces/PMGAnalysisInterfaces/PMGAnalysisInterfaces/PMGAnalysisInterfacesDict.h index d27b286c7a2c7fdfc71ffd851f21af28f40c9302..fc33b8c9a3f74bd9035e39c0c626b84e58b4e68e 100644 --- a/PhysicsAnalysis/Interfaces/PMGAnalysisInterfaces/PMGAnalysisInterfaces/PMGAnalysisInterfacesDict.h +++ b/PhysicsAnalysis/Interfaces/PMGAnalysisInterfaces/PMGAnalysisInterfaces/PMGAnalysisInterfacesDict.h @@ -1,6 +1,6 @@ // Dear emacs, this is -*- c++ -* // -// Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration +// Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration // #ifndef PMGANALYSISINTERFACES_PMGANALYSISINTERFACESDICT_H #define PMGANALYSISINTERFACES_PMGANALYSISINTERFACESDICT_H @@ -10,5 +10,6 @@ #include "PMGAnalysisInterfaces/IPMGCrossSectionTool.h" #include "PMGAnalysisInterfaces/IPMGSherpaVjetsSysTool.h" #include "PMGAnalysisInterfaces/IPMGTruthWeightTool.h" +#include "PMGAnalysisInterfaces/ISysTruthWeightTool.h" #endif // PMGANALYSISINTERFACES_PMGANALYSISINTERFACESDICT_H diff --git a/PhysicsAnalysis/Interfaces/PMGAnalysisInterfaces/PMGAnalysisInterfaces/selection.xml b/PhysicsAnalysis/Interfaces/PMGAnalysisInterfaces/PMGAnalysisInterfaces/selection.xml index 80a21d691ff8f03bd2bae6e054a122f9d8bafb4a..7e750f82586e3489a02ad4821c0b29bb29063f21 100644 --- a/PhysicsAnalysis/Interfaces/PMGAnalysisInterfaces/PMGAnalysisInterfaces/selection.xml +++ b/PhysicsAnalysis/Interfaces/PMGAnalysisInterfaces/PMGAnalysisInterfaces/selection.xml @@ -5,4 +5,5 @@ +