diff --git a/PhysicsAnalysis/AnalysisCommon/PMGTools/CMakeLists.txt b/PhysicsAnalysis/AnalysisCommon/PMGTools/CMakeLists.txt
index 6fa17a167933302bfb7c23e2ecbff2d9a07f30b0..389be61f89bffe56d007ea7a0f760ce04a17c3c0 100644
--- a/PhysicsAnalysis/AnalysisCommon/PMGTools/CMakeLists.txt
+++ b/PhysicsAnalysis/AnalysisCommon/PMGTools/CMakeLists.txt
@@ -1,4 +1,3 @@
-# $Id: CMakeLists.txt 782436 2016-11-04 15:35:58Z krasznaa $
#################################################################################
# Package: PMGTools
#################################################################################
@@ -7,38 +6,56 @@
atlas_subdir( PMGTools )
# Extra dependencies based on the build environment:
-set( extra_deps )
+set( extra_private_deps )
+# ... for AnalysisBase
if( XAOD_STANDALONE )
- set( extra_deps Control/xAODRootAccess )
+ set( extra_private_deps Control/xAODRootAccess )
+# ... for AthAnalysisBase/Athena
else()
- set( extra_deps Control/AthAnalysisBaseComps PhysicsAnalysis/POOLRootAccess
- GaudiKernel )
+ set( extra_private_deps Control/AthAnalysisBaseComps PhysicsAnalysis/POOLRootAccess GaudiKernel Event/EventInfo)
endif()
+# Extra libraries based on the build environment:
+set( xaod_access_lib )
+set( extra_private_libs )
+# ... for AnalysisBase
+if( XAOD_STANDALONE )
+ set( xaod_access_lib xAODRootAccess )
+# ... for AthAnalysisBase (Athena calls this POOLRootAccess)
+else()
+ set( xaod_access_lib POOLRootAccessLib )
+ set( extra_private_libs EventInfo )
+endif()
+
+
# Declare the package's dependencies:
atlas_depends_on_subdirs(
PUBLIC
Control/AthToolSupport/AsgTools
+ Event/xAOD/xAODTruth
+ PhysicsAnalysis/AnalysisCommon/PATCore
PhysicsAnalysis/AnalysisCommon/PATInterfaces
+ PhysicsAnalysis/Interfaces/PMGAnalysisInterfaces
PRIVATE
Event/FourMomUtils
Event/xAOD/xAODEventInfo
Event/xAOD/xAODJet
- Event/xAOD/xAODTruth
+ Event/xAOD/xAODMetaData
+ PhysicsAnalysis/D3PDTools/RootCoreUtils
Tools/PathResolver
- ${extra_deps} )
+ ${extra_private_deps} )
# External(s) used by the package:
find_package( ROOT COMPONENTS Core Hist RIO )
+find_package( Boost )
# Libraries in the package:
atlas_add_library( PMGToolsLib
PMGTools/*.h Root/*.cxx
PUBLIC_HEADERS PMGTools
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
- LINK_LIBRARIES ${ROOT_LIBRARIES} AsgTools PATInterfaces
- PRIVATE_LINK_LIBRARIES FourMomUtils xAODEventInfo xAODJet xAODTruth
- PathResolver )
+ LINK_LIBRARIES ${ROOT_LIBRARIES} AsgTools PATCoreLib PATInterfaces xAODTruth
+ PRIVATE_LINK_LIBRARIES FourMomUtils PathResolver RootCoreUtils xAODEventInfo xAODJet xAODMetaData ${extra_private_libs})
if( NOT XAOD_STANDALONE )
atlas_add_component( PMGTools
@@ -52,28 +69,36 @@ atlas_add_dictionary( PMGToolsDict
LINK_LIBRARIES PMGToolsLib )
# Executable(s) in the package:
+#
+# ... for Athena/AthAnalysis
if( NOT XAOD_STANDALONE )
atlas_add_executable( MyPMGApp
test/MyPMGApp.cxx
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
- LINK_LIBRARIES ${ROOT_LIBRARIES} AthAnalysisBaseCompsLib POOLRootAccess
- AsgTools xAODJet xAODTruth FourMomUtils PATInterfaces )
+ LINK_LIBRARIES ${ROOT_LIBRARIES} AthAnalysisBaseCompsLib ${xaod_access_lib}
+ AsgTools xAODJet xAODTruth FourMomUtils PATInterfaces PMGToolsLib )
endif()
# Test(s) in the package:
+#
+# ... for AnalysisBase
if( XAOD_STANDALONE )
atlas_add_test( ut_PMGSherpa22VJetsWeightTool_test
SOURCES test/ut_PMGSherpa22VJetsWeightTool_test.cxx
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
- LINK_LIBRARIES ${ROOT_LIBRARIES} AsgTools xAODRootAccess PATInterfaces
- PMGToolsLib )
+ LINK_LIBRARIES ${ROOT_LIBRARIES} AsgTools PATInterfaces PMGToolsLib ${xaod_access_lib} )
atlas_add_test( ut_PMGSherpaVjetsSysTool_test
SOURCES test/ut_PMGSherpaVjetsSysTool_test.cxx
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
- LINK_LIBRARIES ${ROOT_LIBRARIES} AsgTools xAODRootAccess PATInterfaces
- PMGToolsLib )
+ LINK_LIBRARIES ${ROOT_LIBRARIES} AsgTools PATInterfaces PMGToolsLib ${xaod_access_lib} )
endif()
+# ... AthAnalysis/AnalysisBase
+atlas_add_test( ut_PMGTruthWeightTool_test
+ SOURCES test/ut_PMGTruthWeightTool_test.cxx
+ INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
+ LINK_LIBRARIES ${ROOT_LIBRARIES} PMGToolsLib ${xaod_access_lib} )
+
# Install files from the package:
atlas_install_data( data/*.txt share/*.txt )
diff --git a/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/PMGCrossSectionTool.h b/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/PMGCrossSectionTool.h
index e475076a9963fe04a9aa65ba6b934f35324dc008..f9af041e40f19e31b6108bcca5e3656dddde263e 100644
--- a/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/PMGCrossSectionTool.h
+++ b/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/PMGCrossSectionTool.h
@@ -8,13 +8,14 @@
#define PMGTOOLS_PMGCROSSSECTIONTOOL_H
// System include(s):
-#include <array>
+#include <map>
+#include <string>
// Infrastructure include(s):
#include "AsgTools/AsgTool.h"
// Interface include(s):
-#include "PMGTools/IPMGCrossSectionTool.h"
+#include "PMGAnalysisInterfaces/IPMGCrossSectionTool.h"
/// Tool providing sample cross-sections and k-factors etc
///
@@ -51,9 +52,21 @@ namespace PMGTools {
std::string getSampleName(const int dsid) const;
/// return the AMI cross-section for DSID
- double getGeneratorXsection(const int dsid) const;
-
+ double getAMIXsection(const int dsid) const;
+
+ /// return the cross-section uncertainty for DSID
+ double getXsectionUncertainty(const int dsid) const;
+
+ /// return the cross-section uncertainty for DSID
+ double getXsectionUncertaintyUP(const int dsid) const;
+
+ /// return the cross-section uncertainty for DSID
+ double getXsectionUncertaintyDOWN(const int dsid) const;
+ // :: below is for future use?
+ /// return the branching ratio for DSID
+ //double getBR(const int dsid) const;
+
/// return the k-factor for DSID
double getKfactor(const int dsid) const;
@@ -66,7 +79,7 @@ namespace PMGTools {
private:
// store vector of structures, each structure contains full info for DSID
- std::map<int,PMGTools::AllSampleInfo> m_storeSampleInfo;
+ std::map<unsigned, PMGTools::AllSampleInfo> m_fStoreSampleInfo;
std::string m_InputFileName;
}; // class PMGCrossSectionTool
diff --git a/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/PMGDecayProductsSelectionTool.h b/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/PMGDecayProductsSelectionTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..13ada8706d38281c9895c011299fcb33fdfe64b3
--- /dev/null
+++ b/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/PMGDecayProductsSelectionTool.h
@@ -0,0 +1,91 @@
+/*
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+ */
+
+/// @author Tadej Novak
+
+
+
+#ifndef PMG_TOOLS__PMG_DECAY_PRODUCTS_SELECTION_TOOL_H
+#define PMG_TOOLS__PMG_DECAY_PRODUCTS_SELECTION_TOOL_H
+
+#include <AsgTools/AsgTool.h>
+#include <PATCore/IAsgSelectionTool.h>
+#include <xAODTruth/TruthParticle.h>
+#include <xAODTruth/TruthParticleContainer.h>
+
+
+namespace PMGTools
+{
+ /// \brief an \ref IAsgSelectionTool that select particles
+ /// based on the allowed decay chain
+
+ class PMGDecayProductsSelectionTool final
+ : public asg::AsgTool, virtual public IAsgSelectionTool
+ {
+ // Create a proper constructor for Athena
+ ASG_TOOL_CLASS (PMGDecayProductsSelectionTool, IAsgSelectionTool)
+
+
+ /// \brief standard constructor
+ public:
+ PMGDecayProductsSelectionTool (const std::string& name);
+
+
+
+
+ //
+ // inherited interface
+ //
+
+ virtual StatusCode initialize () override;
+
+ virtual const asg::AcceptInfo& getAcceptInfo () const override;
+
+ virtual asg::AcceptData accept (const xAOD::IParticle* /*part*/) const override;
+
+
+
+ //
+ // private interface
+ //
+
+ /// \brief Helper function to check for required parent particles
+ private:
+ asg::AcceptData hasRequiredInitialParent (const xAOD::TruthParticle *truthParticle, asg::AcceptData& acceptData) const;
+
+ /// \brief Helper function to get the number of parent particles
+ private:
+ size_t getNParents (const xAOD::TruthParticle *truthParticle) const;
+
+ /// \brief Helper function get a parent by index
+ private:
+ const xAOD::TruthParticle* getParent (const xAOD::TruthParticle *truthParticle,
+ size_t index) const;
+
+ /// tool properties
+ /// \{
+
+ private:
+ std::vector<int> m_requiredParentPDGIDs;
+ std::vector<int> m_allowedIntermediatePDGIDs;
+
+ /// \}
+
+ private:
+ /// \brief Index for the truth particle link
+ int m_truthParticleIndex{ -1 };
+ /// \brief Index for the required parent particles
+ int m_requiredParentIndex{ -1 };
+
+ /// \brief the \ref AcceptInfo we are using
+ private:
+ asg::AcceptInfo m_accept;
+
+ /// \brief common parents accessor
+ private:
+ std::unique_ptr<const SG::AuxElement::Accessor<std::vector<ElementLink<xAOD::TruthParticleContainer>>>> m_parentsAccessor{};
+ };
+}
+
+#endif
diff --git a/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/PMGSherpaVjetsSysTool.h b/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/PMGSherpaVjetsSysTool.h
index 086e21e0b3e149f659b002df8363d0b723f65c1d..7bef3eb60c11125e4b60eae8dd1466bb3a97ff98 100644
--- a/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/PMGSherpaVjetsSysTool.h
+++ b/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/PMGSherpaVjetsSysTool.h
@@ -21,7 +21,7 @@
#include "AsgTools/AsgTool.h"
// Interface include(s):
-#include "PMGTools/IPMGSherpaVjetsSysTool.h"
+#include "PMGAnalysisInterfaces/IPMGSherpaVjetsSysTool.h"
// System include(s):
#include <map>
@@ -29,8 +29,10 @@
#include <memory>
// ROOT include(s):
-#include <TH2F.h>
-#include <TFile.h>
+#include "TFile.h"
+
+// Forward declaration
+class TH2F;
namespace PMGTools {
diff --git a/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/PMGToolsDict.h b/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/PMGToolsDict.h
index 2c833517ea242af8968d435ee1e956fc4b823df0..5066d2443fe33c4049bf2691302086a69439887f 100644
--- a/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/PMGToolsDict.h
+++ b/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/PMGToolsDict.h
@@ -9,7 +9,9 @@
#define PMGTOOLS_PMGTOOLSDICT_H
// Local include(s):
-#include "PMGTools/PMGSherpa22VJetsWeightTool.h"
#include "PMGTools/PMGCrossSectionTool.h"
+#include "PMGTools/PMGDecayProductsSelectionTool.h"
+#include "PMGTools/PMGSherpa22VJetsWeightTool.h"
+#include "PMGTools/PMGTruthWeightTool.h"
#endif // PMGTOOLS_PMGTOOLSDICT_H
diff --git a/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/PMGTruthWeightTool.h b/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/PMGTruthWeightTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..935571ad6e8edc2f9b9f02d831c044dd8d4a6bf3
--- /dev/null
+++ b/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/PMGTruthWeightTool.h
@@ -0,0 +1,150 @@
+/*
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef PMGTOOLS_PMGTRUTHWEIGHTTOOL_H
+#define PMGTOOLS_PMGTRUTHWEIGHTTOOL_H
+
+// EDM include(s):
+#include "AsgTools/AsgMetadataTool.h"
+#include "xAODTruth/TruthEventContainer.h"
+#include "xAODTruth/TruthMetaData.h"
+#include "xAODTruth/TruthMetaDataContainer.h"
+
+// Interface include(s):
+#include "PMGAnalysisInterfaces/IPMGTruthWeightTool.h"
+
+
+namespace PMGTools
+{
+ /// \brief the systematics prefix for generator weights
+ constexpr char generatorSystematicsPrefix[] {"GEN_"};
+
+ /// Implementation for the xAOD truth meta data weight tool
+ ///
+ /// @author James Robinson <james.robinson@cern.ch>
+ ///
+ class PMGTruthWeightTool : public virtual IPMGTruthWeightTool, public asg::AsgMetadataTool
+ {
+ /// Create a proper constructor for Athena
+ ASG_TOOL_CLASS(PMGTruthWeightTool, PMGTools::IPMGTruthWeightTool)
+
+ public:
+ /// Create a constructor for standalone usage
+ PMGTruthWeightTool(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 IPMGTruthWeightTool interface
+ /// @{
+
+ /// Implements interface from IPMGTruthWeightTool
+ virtual const std::vector<std::string>& getWeightNames() const override;
+
+ /// Implements interface from IPMGTruthWeightTool
+ virtual float getWeight(const std::string& weightName) const override;
+
+ /// Implements interface from IPMGTruthWeightTool
+ virtual bool hasWeight(const std::string& weightName) const override;
+
+ /// Implements interface from IPMGTruthWeightTool
+ virtual float getSysWeight() const override;
+
+ /// Implements interface from IPMGTruthWeightTool
+ virtual size_t getSysWeightIndex() const override;
+
+ /// @}
+
+ /// @name Function(s) implementing the ISystematicsTool interface
+ /// @{
+
+ /// Implements interface from ISystematicsTool
+ virtual bool isAffectedBySystematic(const CP::SystematicVariation& systematic) const override;
+
+ /// Implements interface from ISystematicsTool
+ virtual CP::SystematicSet affectingSystematics() const override;
+
+ /// Implements interface from ISystematicsTool
+ virtual CP::SystematicSet recommendedSystematics() const override;
+
+ /// Implements interface from ISystematicsTool
+ virtual CP::SystematicCode applySystematicVariation(const CP::SystematicSet& systConfig) override;
+
+ /// @}
+
+ protected:
+ /// @name Callback function(s) from AsgMetadataTool
+ /// @{
+
+ /// Function called when a new input file is opened
+ virtual StatusCode beginInputFile() override;
+
+ /// Function called when a new event is loaded
+ virtual StatusCode beginEvent() override;
+
+ /// @}
+
+ /// Loads weight information from xAOD::TruthMetaDataContainer
+ StatusCode loadMetaData();
+
+ /// Loads weight information from POOL using HepMCWeightNames
+ StatusCode loadPOOLMetaData();
+
+ /// Validate weight caches
+ StatusCode validateWeightLocationCaches();
+
+ /// Clear caches
+ void clearWeightLocationCaches();
+
+ /// Process the weight name to be used for systematics
+ std::string weightNameToSys(const std::string &name) const;
+
+ /// Stores the meta data record name
+ std::string m_metaName;
+
+ /// Current MC channel number
+ uint32_t m_mcChannelNumber{};
+
+ /// Systematics set of the weight systematics
+ CP::SystematicSet m_systematicsSet;
+
+ /// Ptr to the meta data container for this file
+ const xAOD::TruthMetaDataContainer* m_metaDataContainer{nullptr};
+
+ /// Flag to indicate whether the xAOD::TruthMetaData objects have incorrect McChannelNumber
+ bool m_useChannelZeroInMetaData{false};
+
+ /// Available weight names for this file
+ std::vector<std::string> m_weightNames;
+
+ /// Indices of available weights in this file
+ std::unordered_map<std::string, size_t> m_weightIndices;
+
+ /// Indices of available weights in this file
+ std::unordered_map<std::string, size_t> m_weightIndicesSys;
+
+ /// Values of weights in this event
+ std::vector<float> m_weights;
+
+ /// Weight data cache helper struct
+ struct WeightData
+ {
+ float weight;
+ std::size_t index;
+ };
+
+ /// Systematics cache of available weights in this file
+ std::unordered_map<CP::SystematicSet, WeightData> m_weightData;
+
+ /// Current systematics weight
+ WeightData *m_currentWeightData {nullptr};
+ };
+} // namespace PMGTools
+
+#endif // PMGTOOLS_PMGTRUTHWEIGHTTOOL_H
diff --git a/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/selection.xml b/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/selection.xml
index 624d061a6320169f0c03caa98151870bde4ebf64..7141c2b611a64d6824a1bae5f8cc9c336eb40daa 100644
--- a/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/selection.xml
+++ b/PhysicsAnalysis/AnalysisCommon/PMGTools/PMGTools/selection.xml
@@ -1,9 +1,10 @@
<!-- $Id: selection.xml 772287 2016-09-08 13:53:24Z tripiana $ -->
<lcgdict>
- <class name="PMGTools::PMGSherpa22VJetsWeightTool" />
<class name="PMGTools::PMGCrossSectionTool" />
+ <class name="PMGTools::PMGDecayProductsSelectionTool" />
+ <class name="PMGTools::PMGSherpa22VJetsWeightTool" />
<class name="PMGTools::PMGSherpaVjetsSysTool" />
- <class name="PMGTools::IPMGSherpaVjetsSysTool" />
+ <class name="PMGTools::PMGTruthWeightTool" />
</lcgdict>
diff --git a/PhysicsAnalysis/AnalysisCommon/PMGTools/Root/PMGCrossSectionTool.cxx b/PhysicsAnalysis/AnalysisCommon/PMGTools/Root/PMGCrossSectionTool.cxx
index 19d5c0a402e18b66db0dbf5abafa876b4dd08d3e..aaa346dbd3c2e52b2d1f20cf5b8d8e008a880aa5 100644
--- a/PhysicsAnalysis/AnalysisCommon/PMGTools/Root/PMGCrossSectionTool.cxx
+++ b/PhysicsAnalysis/AnalysisCommon/PMGTools/Root/PMGCrossSectionTool.cxx
@@ -2,195 +2,216 @@
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
-// $Id: PMGCrossSectionTool.cxx 810599 2017-09-22 15:50:07Z cohm $
+// $Id: PMGCrossSectionTool.cxx 764400 2016-07-26 17:47:39Z tripiana $
#include <fstream>
#include <sstream>
#include <string>
#include <TSystem.h>
#include <TString.h>
-#include <TTree.h>
-#include <iostream>
#include <stdlib.h>
+#include <limits>
// Local include(s):
#include "PMGTools/PMGCrossSectionTool.h"
namespace PMGTools {
-
- PMGCrossSectionTool::
- PMGCrossSectionTool( const std::string& name ) : asg::AsgTool( name )
- {
-
- }
-
- StatusCode PMGCrossSectionTool::initialize() {
-
- // Tell the user what's happening:
- ATH_MSG_INFO( "Initializing " << name() << "..." );
- //ATH_MSG_INFO( "Read in all Xsec Info ..." );
- //readInfosFromDir(inputDir.c_str());
+PMGCrossSectionTool::
+PMGCrossSectionTool(const std::string& name) : asg::AsgTool(name) { }
- // Return gracefully:
- return StatusCode::SUCCESS;
- }
- bool PMGCrossSectionTool::readInfosFromFiles(std::vector<std::string> InputFiles)
- {
- for (const auto& currentFileName : InputFiles) {
-
- std::ifstream currentFile;
- currentFile.open(currentFileName.c_str());
-
- // use the TTree functionality for reading in text files - easy way to remove
- // dependency on specific order of the columns in the text file
- TTree* xsecTree = new TTree();
- xsecTree->ReadFile(currentFileName.c_str());
-
- // one branch per variable we care about of the existing ones: DSID/I:DSName/C:GenXsec/D:BR:FilterEff:Xsec:Kfactor:TotSampleXsec
- int dsid(0);
- double genXsec(0);
- double filterEff(0);
- double kFactor(0);
- TBranch* dsidBranch = xsecTree->GetBranch("DSID");
- dsidBranch->SetAddress(&dsid);
- TBranch* genXsecBranch = xsecTree->GetBranch("crossSection");
- genXsecBranch->SetAddress(&genXsec);
- TBranch* filterEffBranch = xsecTree->GetBranch("genFiltEff");
- filterEffBranch->SetAddress(&filterEff);
- TBranch* kFactorBranch = xsecTree->GetBranch("kFactor");
- kFactorBranch->SetAddress(&kFactor);
-
- auto nEntries = xsecTree->GetEntries();
- for (Long64_t i = 0; i < nEntries; i++) {
- xsecTree->GetEntry(i);
- AllSampleInfo sample;
- sample.dsid = dsid;
- sample.genXsec = genXsec;
- sample.filterEff = filterEff;
- sample.kFactor = kFactor;
- //sample.containerName = containerName;
- m_storeSampleInfo[sample.dsid] = sample;
- }
-
- delete xsecTree;
- }
-
- return true;
- }
+StatusCode PMGCrossSectionTool::initialize() {
- bool PMGCrossSectionTool::readInfosFromDir(const std::string& inputDir)
- {
+ // Tell the user what's happening:
+ ATH_MSG_INFO( "Initializing " << name() << "..." );
+ //ATH_MSG_INFO( "Read in all Xsec Info ..." );
- TString mydir = inputDir;
- gSystem->ExpandPathName (mydir);
- void *dirp = 0;
-
- std::vector<std::string> inFiles={};
-
- try
- {
- dirp = gSystem->OpenDirectory (mydir.Data());
- const char *file = 0;
- while ((file = gSystem->GetDirEntry (dirp)))
- {
- std::string myfile = inputDir + "/" + file;
- if (myfile.size() > 4 && myfile.substr (myfile.size()-4) == ".txt")
- inFiles.push_back(myfile);
- }
- gSystem->FreeDirectory (dirp);
- }
- catch (...)
- {
- gSystem->FreeDirectory (dirp);
- //throw;
- return false;
- }
-
- readInfosFromFiles(inFiles);
- return true;
- }
+ //readInfosFromDir(inputDir.c_str());
+ // Return gracefully:
+ return StatusCode::SUCCESS;
+}
- double PMGCrossSectionTool::getFilterEff(const int dsid) const
- {
-
- for (const auto& info : m_storeSampleInfo){
- if(dsid == info.second.dsid)
- return info.second.filterEff;
+
+bool PMGCrossSectionTool::readInfosFromFiles(std::vector<std::string> InputFiles)
+{
+
+ for (const auto& currentFileName : InputFiles) {
+
+ std::ifstream currentFile(currentFileName);
+ if (not currentFile.is_open()) {
+ ATH_MSG_WARNING("cannot open file " << currentFileName);
+ continue;
}
-
- std::cout << "ERROR::getFilterEff --> Sample with DSID " << dsid << " has no info stored!!! ---> EXIT." << std::endl;
-
- return -1;
-
+
+ // skip first line
+ currentFile.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
+
+ int nfound = 0; // for this file
+ for (std::string line; std::getline(currentFile, line); ) {
+ std::stringstream input_line(line);
+ AllSampleInfo help;
+
+ input_line >> help.dsid;
+ input_line >> help.containerName;
+ input_line >> help.amiXSec;
+ input_line >> help.filterEff;
+ input_line >> help.kFactor;
+ input_line >> help.XSecUncUP;
+ input_line >> help.XSecUncDOWN;
+ // :: below is for future use?
+ //input_line >> help.br;
+ //input_line >> help.higherOrderXsecTotal;
+ //input_line >> help.higherOrderXsecSample; // this is hoXsec * filter eff
+
+ if (input_line.fail()) { ATH_MSG_ERROR("cannot parse line '" << line << "' from file " << currentFileName); continue; }
+
+ m_fStoreSampleInfo[help.dsid] = help;
+ ++nfound;
+ }
+
+ if (nfound == 0) { ATH_MSG_WARNING("no sample read from file " << currentFileName); }
}
+ if (m_fStoreSampleInfo.empty()) {
+ ATH_MSG_ERROR("list of sample is empty");
+ return false;
+ }
- std::string PMGCrossSectionTool::getSampleName(const int dsid) const
+ return true;
+}
+
+
+bool PMGCrossSectionTool::readInfosFromDir(const std::string& inputDir)
+{
+
+ TString mydir = inputDir;
+ gSystem->ExpandPathName (mydir);
+ void *dirp = 0;
+
+ std::vector<std::string> inFiles = {};
+
+ try
{
-
- for (const auto& info : m_storeSampleInfo){
- if(dsid == info.second.dsid)
- return info.second.containerName;
+ dirp = gSystem->OpenDirectory (mydir.Data());
+ const char *file = 0;
+ while ((file = gSystem->GetDirEntry (dirp)))
+ {
+ std::string myfile = inputDir + "/" + file;
+ if (myfile.size() > 4 && myfile.substr (myfile.size() - 4) == ".txt")
+ inFiles.push_back(myfile);
}
-
- std::cout << "ERROR::getSampleName --> Sample with DSID " << dsid << " has no info stored!!! ---> EXIT." << std::endl;
-
- return "";
-
+ gSystem->FreeDirectory (dirp);
+ }
+ catch (...)
+ {
+ gSystem->FreeDirectory (dirp);
+ //throw;
+ return false;
}
+ readInfosFromFiles(inFiles);
+ return true;
+}
- double PMGCrossSectionTool::getGeneratorXsection(const int dsid) const
- {
- for (const auto& info : m_storeSampleInfo){
- if(dsid == info.second.dsid)
- return info.second.genXsec;
- }
+double PMGCrossSectionTool::getFilterEff(const int dsid) const
+{
+ const auto it = m_fStoreSampleInfo.find(dsid);
+ if (it != m_fStoreSampleInfo.end()) { return it->second.filterEff; }
- std::cout << "ERROR::getGeneratorXsection --> Sample with DSID " << dsid << " has no info stored!!! ---> EXIT." << std::endl;
-
- return -1;
+ ATH_MSG_ERROR("Sample with DSID " << dsid << " has no info stored!!!");
+ return -1;
+}
- }
- double PMGCrossSectionTool::getKfactor(const int dsid) const
- {
+std::string PMGCrossSectionTool::getSampleName(const int dsid) const
+{
+ const auto it = m_fStoreSampleInfo.find(dsid);
+ if (it != m_fStoreSampleInfo.end()) { return it->second.containerName; }
- for (const auto& info : m_storeSampleInfo){
- if(dsid == info.second.dsid)
- return info.second.kFactor;
- }
+ ATH_MSG_ERROR("Sample with DSID " << dsid << " has no info stored!!!");
+ return std::string("");
+}
- std::cout << "ERROR::getKfactor --> Sample with DSID " << dsid << " has no info stored!!! ---> EXIT." << std::endl;
- return -1;
+double PMGCrossSectionTool::getAMIXsection(const int dsid) const
+{
+ const auto it = m_fStoreSampleInfo.find(dsid);
+ if (it != m_fStoreSampleInfo.end()) { return it->second.amiXSec; }
- }
+ ATH_MSG_ERROR("Sample with DSID " << dsid << " has no info stored!!!");
+ return -1;
+}
- double PMGCrossSectionTool::getSampleXsection(const int dsid) const
- {
- for (const auto& info : m_storeSampleInfo){
- if(dsid == info.second.dsid)
- return info.second.genXsec * info.second.kFactor * info.second.filterEff;
- }
+// :: below is for future use?
+/*double PMGCrossSectionTool::getBR(const int dsid) const
+{
+ const auto it = m_fStoreSampleInfo.find(dsid);
+ if (it != m_fStoreSampleInfo.end()) { return it->second.br; }
- std::cout << "ERROR::getSampleXsection --> Sample with DSID " << dsid << " has no info stored!!! ---> EXIT." << std::endl;
+ ATH_MSG_ERROR("Sample with DSID " << dsid << " has no info stored!!!");
+ return -1;
- return -1;
+}*/
- }
- std::vector<int> PMGCrossSectionTool::getLoadedDSIDs() const {
- std::vector<int> dsids;
- for (const auto& info : m_storeSampleInfo){
- dsids.push_back(info.second.dsid);
- }
- return dsids;
- }
+double PMGCrossSectionTool::getXsectionUncertaintyUP(const int dsid) const
+{
+ const auto it = m_fStoreSampleInfo.find(dsid);
+ if (it != m_fStoreSampleInfo.end()) { return it->second.XSecUncUP; }
+ ATH_MSG_ERROR("Sample with DSID " << dsid << " has no info stored!!!");
+ return -1;
}
+
+double PMGCrossSectionTool::getXsectionUncertaintyDOWN(const int dsid) const
+{
+ const auto it = m_fStoreSampleInfo.find(dsid);
+ if (it != m_fStoreSampleInfo.end()) { return it->second.XSecUncDOWN; }
+
+ ATH_MSG_ERROR("Sample with DSID " << dsid << " has no info stored!!!");
+ return -1;
+}
+
+double PMGCrossSectionTool::getXsectionUncertainty(const int dsid) const
+{
+ //symmetrize the up and down variations
+ return 0.5*( fabs(getXsectionUncertaintyDOWN(dsid)) + fabs(getXsectionUncertaintyUP(dsid)) );
+}
+
+double PMGCrossSectionTool::getKfactor(const int dsid) const
+{
+ const auto it = m_fStoreSampleInfo.find(dsid);
+ if (it != m_fStoreSampleInfo.end()) { return it->second.kFactor; }
+
+ ATH_MSG_ERROR("Sample with DSID " << dsid << " has no info stored!!!");
+ return -1;
+}
+
+
+// :: below is for future use?
+double PMGCrossSectionTool::getSampleXsection(const int dsid) const
+{
+ const auto it = m_fStoreSampleInfo.find(dsid);
+ if (it != m_fStoreSampleInfo.end()) { return it->second.amiXSec * it->second.kFactor * it->second.filterEff; }
+ // :: below is for future use?
+ //return info.higherOrderXsecSample;
+
+ ATH_MSG_ERROR("Sample with DSID " << dsid << " has no info stored!!!");
+ return -1;
+}
+
+
+std::vector<int> PMGCrossSectionTool::getLoadedDSIDs() const {
+ std::vector<int> dsids;
+ dsids.reserve(m_fStoreSampleInfo.size());
+ for (const std::pair<unsigned, AllSampleInfo>& key_info : m_fStoreSampleInfo) {
+ dsids.push_back(key_info.second.dsid);
+ }
+ return dsids;
+}
+
+} // end namespace
diff --git a/PhysicsAnalysis/AnalysisCommon/PMGTools/Root/PMGDecayProductsSelectionTool.cxx b/PhysicsAnalysis/AnalysisCommon/PMGTools/Root/PMGDecayProductsSelectionTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..b737b0446db6c5a4d76685278c0a465930273ecb
--- /dev/null
+++ b/PhysicsAnalysis/AnalysisCommon/PMGTools/Root/PMGDecayProductsSelectionTool.cxx
@@ -0,0 +1,177 @@
+/*
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+/// @author Tadej Novak
+
+
+
+//
+// includes
+//
+
+#include <PMGTools/PMGDecayProductsSelectionTool.h>
+
+#include <xAODTruth/xAODTruthHelpers.h>
+#include <xAODBase/IParticle.h>
+
+#include <sstream>
+
+//
+// method implementations
+//
+
+namespace PMGTools
+{
+ PMGDecayProductsSelectionTool ::
+ PMGDecayProductsSelectionTool (const std::string& name)
+ : AsgTool (name)
+ {
+ declareProperty ("requiredParentPDGIDs", m_requiredParentPDGIDs, "required parent particle PDG IDs (positive only)");
+ declareProperty ("allowedIntermediatePDGIDs", m_allowedIntermediatePDGIDs, "allowed intermediate particles PDG IDs (positive only)");
+ }
+
+
+
+ StatusCode PMGDecayProductsSelectionTool ::
+ initialize ()
+ {
+ if (m_requiredParentPDGIDs.empty())
+ {
+ ATH_MSG_ERROR ("required parent particles need to be set");
+ return StatusCode::FAILURE;
+ }
+
+ m_truthParticleIndex = m_accept.addCut ("truthParticle", "has truth particle cut");
+
+ std::ostringstream particles;
+ for (int pdgId : m_requiredParentPDGIDs)
+ {
+ particles << " " << pdgId;
+ }
+ ATH_MSG_DEBUG ("Performing required parent particle selection with parent PDG IDs:" << particles.str());
+ m_requiredParentIndex = m_accept.addCut ("requiredParents", "required parent particles cut");
+
+ if (!m_allowedIntermediatePDGIDs.empty())
+ {
+ std::ostringstream particles;
+ for (int pdgId : m_allowedIntermediatePDGIDs)
+ {
+ particles << " " << pdgId;
+ }
+ ATH_MSG_DEBUG ("Performing allowed intermediate particle selection with allowed PDG IDs:" << particles.str());
+ }
+
+ m_parentsAccessor = std::make_unique<const SG::AuxElement::Accessor<std::vector<ElementLink<xAOD::TruthParticleContainer>>>>("parentLinks");
+
+ return StatusCode::SUCCESS;
+ }
+
+
+
+ const asg::AcceptInfo& PMGDecayProductsSelectionTool ::
+ getAcceptInfo () const
+ {
+ return m_accept;
+ }
+
+
+
+ asg::AcceptData PMGDecayProductsSelectionTool ::
+ accept (const xAOD::IParticle *particle) const
+ {
+ asg::AcceptData acceptData (&m_accept);
+
+ // Check if xAOD::TruthParticle or if not if it has the TruthParticleLink
+ const xAOD::TruthParticle *truthParticle
+ = dynamic_cast<const xAOD::TruthParticle *> (particle);
+ if (truthParticle == nullptr)
+ {
+ // need to find the truth particle
+ truthParticle = xAOD::TruthHelpers::getTruthParticle(*particle);
+ if (truthParticle == nullptr)
+ {
+ acceptData.setCutResult (m_truthParticleIndex, false);
+ return acceptData;
+ }
+ }
+ acceptData.setCutResult (m_truthParticleIndex, true);
+
+ return hasRequiredInitialParent(truthParticle, acceptData);
+ }
+
+
+
+ asg::AcceptData PMGDecayProductsSelectionTool ::
+ hasRequiredInitialParent (const xAOD::TruthParticle *truthParticle, asg::AcceptData& acceptData) const
+ {
+ size_t nParents = getNParents (truthParticle);
+ for (size_t i = 0; i < nParents; i++)
+ {
+ const xAOD::TruthParticle *parent = getParent(truthParticle, i);
+ if (parent)
+ {
+ if (std::find(m_requiredParentPDGIDs.begin(), m_requiredParentPDGIDs.end(), std::abs(parent->pdgId())) != m_requiredParentPDGIDs.end())
+ {
+ acceptData.setCutResult (m_requiredParentIndex, true);
+ return acceptData;
+ }
+ else if (m_allowedIntermediatePDGIDs.empty() || std::find(m_allowedIntermediatePDGIDs.begin(), m_allowedIntermediatePDGIDs.end(), std::abs(parent->pdgId())) != m_allowedIntermediatePDGIDs.end())
+ {
+ return hasRequiredInitialParent(parent, acceptData);
+ }
+ else
+ {
+ ATH_MSG_VERBOSE("Removing particle as parent is not allowed: " << parent->pdgId());
+ return acceptData;
+ }
+ }
+ else
+ {
+ ATH_MSG_WARNING("Particle parent is not valid");
+ return acceptData;
+ }
+ }
+
+ if (std::find(m_requiredParentPDGIDs.begin(), m_requiredParentPDGIDs.end(), std::abs(truthParticle->pdgId())) != m_requiredParentPDGIDs.end())
+ {
+ acceptData.setCutResult (m_requiredParentIndex, true);
+ return acceptData;
+ }
+
+ return acceptData;
+ }
+
+
+
+ size_t PMGDecayProductsSelectionTool ::
+ getNParents (const xAOD::TruthParticle *truthParticle) const
+ {
+ if (m_parentsAccessor->isAvailable (*truthParticle))
+ {
+ return (*m_parentsAccessor)(*truthParticle).size();
+ }
+ else
+ {
+ return truthParticle->nParents();
+ }
+ }
+
+
+
+ const xAOD::TruthParticle* PMGDecayProductsSelectionTool ::
+ getParent (const xAOD::TruthParticle *truthParticle,
+ size_t index) const
+ {
+ if (m_parentsAccessor->isAvailable (*truthParticle))
+ {
+ ElementLink<xAOD::TruthParticleContainer> parentElementLink
+ = (*m_parentsAccessor)(*truthParticle).at(index);
+ return parentElementLink.isValid() ? *parentElementLink : nullptr;
+ }
+ else
+ {
+ return truthParticle->parent(index);
+ }
+ }
+}
diff --git a/PhysicsAnalysis/AnalysisCommon/PMGTools/Root/PMGSherpaVjetsSysTool.cxx b/PhysicsAnalysis/AnalysisCommon/PMGTools/Root/PMGSherpaVjetsSysTool.cxx
index 08e3c65b0995496dd14dac4bdae6b329d9a59009..73781d56a9401ad4cca43c08a429cd1de9a09a9a 100644
--- a/PhysicsAnalysis/AnalysisCommon/PMGTools/Root/PMGSherpaVjetsSysTool.cxx
+++ b/PhysicsAnalysis/AnalysisCommon/PMGTools/Root/PMGSherpaVjetsSysTool.cxx
@@ -5,8 +5,9 @@
// $Id: PMGSherpaVjetsSysTool.cxx.cxx 2016-07-12 tripiana $
// ROOT include(s):
-#include <TKey.h>
-#include <TIterator.h>
+#include "TKey.h"
+#include "TIterator.h"
+#include "TH2F.h"
// EDM include(s):
#include "xAODEventInfo/EventInfo.h"
@@ -326,6 +327,64 @@ namespace PMGTools {
if( MCID>=361054 && MCID<361057 ){ pType=SysParType::GammaJets; return 6; }
if( MCID>=361057 && MCID<361061 ){ pType=SysParType::GammaJets; return 7; }
+
+ //Check Znunu 2.2.1 samples
+ if( MCID>=364142 && MCID<364145 ){ pType=SysParType::Znunu; return 1; }
+ if( MCID>=364145 && MCID<364148 ){ pType=SysParType::Znunu; return 2; }
+ if( MCID>=364148 && MCID<364151 ){ pType=SysParType::Znunu; return 3; }
+ if( MCID>=364151 && MCID<364154 ){ pType=SysParType::Znunu; return 4; }
+ if( MCID==364154 ){ pType=SysParType::Znunu; return 5; }
+ if( MCID==364155 ){ pType=SysParType::Znunu; return 7; }
+
+ //Check Zee 2.2.1 samples
+ if( MCID>=364114 && MCID<364117 ){ pType=SysParType::Zll; return 1; }
+ if( MCID>=364117 && MCID<364120 ){ pType=SysParType::Zll; return 2; }
+ if( MCID>=364120 && MCID<364123 ){ pType=SysParType::Zll; return 3; }
+ if( MCID>=364123 && MCID<364126 ){ pType=SysParType::Zll; return 4; }
+ if( MCID==364126 ){ pType=SysParType::Zll; return 5; }
+ if( MCID==364127 ){ pType=SysParType::Zll; return 7; }
+
+ //Check Zmumu 2.2.1 samples
+ if( MCID>=364100 && MCID<364103 ){ pType=SysParType::Zll; return 1; }
+ if( MCID>=364103 && MCID<364106 ){ pType=SysParType::Zll; return 2; }
+ if( MCID>=364106 && MCID<364109 ){ pType=SysParType::Zll; return 3; }
+ if( MCID>=364109 && MCID<364112 ){ pType=SysParType::Zll; return 4; }
+ if( MCID==364112 ){ pType=SysParType::Zll; return 5; }
+ if( MCID==364113 ){ pType=SysParType::Zll; return 7; }
+
+ //Check Ztautau 2.2.1 samples
+ if( MCID>=364128 && MCID<364131 ){ pType=SysParType::Zll; return 1; }
+ if( MCID>=364131 && MCID<364134 ){ pType=SysParType::Zll; return 2; }
+ if( MCID>=364134 && MCID<364137 ){ pType=SysParType::Zll; return 3; }
+ if( MCID>=364137 && MCID<364140 ){ pType=SysParType::Zll; return 4; }
+ if( MCID==364140 ){ pType=SysParType::Zll; return 5; }
+ if( MCID==364141 ){ pType=SysParType::Zll; return 7; }
+
+ //Check Wenu 2.2.1 samples
+ if( MCID>=364170 && MCID<364173 ){ pType=SysParType::Wlnu; return 1; }
+ if( MCID>=364173 && MCID<364176 ){ pType=SysParType::Wlnu; return 2; }
+ if( MCID>=364176 && MCID<364179 ){ pType=SysParType::Wlnu; return 3; }
+ if( MCID>=364179 && MCID<364182 ){ pType=SysParType::Wlnu; return 4; }
+ if( MCID==364182 ){ pType=SysParType::Wlnu; return 5; }
+ if( MCID==364183 ){ pType=SysParType::Wlnu; return 7; }
+
+ //Check Wmunu 2.2.1 samples
+ if( MCID>=364156 && MCID<364159 ){ pType=SysParType::Wlnu; return 1; }
+ if( MCID>=364159 && MCID<364162 ){ pType=SysParType::Wlnu; return 2; }
+ if( MCID>=364162 && MCID<364165 ){ pType=SysParType::Wlnu; return 3; }
+ if( MCID>=364165 && MCID<364168 ){ pType=SysParType::Wlnu; return 4; }
+ if( MCID==364168 ){ pType=SysParType::Wlnu; return 5; }
+ if( MCID==364169 ){ pType=SysParType::Wlnu; return 7; }
+
+ //Check Wtaunu 2.2.1 samples
+ if( MCID>=364184 && MCID<364187 ){ pType=SysParType::Wlnu; return 1; }
+ if( MCID>=364187 && MCID<364190 ){ pType=SysParType::Wlnu; return 2; }
+ if( MCID>=364190 && MCID<364193 ){ pType=SysParType::Wlnu; return 3; }
+ if( MCID>=364193 && MCID<364196 ){ pType=SysParType::Wlnu; return 4; }
+ if( MCID==364196 ){ pType=SysParType::Wlnu; return 5; }
+ if( MCID==364197 ){ pType=SysParType::Wlnu; return 7; }
+
+
return -1;
}
diff --git a/PhysicsAnalysis/AnalysisCommon/PMGTools/Root/PMGTruthWeightTool.cxx b/PhysicsAnalysis/AnalysisCommon/PMGTools/Root/PMGTruthWeightTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..fc1f8fccc374d93c89f8ab3785bf064c886e8a89
--- /dev/null
+++ b/PhysicsAnalysis/AnalysisCommon/PMGTools/Root/PMGTruthWeightTool.cxx
@@ -0,0 +1,397 @@
+/*
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+// EDM include(s):
+#ifndef XAOD_STANDALONE
+ #include "AthAnalysisBaseComps/AthAnalysisHelper.h"
+ #include "EventInfo/EventInfo.h"
+ #include "EventInfo/EventType.h"
+#endif
+
+// Local include(s):
+#include <CxxUtils/StringUtils.h>
+#include <RootCoreUtils/StringUtil.h>
+
+#include <PATInterfaces/SystematicRegistry.h>
+#include <xAODEventInfo/EventInfo.h>
+#include <xAODMetaData/FileMetaData.h>
+
+#include <PMGTools/PMGTruthWeightTool.h>
+
+// For replacing substrings
+#include <boost/algorithm/string/replace.hpp>
+#include <boost/algorithm/string/case_conv.hpp>
+
+namespace PMGTools
+{
+ PMGTruthWeightTool::PMGTruthWeightTool(const std::string& name)
+ : asg::AsgMetadataTool(name)
+ {
+ declareProperty("MetaObjectName", m_metaName = "TruthMetaData");
+ }
+
+
+ StatusCode PMGTruthWeightTool::initialize()
+ {
+ ATH_MSG_DEBUG("Initialising...");
+
+ ATH_MSG_INFO("Attempting to retrieve truth meta data from the first file...");
+
+ // Clear cached weights
+ clearWeightLocationCaches();
+
+ // Set the MC channel number to an impossible number (max-uint) to force meta data updating
+ m_mcChannelNumber = uint32_t(-1);
+
+ // Try to load MC channel number from file metadata
+ ATH_MSG_INFO("Attempting to retrieve MC channel number...");
+ const xAOD::FileMetaData *fmd = nullptr;
+ if (inputMetaStore()->contains<xAOD::FileMetaData>("FileMetaData")) {
+ ATH_CHECK(inputMetaStore()->retrieve(fmd, "FileMetaData"));
+ float fltChannelNumber(-1);
+ if (fmd->value(xAOD::FileMetaData::mcProcID, fltChannelNumber)) {
+ m_mcChannelNumber = uint32_t(fltChannelNumber);
+ }
+ }
+ if (m_mcChannelNumber == uint32_t(-1)) {
+ ATH_MSG_WARNING("... MC channel number could not be loaded");
+ } else {
+ ATH_MSG_INFO("... MC channel number identified as " << m_mcChannelNumber);
+ }
+
+ // Start by trying to load metadata from the store
+ m_metaDataContainer = nullptr;
+ if (inputMetaStore()->contains<xAOD::TruthMetaDataContainer>(m_metaName)) {
+ ATH_CHECK(inputMetaStore()->retrieve(m_metaDataContainer, m_metaName));
+ ATH_MSG_INFO("Loaded xAOD::TruthMetaDataContainer");
+
+ // Check for incorrectly stored McChannelNumber
+ m_useChannelZeroInMetaData = true;
+ for (auto truthMetaData : *m_metaDataContainer) {
+ if (truthMetaData->mcChannelNumber() != 0) { m_useChannelZeroInMetaData = false; }
+ }
+ // If we only have one metadata item take MC channel from there if needed
+ if (m_mcChannelNumber == uint32_t(-1) && m_metaDataContainer->size() == 1) {
+ m_mcChannelNumber = m_metaDataContainer->at(0)->mcChannelNumber();
+ ATH_MSG_WARNING("... MC channel number taken from the metadata as " << m_mcChannelNumber);
+ }
+ if (m_useChannelZeroInMetaData) { ATH_MSG_WARNING("MC channel number in TruthMetaData is invalid - assuming that channel 0 has the correct information."); }
+
+ // Load metadata from TruthMetaDataContainer if we have a valid channel number or if we're going to use 0 anyway
+ // ... otherwise wait until we can load a channel number from EventInfo
+ if (m_mcChannelNumber != uint32_t(-1) || m_useChannelZeroInMetaData) {
+ if (loadMetaData().isFailure()) {
+ ATH_MSG_ERROR("Could not load metadata for MC channel number " << m_mcChannelNumber);
+ return StatusCode::FAILURE;
+ }
+ }
+ } else {
+ // ... now try to load the weight container using the POOL metadata (not possible in AnalysisBase)
+ // see https://twiki.cern.ch/twiki/bin/viewauth/AtlasProtected/AthAnalysisBase#How_to_read_the_truth_weight_nam
+ if (loadPOOLMetaData().isFailure()) {
+ ATH_MSG_ERROR("Could not load POOL HepMCWeightNames");
+ return StatusCode::FAILURE;
+ }
+ }
+
+ // Add the affecting systematics to the global registry
+ CP::SystematicRegistry& registry = CP::SystematicRegistry::getInstance();
+ if (!registry.registerSystematics(*this)) {
+ ATH_MSG_ERROR("unable to register the systematics");
+ return StatusCode::FAILURE;
+ }
+
+ ATH_MSG_DEBUG("Successfully initialized!");
+
+ return StatusCode::SUCCESS;
+ }
+
+
+ const std::vector<std::string>& PMGTruthWeightTool::getWeightNames() const {
+ return m_weightNames;
+ }
+
+
+ float PMGTruthWeightTool::getWeight(const std::string& weightName) const {
+ // Check that we have loaded event weights
+ if (m_weights.size() == 0) {
+ ATH_MSG_FATAL("No MC weights found!");
+ throw std::runtime_error(name() + ": No MC weights found!");
+ }
+
+ // Return appropriate weight from EventInfo: this should be identical to the TruthEvent
+ try {
+ return m_weights.at(m_weightIndices.at(weightName));
+ } catch (const std::out_of_range& e) {
+ // Before throwing an exception, try to recover with bad naming conventions
+ std::string strippedName = boost::algorithm::to_lower_copy(boost::replace_all_copy(weightName, " ", ""));
+ for (const auto & weight:m_weightNames){
+ if (strippedName==boost::algorithm::to_lower_copy(boost::replace_all_copy(weight," ", ""))){
+ ATH_MSG_WARNING("Using weight name \"" << weight << "\" instead of requested \"" << weightName << "\"");
+ return getWeight(weight);
+ }
+ }
+ ATH_MSG_FATAL("Weight \"" + weightName + "\" could not be found");
+ throw std::runtime_error(name() + ": Weight \"" + weightName + "\" could not be found");
+ }
+ }
+
+
+ bool PMGTruthWeightTool::hasWeight(const std::string& weightName) const {
+ return (m_weightIndices.count(weightName) > 0);
+ }
+
+
+ float PMGTruthWeightTool::getSysWeight() const
+ {
+ return m_currentWeightData->weight;
+ }
+
+
+ size_t PMGTruthWeightTool::getSysWeightIndex() const
+ {
+ return m_currentWeightData->index;
+ }
+
+
+ bool PMGTruthWeightTool::isAffectedBySystematic(const CP::SystematicVariation& systematic) const
+ {
+ return m_systematicsSet.find( systematic ) != m_systematicsSet.end();
+ }
+
+
+ CP::SystematicSet PMGTruthWeightTool::affectingSystematics() const
+ {
+ return m_systematicsSet;
+ }
+
+
+ CP::SystematicSet PMGTruthWeightTool::recommendedSystematics() const
+ {
+ return affectingSystematics();
+ }
+
+
+ CP::SystematicCode PMGTruthWeightTool::applySystematicVariation(const CP::SystematicSet& systConfig)
+ {
+ auto iter = m_weightData.find (systConfig);
+ if (iter != m_weightData.end())
+ {
+ m_currentWeightData = &iter->second;
+ return CP::SystematicCode::Ok;
+ }
+
+ CP::SystematicSet currentSys;
+ ANA_CHECK_SET_TYPE (CP::SystematicCode);
+ ANA_CHECK (CP::SystematicSet::filterForAffectingSystematics (systConfig, m_systematicsSet, currentSys));
+
+ WeightData currentWeight{};
+ currentWeight.index = m_weightIndicesSys.at(currentSys.name());
+ if (!m_weights.empty())
+ {
+ currentWeight.weight = m_weights.at(currentWeight.index);
+ }
+
+ auto insert = m_weightData.emplace(systConfig, std::move(currentWeight));
+ m_currentWeightData = &insert.first->second;
+
+ return CP::SystematicCode::Ok;
+ }
+
+
+ StatusCode PMGTruthWeightTool::beginInputFile()
+ {
+ // Detect possible MC channel number change
+ uint32_t mcChannelNumber(-1);
+
+ // Try to load MC channel number from file metadata
+ ATH_MSG_DEBUG("Attempting to retrieve MC channel number...");
+ const xAOD::FileMetaData *fmd = nullptr;
+ if (inputMetaStore()->contains<xAOD::FileMetaData>("FileMetaData")) {
+ ATH_CHECK(inputMetaStore()->retrieve(fmd, "FileMetaData"));
+ float fltChannelNumber(-1);
+ if (fmd->value(xAOD::FileMetaData::mcProcID, fltChannelNumber)) {
+ mcChannelNumber = uint32_t(fltChannelNumber);
+ }
+ }
+ if (m_mcChannelNumber != uint32_t(-1) && mcChannelNumber != uint32_t(-1) && mcChannelNumber != m_mcChannelNumber) {
+ ATH_MSG_ERROR("MC channel number from a new file does not match the previously processed files.");
+ return StatusCode::FAILURE;
+ }
+
+ return StatusCode::SUCCESS;
+ }
+
+
+ StatusCode PMGTruthWeightTool::beginEvent() {
+ // Clear weights and channel number from previous event
+ m_weights.clear();
+ m_weightData.clear();
+
+ // Try to read information about weights and channel number from EventInfo
+ // 1. try the xAOD::EventInfo object
+ if (evtStore()->contains<xAOD::EventInfo>("EventInfo")) {
+ const xAOD::EventInfo* evtInfo;
+ ATH_CHECK(evtStore()->retrieve(evtInfo, "EventInfo"));
+ m_weights = evtInfo->mcEventWeights();
+
+ // 2. otherwise, if we're not in AnalysisBase, see if there's an EventInfo object
+#ifndef XAOD_STANDALONE
+ } else if (evtStore()->contains<EventInfo>("McEventInfo")) {
+ const EventInfo* evtInfo;
+ ATH_CHECK(evtStore()->retrieve(evtInfo, "McEventInfo"));
+ for (unsigned int idx = 0; idx < evtInfo->event_type()->n_mc_event_weights(); ++idx) {
+ m_weights.push_back(evtInfo->event_type()->mc_event_weight(idx));
+ }
+#endif // not XAOD_STANDALONE
+
+ // 3. otherwise let's bail here
+ } else {
+ ATH_MSG_ERROR("No event information is available in this file!");
+ return StatusCode::FAILURE;
+ }
+
+ // Validate weight caches against event information
+ if (m_weightNames.size() != m_weights.size()) {
+ // Special case to allow nominal weight to be used when this is the only weight in the event
+ if ((m_weightNames.size() == 0) && (m_weights.size() == 1)) {
+ ATH_MSG_WARNING("No weight names were available in this sample! Proceeding under the assumption that the single available weight should be 'nominal'");
+ m_weightNames.push_back("nominal");
+ m_weightIndices["nominal"] = 0;
+ m_weightIndicesSys[""] = 0;
+ } else {
+ ATH_MSG_ERROR("Expected " << m_weightNames.size() << " weights from the metadata but found " << m_weights.size() << " in this event");
+ ATH_MSG_ERROR("Perhaps this sample was made using a release which did not correctly propagate the event weights.");
+ }
+ }
+
+ // Apply nominal systematics by default
+ ANA_CHECK (applySystematicVariation (CP::SystematicSet()));
+
+ return StatusCode::SUCCESS;
+ }
+
+
+ StatusCode PMGTruthWeightTool::loadMetaData() {
+ // Clear cached weight names and positions
+ this->clearWeightLocationCaches();
+
+ // Find the correct truth meta data object
+ uint32_t targetChannelNumber = (m_useChannelZeroInMetaData ? 0 : m_mcChannelNumber);
+ ATH_MSG_INFO("Attempting to load weight meta data from xAOD TruthMetaData for channel " << targetChannelNumber);
+ auto itTruthMetaDataPtr = std::find_if(m_metaDataContainer->begin(), m_metaDataContainer->end(),
+ [targetChannelNumber] (const auto& it) { return it->mcChannelNumber() == targetChannelNumber; }
+ );
+
+ // If no such object is found then return
+ if (itTruthMetaDataPtr == m_metaDataContainer->end()) {
+ ATH_MSG_ERROR("Could not load weight meta data!");
+ return StatusCode::FAILURE;
+ }
+
+ // Update cached weight data
+ const std::vector<std::string> &truthWeightNames = (*itTruthMetaDataPtr)->weightNames();
+ for(std::size_t idx = 0; idx < truthWeightNames.size(); ++idx ) {
+ ANA_MSG_VERBOSE(" " << truthWeightNames.at(idx));
+ m_weightNames.push_back(truthWeightNames.at(idx));
+ m_weightIndices[truthWeightNames.at(idx)] = idx;
+
+ std::string sysName = weightNameToSys(truthWeightNames.at(idx));
+ if (!sysName.empty()) {
+ m_systematicsSet.insert(CP::SystematicVariation(sysName));
+ }
+ m_weightIndicesSys[sysName] = idx;
+ }
+ return this->validateWeightLocationCaches();
+ }
+
+
+ StatusCode PMGTruthWeightTool::loadPOOLMetaData() {
+ // AnalysisBase can only use the xAOD::TruthMetaDataContainer, so skip this
+#ifdef XAOD_STANDALONE
+ return StatusCode::SUCCESS;
+#else
+ ATH_MSG_INFO("Looking for POOL HepMC IOVMetaData...");
+ std::map<std::string, int> hepMCWeightNamesMap;
+ if (AAH::retrieveMetadata("/Generation/Parameters", "HepMCWeightNames", hepMCWeightNamesMap, inputMetaStore()).isFailure()) {
+ ATH_MSG_FATAL("Cannot access metadata " << m_metaName << " and failed to get names from IOVMetadata");
+ return StatusCode::FAILURE;
+ }
+
+ // Clear cached weight names and positions
+ this->clearWeightLocationCaches();
+
+ // Use input map to fill the index map and the weight names
+ ATH_MSG_INFO("Attempting to load weight meta data from HepMC IOVMetaData container");
+ for (auto& kv : hepMCWeightNamesMap) {
+ ANA_MSG_VERBOSE(" " << kv.first);
+ m_weightNames.push_back(kv.first);
+ m_weightIndices[kv.first] = kv.second;
+
+ std::string sysName = weightNameToSys(kv.first);
+ if (!sysName.empty()) {
+ m_systematicsSet.insert(CP::SystematicVariation(sysName));
+ }
+ m_weightIndicesSys[sysName] = kv.second;
+ }
+ return this->validateWeightLocationCaches();
+#endif // XAOD_STANDALONE
+ }
+
+
+ StatusCode PMGTruthWeightTool::validateWeightLocationCaches() {
+ // Validate weight caches against one another
+ if (m_weightNames.size() != m_weightIndices.size()) {
+ ATH_MSG_ERROR("Found " << m_weightNames.size() << " but " << m_weightIndices.size() << " weight indices!");
+ return StatusCode::FAILURE;
+ }
+ // Check if we can work with systematics
+ auto it = m_weightIndicesSys.find("");
+ if (it == m_weightIndicesSys.end()) {
+ ATH_MSG_WARNING("Could not detect nominal weight automatically. The first weight will also be considered nominal.");
+ m_weightIndicesSys[""] = 0;
+ }
+
+ ATH_MSG_INFO("Successfully loaded information about " << m_weightNames.size() << " weights");
+ return StatusCode::SUCCESS;
+ }
+
+
+ void PMGTruthWeightTool::clearWeightLocationCaches() {
+ m_weightNames.clear();
+ m_weightIndices.clear();
+ m_weightIndicesSys.clear();
+ m_systematicsSet.clear();
+ }
+
+
+ std::string PMGTruthWeightTool::weightNameToSys (const std::string &name) const
+ {
+ if (name.empty()) // empty weight is nominal
+ {
+ return "";
+ }
+
+ // Trim trailing whitespace
+ std::string sys = CxxUtils::StringUtils::trim (name);
+ if (sys == "nominal" // Powheg calls it "nominal"
+ || sys == "Weight") // Sherpa names the nominal weight just "Weight"
+ {
+ return "";
+ }
+
+ sys = RCU::substitute (sys, " set = ", "_"); // Powheg
+ sys = RCU::substitute (sys, " = ", "_"); // Powheg
+ sys = RCU::substitute (sys, "=", "");
+ sys = RCU::substitute (sys, ",", "");
+ sys = RCU::substitute (sys, ".", "");
+ sys = RCU::substitute (sys, ":", "");
+ sys = RCU::substitute (sys, " ", "_");
+ sys = RCU::substitute (sys, "#", "num");
+ sys = RCU::substitute (sys, "\n", "_");
+ sys = RCU::substitute (sys, "/", "over"); // MadGraph
+
+ return generatorSystematicsPrefix + sys;
+ }
+} // namespace PMGTools
diff --git a/PhysicsAnalysis/AnalysisCommon/PMGTools/python/test_PMGCrossSectionTool.py b/PhysicsAnalysis/AnalysisCommon/PMGTools/python/test_PMGCrossSectionTool.py
index 5de2f5f4e487ee0c049cc8d7ad6e66fa2c6847d0..3b4b493acf491115fe49df4962b4bf5e51278df1 100644
--- a/PhysicsAnalysis/AnalysisCommon/PMGTools/python/test_PMGCrossSectionTool.py
+++ b/PhysicsAnalysis/AnalysisCommon/PMGTools/python/test_PMGCrossSectionTool.py
@@ -1,24 +1,34 @@
#!/bin/env python
# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-import ROOT, os
+import ROOT
-#def main(filename, **args):
-def main(**args):
+def main():
+ from PathResolver import PathResolver
tool = ROOT.PMGTools.PMGCrossSectionTool('MyXSectionTool')
- tool.readInfosFromDir(os.getenv('ROOTCOREBIN')+'/data/PMGTools/')
-
- #take a ttbar sample as example ( user should get this from the EventInfo of course )
+ fn = '/eos/atlas/atlascerngroupdisk/asg-calib/dev/PMGTools/PMGxsecDB_mc16.txt'
+ fn = PathResolver.FindCalibFile(fn)
+ vv = ROOT.std.vector('std::string')()
+ vv.push_back(fn)
+
+ tool.readInfosFromFiles(vv)
+
+ # take a ttbar sample as example ( users should get this from the EventInfo )
sample_id = 410000
+ print '%d sample loaded' % tool.getLoadedDSIDs().size()
print
- print 'sample name = ',tool.getSampleName(sample_id)
- print 'xsection [pb] = ',tool.getSampleXsection(sample_id)
- print 'filter eff = ',tool.getFilterEff(sample_id)
- print 'branching ratio = ',tool.getBR(sample_id)
- print 'k factor = ',tool.getKfactor(sample_id)
+ print 'Sample dsid = ', sample_id
+ print 'Sample name = ', tool.getSampleName(sample_id)
+ print 'xsection [pb] = ', tool.getSampleXsection(sample_id)
+ print 'filter eff = ', tool.getFilterEff(sample_id)
+ print 'k factor = ', tool.getKfactor(sample_id)
+ print 'xsection uncertainty = ', tool.getXsectionUncertainty(sample_id)
+ print 'xsection uncertainty up = ', tool.getXsectionUncertaintyUP(sample_id)
+ print 'xsection uncertainty down = ', tool.getXsectionUncertaintyDOWN(sample_id)
+
print
if __name__ == '__main__':
diff --git a/PhysicsAnalysis/AnalysisCommon/PMGTools/src/components/PMGTools_entries.cxx b/PhysicsAnalysis/AnalysisCommon/PMGTools/src/components/PMGTools_entries.cxx
index 11d64b3b214502c0077ce62006ec844809e3dda0..f7ce363ba052d72156fb7efeb01ffe06453cb1d8 100644
--- a/PhysicsAnalysis/AnalysisCommon/PMGTools/src/components/PMGTools_entries.cxx
+++ b/PhysicsAnalysis/AnalysisCommon/PMGTools/src/components/PMGTools_entries.cxx
@@ -1,10 +1,14 @@
+
+#include "PMGTools/PMGCrossSectionTool.h"
+#include "PMGTools/PMGDecayProductsSelectionTool.h"
#include "PMGTools/PMGSherpa22VJetsWeightTool.h"
#include "PMGTools/PMGSherpaVjetsSysTool.h"
-#include "PMGTools/PMGCrossSectionTool.h"
+#include "PMGTools/PMGTruthWeightTool.h"
using namespace PMGTools;
+DECLARE_COMPONENT( PMGCrossSectionTool )
+DECLARE_COMPONENT( PMGDecayProductsSelectionTool )
DECLARE_COMPONENT( PMGSherpa22VJetsWeightTool )
DECLARE_COMPONENT( PMGSherpaVjetsSysTool )
-DECLARE_COMPONENT( PMGCrossSectionTool )
-
+DECLARE_COMPONENT( PMGTruthWeightTool )
diff --git a/PhysicsAnalysis/AnalysisCommon/PMGTools/src/components/PMGTools_load.cxx b/PhysicsAnalysis/AnalysisCommon/PMGTools/src/components/PMGTools_load.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..cec1175f2c6a4fa2044b0eb24c6aa440c489386b
--- /dev/null
+++ b/PhysicsAnalysis/AnalysisCommon/PMGTools/src/components/PMGTools_load.cxx
@@ -0,0 +1,3 @@
+
+#include "GaudiKernel/LoadFactoryEntries.h"
+LOAD_FACTORY_ENTRIES( PMGTools )
diff --git a/PhysicsAnalysis/AnalysisCommon/PMGTools/test/MyPMGApp.cxx b/PhysicsAnalysis/AnalysisCommon/PMGTools/test/MyPMGApp.cxx
index bcd170a2b58f05945351b142a31583cc88ce36bf..7398015984b708c78e708ffe948ed2c8bd41cbf9 100644
--- a/PhysicsAnalysis/AnalysisCommon/PMGTools/test/MyPMGApp.cxx
+++ b/PhysicsAnalysis/AnalysisCommon/PMGTools/test/MyPMGApp.cxx
@@ -2,183 +2,175 @@
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
-#ifndef MYPACKAGE_MYAPP_H
-#define MYPACKAGE_MYAPP_H
-#ifndef ROOTCORE
+/// Example standalone executable using POOL to read an xAOD
+/// Tests several tools, showing example usage
-#include "AthAnalysisBaseComps/AthAnalysisHelper.h" //tool creation and configuration
+#ifndef PMGTOOLS_MYPMGAPP_H
+#define PMGTOOLS_MYPMGAPP_H
-
-#include "POOLRootAccess/TEvent.h" //event looping
-#include "GaudiKernel/ToolHandle.h" //for better working with tools
-
-// Tool include(s):
+// EDM includes
#include "AsgTools/AnaToolHandle.h"
-
-#include "AsgTools/MessageCheck.h" //messaging
-using namespace asg::msgUserCode; //messaging
-
-// added for truth studies
-#include "xAODJet/JetContainer.h"
-
+#include "AsgTools/MessageCheck.h" // for messaging
+#include "POOLRootAccess/TEvent.h" // event looping
+#include "xAODJet/JetContainer.h" // for jet studies
+#include "GaudiKernel/ToolHandle.h" // for better working with tools
#include "PATInterfaces/IWeightTool.h"
-#include "../Root/PMGSherpa22VJetsWeightTool.cxx"
-//#PMGSherpa22VJetsWeightTooll
-//ROOT includes
+// ROOT includes
#include "TString.h"
#include "TSystem.h"
#include "TH1F.h"
#include "TFile.h"
-int main( int argc, char* argv[] ) {
+// Local includes
+#include "PMGTools/PMGSherpa22VJetsWeightTool.h"
+#include "PMGTools/PMGTruthWeightTool.h"
+
+// For convenience messaging macros
+using namespace asg::msgUserCode;
+int main(int argc, char *argv[])
+{
+ ANA_CHECK_SET_TYPE (int); // makes ANA_CHECK return ints if exiting function
+
// for access to everything including EVNT
- IAppMgrUI* app = POOL::Init(); //important to do this first!
- // just for xAOD file type access - improved performance due to
+ IAppMgrUI *app = POOL::Init(); // important to do this first!
+ // just for xAOD file type access - improved performance due to
// fast xAOD reading mechanism rather than generic POOL reading mechanism
- //IAppMgrUI* app = POOL::Init("POOLRootAccess/basixAOD.opts");
-
-
- ANA_MSG_INFO("Making the tool");
-
- // test the interface
- asg::AnaToolHandle< IWeightTool > pmgTool( "PMGSherpa22VJetsWeightTool" );
- ANA_CHECK( pmgTool.setProperty( "TruthJetContainer", "AntiKt4TruthWZJets" ) ); // default
- ANA_CHECK( pmgTool.initialize() );
-
- // can make the subtool directly or via this cast
- ANA_MSG_INFO("Making the sub tool");
- // PMGTools::PMGSherpa22VJetsWeightTool sherpatool;
- //sherpatool.initialize();
- PMGTools::PMGSherpa22VJetsWeightTool* sherpatool = dynamic_cast<PMGTools::PMGSherpa22VJetsWeightTool*>(&*pmgTool);
- // PMGCrossSectionTool* xsectool = dynamic_cast<PMGCrossSectionTool*>(&*pmgXsecTool);
-
- // Open the input file:
- TString fileName = "$ASG_TEST_FILE_MC";
- if( argc < 2 ) {
- ANA_MSG_WARNING( "No file name received, using $ASG_TEST_FILE_MC" );
- } else {
- fileName = argv[1]; //use the user provided file
- }
-
- ANA_MSG_INFO("Opening file: " << gSystem->ExpandPathName(fileName.Data()) );
-
-
- //loop over input file with POOL
- POOL::TEvent evt;
- evt.readFrom( fileName );
-
- ANA_MSG_INFO("Opening file: " << evt.getEntries() );
-
- // book some histograms
- TH1F* h_njetTruth = new TH1F("jetmult_AntiKt4TruthJets","jetmult_AntiKt4TruthJets",10,-0.5,9.5);
- TH1F* h_njetTruthWZ = new TH1F("jetmult_AntiKt4TruthWZJets","jetmult_AntiKt4TruthWZJets",10,-0.5,9.5);
-
- //
- TH1F* h_njetTruth_Tool = new TH1F("jetmult_AntiKt4TruthJets_Tool","jetmult_AntiKt4TruthJets_Tool",10,-0.5,9.5);
- TH1F* h_njetTruthWZ_Tool = new TH1F("jetmult_AntiKt4TruthWZJets_Tool","jetmult_AntiKt4TruthWZJets_Tool",10,-0.5,9.5);
-
-
- TFile* histfile=TFile::Open("hists.root","RECREATE");
-
- bool debug=false;
-
- for(int i=0;i < evt.getEntries(); i++) {
- // for(int i=0;i < 200; i++) {
- if( evt.getEntry(i) < 0) { ANA_MSG_ERROR("Failed to read event " << i); continue; }
-
- const xAOD::JetContainer* truthWZJets=0;
- evt.retrieve(truthWZJets, "AntiKt4TruthWZJets" );
-
- int nTruthWZJets20 = 0;
- int nTruthWZJets30 = 0;
- //loop over all AntiKtTruthWZ jets
- for (unsigned i=0; i!=truthWZJets->size(); i++) {
- const xAOD::Jet* truthjet = (*truthWZJets)[i];
- //if(m_debug) {std::cout<<"Index "<<i<<" truthjet pointer = "<<truthjet<<std::endl;}
- double Pt=truthjet->pt();
- double AbsEta=fabs(truthjet->eta());
- if (Pt>20000.0 && AbsEta<4.5) {nTruthWZJets20++;}
- if (Pt>30000.0 && AbsEta<4.5) {nTruthWZJets30++;}
- }//end loop over truth jets
- if (debug) {
- std::cout<<"nTruthWZJets20="<<nTruthWZJets20<<std::endl;
- std::cout<<"nTruthWZJets30="<<nTruthWZJets30<<std::endl;
- }
- h_njetTruthWZ->Fill(nTruthWZJets20);
-
- // same for AntiKt4TruthJets
- const xAOD::JetContainer* truthJets=0;
- evt.retrieve(truthJets, "AntiKt4TruthJets" );
-
- int nTruthJets20 = 0;
- //loop over all truth jets
- for (unsigned i=0; i!=truthJets->size(); i++) {
- const xAOD::Jet* truthjet = (*truthJets)[i];
- double Pt=truthjet->pt();
- double AbsEta=fabs(truthjet->eta());
- if (Pt>20000.0 && AbsEta<4.5) {nTruthJets20++;}
- }//end loop over truth jets
- h_njetTruth->Fill(nTruthJets20);
-
- // Now use the PMG tool directly over AntiKt4TruthWZJets
- int njet=sherpatool->getSherpa22VJets_NJet("AntiKt4TruthWZJets");
-
- if(debug)
- std::cout<<"njets TruthWZJets "<< njet << std::endl;
-
- h_njetTruthWZ_Tool->Fill(njet);
-
- double reweight(0);
-
- int nuse=njet;
- reweight=sherpatool->getSherpa22VJets_NJetCorrection(nuse) ;
- if (debug)
- std::cout<<"correction with njet input(check) "<< nuse << " " << reweight << std::endl;
-
- // check
- reweight=sherpatool->getSherpa22VJets_NJetCorrection("AntiKt4TruthWZJets");
- if (debug)
- std::cout<<"correction with string "<< reweight << std::endl;
-
-
- // and do again for TruthJets
- njet=sherpatool->getSherpa22VJets_NJet("AntiKt4TruthJets");
- if (debug)
- std::cout<<"njets (TruthJets) "<< njet << std::endl;
-
- h_njetTruth_Tool->Fill(njet);
-
- nuse=njet;
- reweight=sherpatool->getSherpa22VJets_NJetCorrection(nuse) ;
- if(debug)
- std::cout<<"correction with njet input (TruthJets) (check) "<< nuse << " " << reweight << std::endl;
-
- // check
- reweight=sherpatool->getSherpa22VJets_NJetCorrection("AntiKt4TruthJets");
- if(debug)
- std::cout<<"correction with string (TruthJets) "<< reweight << std::endl;
-
- // and test the interface itself
- reweight=sherpatool->getWeight();
- if(debug)
- std::cout<<"correction from Interface" << reweight << std::endl;
-
- }
-
- app->finalize(); //trigger finalization of all services and tools created by the Gaudi Application
-
- histfile->cd();
- h_njetTruthWZ->Write();
- h_njetTruthWZ_Tool->Write();
- h_njetTruth->Write();
- h_njetTruth_Tool->Write();
- histfile->Close();
-
- return 0;
+ // IAppMgrUI* app = POOL::Init("POOLRootAccess/basixAOD.opts");
+
+ // test the interface
+ ANA_MSG_INFO("Creating the PMGSherpa22VJetsWeightTool...");
+ asg::AnaToolHandle< IWeightTool > pmgTool("PMGTools::PMGSherpa22VJetsWeightTool/PMGSherpa22VJetsWeightTool");
+ ANA_CHECK(pmgTool.setProperty("TruthJetContainer", "AntiKt4TruthWZJets")); // default
+ ANA_CHECK(pmgTool.initialize());
+ // can make the subtool directly or via this cast
+ ANA_MSG_INFO("Casting to PMGSherpa22VJetsWeightTool...");
+ PMGTools::PMGSherpa22VJetsWeightTool* sherpaTool = dynamic_cast<PMGTools::PMGSherpa22VJetsWeightTool*>(&*pmgTool);
+
+ // Create the truth weight tool:
+ // ... could create directly with
+ // ANA_MSG_INFO("Creating the PMGTruthWeightTool...");
+ // PMGTools::PMGTruthWeightTool truthWeightTool("PMGTruthWeightTool");
+ // ANA_CHECK(truthWeightTool.setProperty("OutputLevel", MSG::INFO));
+ // ANA_CHECK(truthWeightTool.sysInitialize()); // must call sysInitialize to get the callbacks registered properly for an AsgMetadataTool
+ // ... but better to do this through a ToolHandle
+ asg::AnaToolHandle< PMGTools::IPMGTruthWeightTool > truthWeightTool("PMGTools::PMGTruthWeightTool/PMGTruthWeightTool");
+ ANA_CHECK(truthWeightTool.initialize());
+
+ // Open the input file:
+ TString fileName = "$ASG_TEST_FILE_MC";
+ if (argc < 2) {
+ ANA_MSG_WARNING("No file name received, using $ASG_TEST_FILE_MC");
+ } else {
+ fileName = argv[1]; // use the user provided file
+ }
+ ANA_MSG_INFO("Opening file: " << gSystem->ExpandPathName(fileName.Data()));
+
+
+ // loop over input file with POOL
+ POOL::TEvent evt;
+ ANA_CHECK(evt.readFrom(fileName));
+ ANA_MSG_INFO("... with " << evt.getEntries() << " entries");
+
+ // book some histograms
+ TH1F *h_njetTruth = new TH1F("jetmult_AntiKt4TruthJets", "jetmult_AntiKt4TruthJets", 10, -0.5, 9.5);
+ TH1F *h_njetTruthWZ = new TH1F("jetmult_AntiKt4TruthWZJets", "jetmult_AntiKt4TruthWZJets", 10, -0.5, 9.5);
+ TH1F *h_njetTruth_Tool = new TH1F("jetmult_AntiKt4TruthJets_Tool", "jetmult_AntiKt4TruthJets_Tool", 10, -0.5, 9.5);
+ TH1F *h_njetTruthWZ_Tool = new TH1F("jetmult_AntiKt4TruthWZJets_Tool", "jetmult_AntiKt4TruthWZJets_Tool", 10, -0.5, 9.5);
+ TFile *histfile = TFile::Open("hists.root", "RECREATE");
+
+ bool debug = false;
+
+ for (int i = 0; i < evt.getEntries(); i++) {
+ if (evt.getEntry(i) < 0) { ANA_MSG_ERROR("Failed to read event " << i); continue; }
+
+ const xAOD::JetContainer *truthWZJets = 0;
+ ANA_CHECK(evt.retrieve(truthWZJets, "AntiKt4TruthWZJets"));
+ int nTruthWZJets20 = 0;
+ int nTruthWZJets30 = 0;
+
+ // loop over all AntiKtTruthWZ jets
+ for (unsigned i = 0; i != truthWZJets->size(); i++) {
+ const xAOD::Jet *truthjet = (*truthWZJets)[i];
+ // if(m_debug) { std::cout << "Index " << i << " truthjet pointer = " << truthjet << std::endl; }
+ double Pt = truthjet->pt();
+ double AbsEta = fabs(truthjet->eta());
+
+ if (Pt > 20000.0 && AbsEta < 4.5) {nTruthWZJets20++;}
+ if (Pt > 30000.0 && AbsEta < 4.5) {nTruthWZJets30++;}
+ } //end loop over truth jets
+
+ if (debug) {
+ std::cout << "nTruthWZJets20=" << nTruthWZJets20 << std::endl;
+ std::cout << "nTruthWZJets30=" << nTruthWZJets30 << std::endl;
+ }
+ h_njetTruthWZ->Fill(nTruthWZJets20);
+
+ // same for AntiKt4TruthJets
+ const xAOD::JetContainer *truthJets = 0;
+ ANA_CHECK(evt.retrieve(truthJets, "AntiKt4TruthJets"));
+ int nTruthJets20 = 0;
+
+ // loop over all truth jets
+ for (unsigned i = 0; i != truthJets->size(); i++) {
+ const xAOD::Jet *truthjet = (*truthJets)[i];
+ double Pt = truthjet->pt();
+ double AbsEta = fabs(truthjet->eta());
+
+ if (Pt > 20000.0 && AbsEta < 4.5) {nTruthJets20++;}
+ } // end loop over truth jets
+ h_njetTruth->Fill(nTruthJets20);
+
+ // Now use the PMG tool directly over AntiKt4TruthWZJets
+ int njet = sherpaTool->getSherpa22VJets_NJet("AntiKt4TruthWZJets");
+ if (debug) { std::cout << "njets TruthWZJets " << njet << std::endl; }
+ h_njetTruthWZ_Tool->Fill(njet);
+
+ double reweight(0);
+ int nuse = njet;
+ reweight = sherpaTool->getSherpa22VJets_NJetCorrection(nuse);
+ if (debug) { std::cout << "correction with njet input(check) " << nuse << " " << reweight << std::endl; }
+
+ // check
+ reweight = sherpaTool->getSherpa22VJets_NJetCorrection("AntiKt4TruthWZJets");
+ if (debug) { std::cout << "correction with string " << reweight << std::endl; }
+
+ // and do again for TruthJets
+ njet = sherpaTool->getSherpa22VJets_NJet("AntiKt4TruthJets");
+ if (debug) { std::cout << "njets (TruthJets) " << njet << std::endl; }
+ h_njetTruth_Tool->Fill(njet);
+
+ nuse = njet;
+ reweight = sherpaTool->getSherpa22VJets_NJetCorrection(nuse) ;
+ if (debug) { std::cout << "correction with njet input (TruthJets) (check) " << nuse << " " << reweight << std::endl; }
+
+ // check
+ reweight = sherpaTool->getSherpa22VJets_NJetCorrection("AntiKt4TruthJets");
+
+ if (debug) { std::cout << "correction with string (TruthJets) " << reweight << std::endl; }
+
+ // and test the interface itself
+ reweight = sherpaTool->getWeight();
+ if (debug) { std::cout << "correction from Interface" << reweight << std::endl; }
+
+ // Test the PMGTruthWeightTool interface
+ auto weightNames = truthWeightTool->getWeightNames();
+ ANA_MSG_INFO("Event #" << i << ": found " << weightNames.size() << " weights for this event");
+
+ }
+
+ ANA_CHECK(app->finalize()); // trigger finalization of all services and tools created by the Gaudi Application
+
+ histfile->cd();
+ h_njetTruthWZ->Write();
+ h_njetTruthWZ_Tool->Write();
+ h_njetTruth->Write();
+ h_njetTruth_Tool->Write();
+ histfile->Close();
+
+ return 0;
}
-#endif // ROOTCORE
-#endif //> !MYPACKAGE_MYAPP_H
+#endif //> !PMGTOOLS_MYPMGAPP_H
diff --git a/PhysicsAnalysis/AnalysisCommon/PMGTools/test/ut_PMGTruthWeightTool_test.cxx b/PhysicsAnalysis/AnalysisCommon/PMGTools/test/ut_PMGTruthWeightTool_test.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..a4bfec7a5ca437e7596170ea065b1ed51408de9e
--- /dev/null
+++ b/PhysicsAnalysis/AnalysisCommon/PMGTools/test/ut_PMGTruthWeightTool_test.cxx
@@ -0,0 +1,160 @@
+/*
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+/// Example standalone executable using TEvent (from POOL or xAODRootAccess) to read an xAOD
+/// Tests the PMGTruthWeightTool, showing example usage
+
+/// @author: James Robinson
+
+// EDM include(s):
+#include "AsgTools/MessageCheck.h"
+#include "AsgTools/AnaToolHandle.h"
+
+// Project dependent include(s)
+#ifdef XAOD_STANDALONE
+#include "xAODRootAccess/Init.h"
+#include "xAODRootAccess/TEvent.h"
+#else
+#include "POOLRootAccess/TEvent.h"
+#endif
+
+// Local include(s):
+#include "PMGTools/PMGTruthWeightTool.h"
+#include "PMGTools/PMGSherpaVjetsSysTool.h"
+
+// ROOT include(s):
+#include "TFile.h"
+#include <chrono>
+
+using namespace asg::msgUserCode; // messaging
+
+
+// Check for exceptions matching known text
+#define EXPECT_THROW(x,text) \
+{ std::string internalTestMessage; \
+ try { \
+ x; internalTestMessage = std::string ("expected statement ") + #x " to throw, but it didn't"; \
+ } catch (std::exception& e) { \
+ if (std::string(e.what()).find(std::string(text)) == std::string::npos) { \
+ internalTestMessage = std::string ("expected statement ") + #x " to throw message matching " + (text) + ", but actual message didn't match: " + e.what(); \
+ } \
+ } catch (...) { \
+ internalTestMessage = std::string ("statement ") + #x " threw an exception that didn't derive from std::exception"; \
+ } \
+ if (internalTestMessage.empty()) { ANA_MSG_INFO("Success!"); \
+ } else { ANA_MSG_FATAL(internalTestMessage); return -1; } \
+}
+
+// Check for functions returning false
+#define EXPECT_FALSE(x) \
+{ if (x == false) { ANA_MSG_INFO("Success!"); \
+ } else { ANA_MSG_FATAL("Expected " << #x << " to return false"); return -1; } \
+}
+
+
+int main(int argc, char *argv[])
+{
+ ANA_CHECK_SET_TYPE (int); // makes ANA_CHECK return ints if exiting function
+
+ // Initialise the application:
+#ifdef XAOD_STANDALONE
+ StatusCode::enableFailure();
+ ANA_CHECK(xAOD::Init());
+#else
+ IAppMgrUI *app = POOL::Init();
+#endif
+
+ // Use a default MC file for testing if none is provided
+ TString fileName = "$ASG_TEST_FILE_MC";
+ if (argc < 2) {
+ ANA_MSG_WARNING("No file name received, using ASG_TEST_FILE_MC");
+ } else {
+ fileName = argv[1]; // use the user provided file
+ }
+
+ // Initialise TEvent reading
+#ifdef XAOD_STANDALONE
+ std::unique_ptr<TFile> ptrFile(TFile::Open(fileName, "READ"));
+ xAOD::TEvent event(xAOD::TEvent::kClassAccess);
+ ANA_CHECK(event.readFrom(ptrFile.get()));
+#else
+ POOL::TEvent event(POOL::TEvent::kClassAccess);
+ ANA_CHECK(event.readFrom(fileName));
+#endif
+ ANA_MSG_INFO("Opened file: " << fileName);
+
+ // Hack to load the file
+ event.getEntries();
+
+ // Create the truth weight tool:
+ ANA_MSG_INFO("Creating PMGTruthWeightTool...");
+#ifdef XAOD_STANDALONE
+ asg::AnaToolHandle< PMGTools::IPMGTruthWeightTool > weightTool;
+ ASG_SET_ANA_TOOL_TYPE(weightTool, PMGTools::PMGTruthWeightTool);
+ weightTool.setName("PMGTruthWeightTool");
+ ANA_CHECK(weightTool.initialize());
+#else
+ asg::AnaToolHandle< PMGTools::IPMGTruthWeightTool > weightTool("PMGTools::PMGTruthWeightTool/PMGTruthWeightTool");
+ ANA_CHECK(weightTool.retrieve());
+#endif
+
+ // Loop over a few events:
+ ANA_MSG_INFO("Preparing to loop over events...");
+ const Long64_t nEntries = 5;
+ // double retrievalTimeNanoSeconds = 0;
+ for(Long64_t entry = 0; entry < nEntries; ++entry) {
+ if (event.getEntry(entry) < 0) { ANA_MSG_ERROR("Failed to read event " << entry); continue; }
+
+ auto weightNames = weightTool->getWeightNames();
+ ANA_MSG_INFO("Event #" << entry << ": found " << weightNames.size() << " weights for this event");
+
+ // Print out the first weight for every event
+ if (weightNames.size() > 1) {
+ static std::string weight_name_to_test = weightNames.at(1);
+ ANA_MSG_INFO("Event #" << entry << ": weight called '" << weight_name_to_test << "' = " << weightTool->getWeight(weight_name_to_test));
+ }
+
+ // Full print out for some of the events
+ if ((entry + 1) % 5 == 0) {
+ // Print out all weights and names
+ ANA_MSG_INFO("Printing all " << weightNames.size() << " weights for this event...");
+ unsigned int idx(0);
+ for (auto weight : weightNames) {
+ ANA_MSG_INFO("... weight " << idx++ << " has name \"" << weight << "\" and value " << weightTool->getWeight(weight));
+ }
+
+ // Give feedback about where we are
+ ANA_MSG_INFO("Processed " << entry + 1 << " / " << nEntries << " events");
+ }
+
+ // Check retrieval speed
+ ANA_MSG_INFO("Check average time taken to retrieve a weight");
+ if (entry == nEntries - 1) {
+ int nWeightCalls = 1000000;
+ static std::string weight_name_to_test = weightNames.back();
+ auto start = std::chrono::steady_clock::now();
+ for (int i = 0; i < nWeightCalls; ++i) {
+ weightTool->getWeight(weight_name_to_test);
+ }
+ auto elapsed = std::chrono::steady_clock::now() - start;
+ double retrievalTimeNanoSeconds = std::chrono::duration <double, std::nano> (elapsed).count() / nWeightCalls;
+ ANA_MSG_INFO("Retrieving " << nWeightCalls << " weights took " << retrievalTimeNanoSeconds << " ns per call");
+ }
+
+ // Check that retrieving a non-existent weight throws an exception
+ ANA_MSG_INFO("Check that retrieving a non-existent weight throws an exception");
+ EXPECT_THROW(weightTool->getWeight("non-existent-name"), "Weight \"non-existent-name\" could not be found");
+
+ // Check that asking for a non-existent weight returns false
+ ANA_MSG_INFO("Check that asking for a non-existent weight returns false");
+ EXPECT_FALSE(weightTool->hasWeight("non-existent-name"));
+ }
+
+ // trigger finalization of all services and tools created by the Gaudi Application
+#ifndef XAOD_STANDALONE
+ ANA_CHECK(app->finalize());
+#endif
+
+ return 0;
+}