diff --git a/Reconstruction/Jet/JetMomentTools/CMakeLists.txt b/Reconstruction/Jet/JetMomentTools/CMakeLists.txt index 96433ecda5cafeebe3b3fd5fbdfcc8f50b5814a9..8946417b2a89a3be4c38c9e9e874697a7e18091d 100644 --- a/Reconstruction/Jet/JetMomentTools/CMakeLists.txt +++ b/Reconstruction/Jet/JetMomentTools/CMakeLists.txt @@ -26,12 +26,14 @@ if( XAOD_STANDALONE OR XAOD_ANALYSIS ) Reconstruction/Jet/JetInterface Reconstruction/Jet/JetRec Reconstruction/Jet/JetUtils + Reconstruction/Jet/JetCalibTools Reconstruction/PFlow/PFlowUtils PRIVATE Calorimeter/CaloGeoHelpers Event/xAOD/xAODEventInfo Event/xAOD/xAODPFlow Tools/PathResolver + Event/xAOD/xAODMetaData ${extra_deps} ) else() atlas_depends_on_subdirs( @@ -48,6 +50,7 @@ else() Reconstruction/Jet/JetEDM Reconstruction/Jet/JetInterface Reconstruction/Jet/JetRec + Reconstruction/Jet/JetCalibTools Reconstruction/Jet/JetRecCalo Reconstruction/Jet/JetUtils PRIVATE @@ -59,12 +62,13 @@ else() Event/xAOD/xAODEventInfo Event/xAOD/xAODPFlow Tools/PathResolver + Event/xAOD/xAODMetaData GaudiKernel ) endif() # External dependencies: find_package( Boost ) -find_package( ROOT COMPONENTS Core Hist RIO ) +find_package( ROOT COMPONENTS Core Hist RIO TMVA) # Component(s) in the package: atlas_add_library( JetMomentToolsLib @@ -73,7 +77,8 @@ atlas_add_library( JetMomentToolsLib PRIVATE_INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${Boost_INCLUDE_DIRS} LINK_LIBRARIES ${ROOT_LIBRARIES} ${Boost_LIBRARIES} AsgTools xAODCaloEvent xAODJet xAODMissingET xAODTracking xAODTruth JetEDM JetInterface TrackVertexAssociationToolLib - JetRecLib JetUtils PFlowUtilsLib + JetRecLib JetUtils PFlowUtilsLib JetCalibToolsLib + xAODMetaData PRIVATE_LINK_LIBRARIES CaloGeoHelpers xAODEventInfo xAODPFlow PathResolver ) if( NOT XAOD_STANDALONE ) @@ -84,8 +89,8 @@ if( NOT XAOD_STANDALONE ) atlas_add_component( JetMomentTools src/*.h src/*.cxx src/components/*.cxx INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} - LINK_LIBRARIES ${ROOT_LIBRARIES} CaloIdentifier xAODCaloEvent xAODJet - GaudiKernel JetRecLib JetUtils PathResolver JetMomentToolsLib + LINK_LIBRARIES ${ROOT_LIBRARIES} CaloIdentifier xAODCaloEvent xAODJet xAODMetaData + GaudiKernel JetRecLib JetUtils PathResolver JetMomentToolsLib JetCalibToolsLib ${extra_libs} ) endif() diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetCaloEnergies.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetCaloEnergies.h index 3c3236f756fbf2c86fdf66b7203a189e56eafca9..dd89073614329666af3b2922d78be9e482f75211 100644 --- a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetCaloEnergies.h +++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetCaloEnergies.h @@ -1,29 +1,43 @@ // this file is -*- C++ -*- /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration */ #ifndef JETMOMENTTOOLS_JETCALOENERGIES_H #define JETMOMENTTOOLS_JETCALOENERGIES_H -#include "JetRec/JetModifierBase.h" +#include "AsgTools/AsgTool.h" +#include "JetInterface/IJetDecorator.h" +#include "StoreGate/ReadDecorHandleKey.h" +#include "StoreGate/ReadDecorHandle.h" +#include "StoreGate/WriteDecorHandleKey.h" +#include "StoreGate/WriteDecorHandle.h" +#include "xAODJet/JetContainer.h" #include <vector> -class JetCaloEnergies : public JetModifierBase { +class JetCaloEnergies : public asg::AsgTool, + virtual public IJetDecorator { ASG_TOOL_CLASS0(JetCaloEnergies) public: JetCaloEnergies(const std::string & t); - - virtual int modifyJet(xAOD::Jet& ) const ; + virtual StatusCode initialize() override; + virtual StatusCode decorate(const xAOD::JetContainer& jets) const override; protected: - void fillEperSamplingCluster(xAOD::Jet &jet, std::vector<float> & ePerSampling ) const ; - void fillEperSamplingPFO(xAOD::Jet &jet, std::vector<float> & ePerSampling ) const ; + void fillEperSamplingCluster(const xAOD::Jet &jet, std::vector<float> & ePerSampling ) const ; + void fillEperSamplingPFO(const xAOD::Jet &jet, std::vector<float> & ePerSampling ) const ; + +private: + Gaudi::Property<std::string> m_jetContainerName{this, "JetContainer", "", "SG key for the input jet container"}; - + SG::WriteDecorHandleKey<xAOD::JetContainer> m_ePerSamplingKey{this, "EPerSamplingName", "EnergyPerSampling", "SG key for the EnergyPerSampling attribute"}; + SG::WriteDecorHandleKey<xAOD::JetContainer> m_emFracKey{this, "EMFracName", "EMFrac", "SG key for the EMFrac attribute"}; + SG::WriteDecorHandleKey<xAOD::JetContainer> m_hecFracKey{this, "HECFracName", "HECFrac", "SG key for the HECFrac attribute"}; + SG::WriteDecorHandleKey<xAOD::JetContainer> m_psFracKey{this, "PSFracName", "PSFrac", "SG key for the PSFrac attribute"}; + }; diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetForwardJvtTool.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetForwardJvtTool.h index 0af9b3aebe1c2a1bbcbcede01ebdc911f9015a2a..128e975e1737383f143e8550728e83e4b8c32049 100644 --- a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetForwardJvtTool.h +++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetForwardJvtTool.h @@ -1,7 +1,7 @@ ///////////////////////// -*- C++ -*- ///////////////////////////// /* - Copyright (C) 2002-2017, 2019 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration */ // JetForwardJvtTool.h @@ -18,20 +18,24 @@ // FrameWork includes #include "AsgTools/ToolHandle.h" #include "AsgTools/AsgTool.h" +#include "StoreGate/ReadDecorHandleKey.h" +#include "StoreGate/ReadDecorHandle.h" +#include "StoreGate/WriteDecorHandleKey.h" +#include "StoreGate/WriteDecorHandle.h" // EDM includes #include "xAODJet/JetContainer.h" #include "xAODMissingET/MissingETContainer.h" #include "xAODCaloEvent/CaloCluster.h" -#include "JetInterface/IJetModifier.h" +#include "JetInterface/IJetDecorator.h" #include "AsgTools/IAsgTool.h" //ASG_TOOL_INTERFACE(IJetUpdateJvt) class JetForwardJvtTool : public asg::AsgTool, - virtual public IJetModifier{ - ASG_TOOL_CLASS(JetForwardJvtTool,IJetModifier) + virtual public IJetDecorator{ + ASG_TOOL_CLASS(JetForwardJvtTool,IJetDecorator) // This macro defines the constructor with the interface declaration //ASG_TOOL_CLASS(JetForwardJvtTool, IJetForwardJvtTool) @@ -50,10 +54,10 @@ virtual ~JetForwardJvtTool(); // Athena algtool's Hooks - StatusCode initialize(); - StatusCode finalize(); + virtual StatusCode initialize() override; + virtual StatusCode finalize() override; - virtual StatusCode modify(xAOD::JetContainer& jetCont) const; + virtual StatusCode decorate(const xAOD::JetContainer& jetCont) const override; float getFJVT(const xAOD::Jet *jet) const; bool forwardJet(const xAOD::Jet *jet) const; @@ -61,37 +65,43 @@ float getDrpt(const xAOD::Jet *jet) const; int getJetVertex(const xAOD::Jet *jet) const; - static StatusCode tagTruth(const xAOD::JetContainer *jets,const xAOD::JetContainer *truthJets); + StatusCode tagTruth(const xAOD::JetContainer *jets,const xAOD::JetContainer *truthJets); void calculateVertexMomenta(const xAOD::JetContainer *jets) const; float getCombinedWidth(const xAOD::Jet *jet) const; private: - std::string m_orLabel; - std::string m_outLabel; - double m_etaThresh; - double m_forwardMinPt; - double m_forwardMaxPt; - double m_centerMinPt; - double m_centerMaxPt; - double m_centerJvtThresh; - std::string m_jvtMomentName; - double m_centerDrptThresh; - double m_maxStochPt; - double m_jetScaleFactor; - double m_fjvtThresh; - bool m_tightOP; + Gaudi::Property<double> m_etaThresh{this, "EtaThresh", 2.5, "Eta threshold"}; + Gaudi::Property<double> m_forwardMinPt{this, "ForwardMinPt", 20e3, "Forward minimum Pt"}; + Gaudi::Property<double> m_forwardMaxPt{this, "ForwardMaxPt", 50e3, "Forward maximum Pt"}; + Gaudi::Property<double> m_centerMinPt{this, "CentralMinPt", 20e3, "Central minimum Pt"}; + Gaudi::Property<double> m_centerMaxPt{this, "CentralMaxPt", -1, "Central maximum Pt"}; + Gaudi::Property<double> m_centerJvtThresh{this, "CentralJvtThresh", 0.14, "Central JVT threshold"}; + Gaudi::Property<double> m_centerDrptThresh{this, "CentralDrptThresh", 0.2, "Central drpt threshold"}; + Gaudi::Property<double> m_maxStochPt{this, "CentralMaxStochPt", 35e3, "Central maximum StochPt"}; + Gaudi::Property<double> m_jetScaleFactor{this, "JetScaleFactor", 0.4, "Jet scale factor"}; + Gaudi::Property<double> m_fjvtThresh{this, "FjvtThresh", 15e3, "FJVT threshold"}; //15GeV->92%,11GeV->85% + Gaudi::Property<bool> m_tightOP{this, "UseTightOP", false, "Use tight (true) or loose (false)"}; mutable std::vector<TVector2> m_pileupMomenta; mutable size_t m_pvind; - SG::AuxElement::Decorator<char>* m_Dec_OR = NULL; - SG::AuxElement::Decorator<char>* m_Dec_out = NULL; void getPV() const; /// Default constructor: JetForwardJvtTool(); - SG::ReadHandleKey<xAOD::VertexContainer> m_vertexContainer_key; - SG::ReadHandleKey<xAOD::MissingETContainer> m_trkMET_key; + Gaudi::Property<std::string> m_jetContainerName{this, "JetContainer", "", "SG key for the input jet container"}; + SG::ReadDecorHandleKey<xAOD::JetContainer> m_orKey{this, "OverlapDec", "", "Overlap decoration key"}; + SG::WriteDecorHandleKey<xAOD::JetContainer> m_outKey{this, "OutputDec", "passFJVT", "Output decoration key"}; + SG::WriteDecorHandleKey<xAOD::JetContainer> m_isHSKey{this, "IsJvtHSName", "isJvtHS", "Decoration key for isJvtHS"}; + SG::WriteDecorHandleKey<xAOD::JetContainer> m_isPUKey{this, "IsJvtPUName", "isJvtPU", "Decoration key for isJvtPU"}; + SG::WriteDecorHandleKey<xAOD::JetContainer> m_fjvtDecKey{this, "FJVTName", "fJvt", "Decoration key for fJvt"}; + + SG::ReadHandleKey<xAOD::VertexContainer> m_vertexContainerName{this, "VertexContainerName", "PrimaryVertices", "SG key for vertex container"}; + SG::ReadHandleKey<xAOD::MissingETContainer> m_trkMETName{this, "Met_TrackName", "Met_Track", "SG key for MET track container"}; + + SG::ReadDecorHandleKey<xAOD::JetContainer> m_widthKey{this, "WidthName", "Width", "SG key for jet width"}; + SG::ReadDecorHandleKey<xAOD::JetContainer> m_jvtMomentKey{this, "JvtMomentName", "Jvt", "JVT moment name"}; + SG::ReadDecorHandleKey<xAOD::JetContainer> m_sumPtsKey{this, "SumPtsName", "SumPtTrkPt500", "SG key for SumPt vector"}; }; #endif //> !FORWARDJVTTOOL_JVT_FORWARDJVTTOOL_H diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetForwardJvtToolBDT.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetForwardJvtToolBDT.h new file mode 100644 index 0000000000000000000000000000000000000000..ec341e1086c634cd3a00c4c421ede91feddad6b4 --- /dev/null +++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetForwardJvtToolBDT.h @@ -0,0 +1,151 @@ +///////////////////////// -*- C++ -*- ///////////////////////////// + +/* + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration +*/ + +// JetForwardJvtToolBDT.h +// Header file for class JetForwardJvtToolBDT +// Author: Louis Portales<louis.portales@cern.ch> +/////////////////////////////////////////////////////////////////// +#ifndef FORWARDJVTTOOLBDT_JVT_FORWARDJVTTOOLBDT_H +#define FORWARDJVTTOOLBDT_JVT_FORWARDJVTTOOLBDT_H 1 + +// STL includes +#include <string> + +// Root includes +#include "TMVA/Reader.h" +#include <TH3.h> +#include <TString.h> +#include <TFile.h> + +// FrameWork includes +#include "AsgTools/ToolHandle.h" +#include "AsgTools/AnaToolHandle.h" + +// EDM includes +#include "xAODEventInfo/EventInfo.h" +#include "xAODJet/JetContainer.h" +#include "xAODMissingET/MissingETContainer.h" +#include "xAODCaloEvent/CaloCluster.h" +#include "xAODCaloEvent/CaloClusterContainer.h" +#include "JetInterface/IJetDecorator.h" +#include "AsgTools/IAsgTool.h" +#include "StoreGate/ReadDecorHandleKey.h" +#include "StoreGate/ReadDecorHandle.h" +#include "StoreGate/WriteDecorHandleKey.h" +#include "StoreGate/WriteDecorHandle.h" + + +namespace TMVA{ class Reader; } + +class JetForwardJvtToolBDT + : public asg::AsgTool, + virtual public IJetDecorator{ + ASG_TOOL_CLASS(JetForwardJvtToolBDT,IJetDecorator) + + /////////////////////////////////////////////////////////////////// + // Public methods: + /////////////////////////////////////////////////////////////////// + public: + + + /// Constructor with parameters: + JetForwardJvtToolBDT(const std::string& name); + + /// Destructor: + virtual ~JetForwardJvtToolBDT(); + + // Athena algtool's Hooks + virtual StatusCode initialize() override; + + virtual StatusCode decorate(const xAOD::JetContainer& jetCont) const override; + + // MVfJVT + StatusCode getInputs ( const xAOD::Jet *jet ) const; + float getMVfJVT ( const xAOD::Jet *jet, int pvind , std::vector<TVector2> pileupMomenta ) const; + float getFJVT ( const xAOD::Jet *jet, int pvind , std::vector<TVector2> pileupMomenta ) const; + bool forwardJet ( const xAOD::Jet *jet ) const; + bool centralJet ( const xAOD::Jet *jet ) const; + float getDrpt ( const xAOD::Jet *jet ) const; + int getJetVertex ( const xAOD::Jet *jet ) const; + bool passMVfJVT ( float mvfjvt , float pt, float eta ) const; + + StatusCode tagTruth(const xAOD::JetContainer *jets,const xAOD::JetContainer *truthJets); + std::vector<TVector2> calculateVertexMomenta(const xAOD::JetContainer *jets, int pvind) const; + int getPV() const; + +private: + + Gaudi::Property<std::string> m_configDir{this, "configDir", "JetPileupTag/MVfJVT/", "Configuration directory"}; + Gaudi::Property<std::vector<std::string> > m_MVconfig{this, "ConfigFiles", { + "weights/MVfJVT_pt2030_etaHigh_muHigh.May2019.weights.xml", + "weights/MVfJVT_pt2030_etaLow_muHigh.May2019.weights.xml", + "weights/MVfJVT_pt3040_etaHigh_muHigh.May2019.weights.xml", + "weights/MVfJVT_pt3040_etaLow_muHigh.May2019.weights.xml", + "weights/MVfJVT_pt4050_etaHigh_muHigh.May2019.weights.xml", + "weights/MVfJVT_pt4050_etaLow_muHigh.May2019.weights.xml", + "weights/MVfJVT_pt50plus_etaHigh_muHigh.May2019.weights.xml", + "weights/MVfJVT_pt50plus_etaLow_muHigh.May2019.weights.xml", + "weights/MVfJVT_pt2030_etaHigh_muLow.May2019.weights.xml", + "weights/MVfJVT_pt2030_etaLow_muLow.May2019.weights.xml", + "weights/MVfJVT_pt3040_etaHigh_muLow.May2019.weights.xml", + "weights/MVfJVT_pt3040_etaLow_muLow.May2019.weights.xml", + "weights/MVfJVT_pt4050_etaHigh_muLow.May2019.weights.xml", + "weights/MVfJVT_pt4050_etaLow_muLow.May2019.weights.xml", + "weights/MVfJVT_pt50plus_etaHigh_muLow.May2019.weights.xml", + "weights/MVfJVT_pt50plus_etaLow_muLow.May2019.weights.xml" + }, "List of config file names"}; // pt [20,30,40,50,120] || |eta| [2.5,3.2,4.5] || mu [0,50,inf.]; + Gaudi::Property<std::string> m_wpFile{this, "WPfile", "MVfJVT_WPs.Nov2019.root", "WP file"}; + + std::unique_ptr< TFile > m_wpFileIn; + std::unique_ptr< TMVA::Reader > m_MVreader; + std::unique_ptr< TH3D > m_mvfjvtThresh; + + + Gaudi::Property<std::string> m_orLabel{this, "OverlapDec", "", "SG key for the overlap decoration"}; + + Gaudi::Property<double> m_etaThresh{this, "EtaThresh", 2.5, "Eta threshold"}; + Gaudi::Property<double> m_forwardMinPt{this, "ForwardMinPt", 20e3, "Forward minimum pt"}; + Gaudi::Property<double> m_forwardMaxPt{this, "ForwardMaxPt", 120e3, "Forward maximum pt"}; + Gaudi::Property<double> m_centerMinPt{this, "CentralMinPt", 20e3, "Central minimum pt"}; + Gaudi::Property<double> m_centerMaxPt{this, "CentralMaxPt", -1, "Central maximum pt (set to -1 for no limit)"}; + Gaudi::Property<double> m_centerJvtThresh{this, "CentralJvtThresh", 0.11, "Central JVT threshold"}; + Gaudi::Property<std::string> m_jvtMomentName{this, "JvtMomentName", "Jvt", "SG key for JVT moment"}; + Gaudi::Property<double> m_centerDrptThresh{this, "CentralDrptThresh", 0.2, "Central drpt threshold"}; + Gaudi::Property<double> m_maxStochPt{this, "CentralMaxStochPt", 35e3, "Central maximum stochpt"}; + Gaudi::Property<double> m_jetScaleFactor{this, "JetScaleFactor", 0.4, "Jet scale factor"}; + Gaudi::Property<std::string> m_OP{this, "OperatingPoint", "DEFAULT", "Selected operating point, can be 'LOOSE', 'TIGHT' or 'TIGHTER'"}; + Gaudi::Property<bool> m_getTagger{this, "retrieveTagger", false, "Whether to retrieve the tagger"}; + Gaudi::Property<bool> m_isAna{this, "AnaToolMode", false, "True if running in AnaTool mode"}; + Gaudi::Property<int> m_pvind{this, "PVIndexHS", -1, ""}; + + Gaudi::Property<std::string> m_jetContainerName{this, "JetContainer", "", "SG key of input jet container"}; + + SG::WriteDecorHandleKey<xAOD::JetContainer> m_outMVKey{this, "OutputDecMV", "passMVFJVT", "SG key for the output MV decoration"}; + + SG::WriteDecorHandleKey<xAOD::JetContainer> m_mvfjvtKey{this, "MVFJVTName", "MVfJVT", "SG key for the output MVfJVT decoration"}; + SG::WriteDecorHandleKey<xAOD::JetContainer> m_cllambda2Key{this, "cllambda2Name", "LeadclSecondLambda", "SG key for the LeadclSecondLambda decoration"}; + SG::WriteDecorHandleKey<xAOD::JetContainer> m_clwidthKey{this, "cletawidthName", "LeadclWidth", "SG key for the cluster eta width decoration"}; + SG::WriteDecorHandleKey<xAOD::JetContainer> m_clisoKey{this, "clisoName", "SumclIso", "SG key for the cluster isolation decoration"}; + SG::WriteDecorHandleKey<xAOD::JetContainer> m_clemprobKey{this, "clemprobName", "SumclEMprob", "SG key for the cluster EMprob decoration"}; + SG::WriteDecorHandleKey<xAOD::JetContainer> m_cleKey{this, "cleName", "Sumcle", "SG key for the cluster energy decoration"}; + + SG::ReadDecorHandleKey<xAOD::JetContainer> m_lcllambda2NTKey{this, "lcllambda2NTName", "LeadingClusterSecondLambda", "Leading cluster second lambda to use if getTagger is false"}; + SG::ReadDecorHandleKey<xAOD::JetContainer> m_lcllambda2Key{this, "lcllambda2Name", "DFCommonJets_MVfJVT_LeadclSecondLambda", "Leading cluster second lambda to use if getTagger is true"}; + SG::ReadDecorHandleKey<xAOD::JetContainer> m_lclwidthKey{this, "lclwidthName", "DFCommonJets_MVfJVT_LeadclWidth", "Leading cluster width to use if getTagger is true"}; + SG::ReadDecorHandleKey<xAOD::JetContainer> m_lclisoKey{this, "lclisoName", "DFCommonJets_MVfJVT_SumclIso", "Leading cluster isolation to use if getTagger is true"}; + SG::ReadDecorHandleKey<xAOD::JetContainer> m_lclemprobKey{this, "lclemprobName", "DFCommonJets_MVfJVT_SumclEMprob", "Leading cluster EMprob to use if getTagger is true"}; + SG::ReadDecorHandleKey<xAOD::JetContainer> m_lcleKey{this, "lcleName", "DFCommonJets_MVfJVT_Sumcle", "Leading cluster energy to use if getTagger is true"}; + + SG::ReadHandleKey<xAOD::EventInfo> m_eventInfoKey{this, "EventInfoName", "EventInfo", "SG key for input EventInfo"}; + SG::ReadHandleKey<xAOD::VertexContainer> m_vertexContainerKey{this, "VertexContainerName", "PrimaryVertices", "SG key for input vertex container"}; + SG::ReadHandleKey<xAOD::CaloClusterContainer> m_caloClusterContainerKey{this, "CaloClusterContainerName" "CaloCalTopoClusters", "SG key for input calo cluster container"}; + SG::ReadHandleKey<xAOD::MissingETContainer> m_trkMetKey{this, "TrackMetName", "MET_Track", "SG key for input track MET container"}; + + SG::WriteDecorHandleKey<xAOD::JetContainer> m_isHSKey{this, "isHSName", "isJvtHS", "SG key for output isHS decoration"}; + SG::WriteDecorHandleKey<xAOD::JetContainer> m_isPUKey{this, "isPUName", "isJvtPU", "SG key for output isPU decoration"}; + +}; +#endif //> !FORWARDJVTTOOLBDT_JVT_FORWARDJVTTOOLBDT_H diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetIsolationTool.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetIsolationTool.h index c3855d1fc3a391f879c37d05504cb05a545adb1d..3d6cd25d95c78a7b5a7aba6a405c6423814afee2 100644 --- a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetIsolationTool.h +++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetIsolationTool.h @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration */ // JetIsolationTool.h @@ -71,11 +71,17 @@ /////////////////////////////////////////////////////////////////////// #include "AsgTools/AsgTool.h" -#include "JetRec/JetModifierBase.h" #include "JetUtils/TiledEtaPhiMap.h" +#include "JetInterface/IJetDecorator.h" #include "fastjet/PseudoJet.hh" #include "JetRec/PseudoJetContainer.h" #include "JetEDM/IConstituentUserInfo.h" +#include "StoreGate/ReadDecorHandleKey.h" +#include "StoreGate/ReadDecorHandle.h" +#include "StoreGate/WriteDecorHandleKey.h" +#include "StoreGate/WriteDecorHandleKeyArray.h" +#include "StoreGate/WriteDecorHandle.h" +#include "xAODJet/JetContainer.h" namespace jet { @@ -120,7 +126,8 @@ protected: //********************************************************************** -class JetIsolationTool : public JetModifierBase { +class JetIsolationTool : public asg::AsgTool, + virtual public IJetDecorator { ASG_TOOL_CLASS0(JetIsolationTool) public: @@ -134,17 +141,22 @@ public: ~JetIsolationTool(); // Athena algtool Hooks - StatusCode initialize(); + virtual StatusCode initialize() override; StatusCode finalize(); // Jet Modifier methods. - StatusCode modify(xAOD::JetContainer& jets) const; - int modifyJet(xAOD::Jet& ) const{return 0;}; + StatusCode decorate(const xAOD::JetContainer& jets) const override; private: + Gaudi::Property<std::vector<std::string>> m_isolationCodes{this, "IsolationCalculations", {}, "Isolation calculation data vector"}; + Gaudi::Property<std::string> m_jetContainerName{this, "JetContainer", "", "SG key for the input jet container"}; - std::vector<std::string> m_isolationCodes; SG::ReadHandleKey<PseudoJetContainer> m_pjsin{this, "PseudoJetsIn", "", "PseudoJetContainer to read"}; + SG::ReadDecorHandleKey<xAOD::JetContainer> m_inputTypeKey{this, "InputTypeName", "InputType", "Key for the InputType field of a jet"}; + SG::WriteDecorHandleKeyArray<xAOD::JetContainer> m_perpKeys{this, "PerpName", {}, "SG key for output perpendicular momentum component decoration (not to be configured manually!)"}; + SG::WriteDecorHandleKeyArray<xAOD::JetContainer> m_sumPtKeys{this, "SumPtName", {}, "SG key for output SumPt decoration (not to be configured manually!)"}; + SG::WriteDecorHandleKeyArray<xAOD::JetContainer> m_parKeys{this, "ParName", {}, "SG key for output parallel momentum component decoration (not to be configured manually!)"}; + SG::WriteDecorHandleKeyArray<xAOD::JetContainer> m_pKeys{this, "PName", {}, "SG key for output momentum decoration (not to be configured manually!)"}; /// the list of isolation calculation objects (they are actually used /// only as template objects from which the actual calculators are build diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetMuonSegmentMomentsTool.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetMuonSegmentMomentsTool.h deleted file mode 100644 index 1f7d1e39c14ead6730a1b547525d804c1cfc2ed3..0000000000000000000000000000000000000000 --- a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetMuonSegmentMomentsTool.h +++ /dev/null @@ -1,38 +0,0 @@ -/* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -*/ - -// JetMuonSegmentMomentsTool.h - -#ifndef JETMOMENTTOOLS_JETMUONSEGMENTMOMENTSTOOL_H -#define JETMOMENTTOOLS_JETMUONSEGMENTMOMENTSTOOL_H - -/// Steven Schramm \n -/// February 2014 -/// -/// Tool to calculate general muon segment based jet moments -/// Primarily intended for punch-through studies and corrections -/// -/// Moments to calculate: -/// None yet... - -#include "JetRec/JetModifierBase.h" - -class JetMuonSegmentMomentsTool : public JetModifierBase -{ - ASG_TOOL_CLASS(JetMuonSegmentMomentsTool,IJetModifier) - - public: - // Constructor from tool name - JetMuonSegmentMomentsTool(const std::string& name); - - // Inherited methods to modify a jet - virtual int modifyJet(xAOD::Jet& jet) const; - - private: - // Configurable parameters - std::string m_assocMuonSegName; -}; - -#endif - diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetOriginCorrectionTool.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetOriginCorrectionTool.h index 7432d9d87fe848df011c07cbec1ca587588ed010..3b3cf4f678b946988120f0e4ca9eaa35f68c1752 100644 --- a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetOriginCorrectionTool.h +++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetOriginCorrectionTool.h @@ -1,7 +1,7 @@ // JetOriginCorrectionTool.h -*- C++ -*- /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration */ #ifndef JetMomentTools_JetOriginCorrectionTool_H @@ -24,12 +24,15 @@ //////////////////////////////////////////////////////////// #include "AsgTools/AsgTool.h" -#include "JetInterface/IJetModifier.h" +#include "JetInterface/IJetDecorator.h" +#include "StoreGate/WriteDecorHandleKey.h" +#include "StoreGate/WriteDecorHandle.h" #include "xAODEventInfo/EventInfo.h" +#include "xAODJet/JetContainer.h" class JetOriginCorrectionTool : public asg::AsgTool, - virtual public IJetModifier { - ASG_TOOL_CLASS(JetOriginCorrectionTool, IJetModifier) + virtual public IJetDecorator { + ASG_TOOL_CLASS(JetOriginCorrectionTool, IJetDecorator) public: @@ -38,19 +41,26 @@ public: /// Inherited method to modify a jet container. Compute the origin-corrected /// momentum and put it in the jets - StatusCode modify(xAOD::JetContainer& jet) const override; - StatusCode initialize() override; + StatusCode decorate(const xAOD::JetContainer& jet) const override; + virtual StatusCode initialize() override; protected: - std::string m_correctionName; - bool m_onlyAssignPV; + Gaudi::Property<std::string> m_correctionName{this, "OriginCorrectedName", "JetOriginConstitScaleMomentum", "Origin corrected name"}; + Gaudi::Property<std::string> m_jetContainerName{this, "JetContainer", "", "SG key for the input jet container"}; + Gaudi::Property<bool> m_onlyAssignPV{this, "OnlyAssignPV", false, "Only write out PV information"}; private: - SG::ReadHandleKey< xAOD::VertexContainer> m_vertexContainer_key; - SG::ReadHandleKey<xAOD::EventInfo> m_eventInfo_key; + SG::ReadHandleKey<xAOD::VertexContainer> m_vertexContainer_key{this, "VertexContainer", "PrimaryVertices", "Primary vertex container name"}; + SG::ReadHandleKey<xAOD::EventInfo> m_eventInfo_key{this, "EventInfoName", "EventInfo", "Event info name"}; + Gaudi::Property<std::string> m_scaleMomentumName{this, "jetScaleMomentName", "JetOriginConstitScaleMomentum", "SG key for JetScaleMomentum components"}; + SG::WriteDecorHandleKey<xAOD::JetContainer> m_scaleMomentumPtKey{this, "JetScaleMomentumPtName", "pt", "SG suffix for output JetScaleMomentum (pt) decorator"}; + SG::WriteDecorHandleKey<xAOD::JetContainer> m_scaleMomentumPhiKey{this, "JetScaleMomentumPhiName", "phi", "SG suffix for output JetScaleMomentum (phi) decorator"}; + SG::WriteDecorHandleKey<xAOD::JetContainer> m_scaleMomentumEtaKey{this, "JetScaleMomentumEtaName", "eta", "SG suffix for output JetScaleMomentum (eta) decorator"}; + SG::WriteDecorHandleKey<xAOD::JetContainer> m_scaleMomentumMKey{this, "JetScaleMomentumMName", "m", "SG suffix for output JetScaleMomentum (mass) decorator"}; + SG::WriteDecorHandleKey<xAOD::JetContainer> m_originVertexKey{this, "OriginVertexName", "OriginVertex", "SG key for output OriginVertex decorator"}; }; #endif diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetPtAssociationTool.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetPtAssociationTool.h index c4fe99d58f6d9f3f689a35f0458d9d50fb574904..431f27d0a9cba23735483efad39c4c0a4c5aacc2 100644 --- a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetPtAssociationTool.h +++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetPtAssociationTool.h @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration */ // JetPtAssociationTool.h @@ -16,11 +16,17 @@ /// AssociationName - Name of the association /// InputContainer - Name of the matching jet container -#include "JetRec/JetModifierBase.h" +#include "AsgTools/AsgTool.h" +#include "JetInterface/IJetDecorator.h" +#include "StoreGate/ReadHandleKey.h" +#include "StoreGate/ReadHandle.h" +#include "StoreGate/WriteDecorHandleKey.h" +#include "StoreGate/WriteDecorHandle.h" #include "xAODJet/JetContainer.h" -class JetPtAssociationTool : public JetModifierBase { - ASG_TOOL_CLASS(JetPtAssociationTool, IJetModifier) +class JetPtAssociationTool : public asg::AsgTool, + virtual public IJetDecorator { + ASG_TOOL_CLASS(JetPtAssociationTool, IJetDecorator) public: @@ -34,12 +40,12 @@ public: JetPtAssociationTool(std::string myname); /// Initialization. - StatusCode initialize(); + virtual StatusCode initialize() override; /// Inherited method to modify a jet. /// Extract and record the element link and pT fraction /// for the jet with the highest such fraction. - int modifyJet(xAOD::Jet& jet) const; + virtual StatusCode decorate(const xAOD::JetContainer& jets) const override; private: // data @@ -59,8 +65,11 @@ private: // data int match(const APVector& aps, const xAOD::Jet& jet, APVector& apvs, double& ptsum_constituents) const; /// Properties. - std::string m_aname; - SG::ReadHandleKey<xAOD::JetContainer> m_jetContainer_key; + Gaudi::Property<std::string> m_aname{this, "AssociationName", "", "Key for association vector"}; + Gaudi::Property<std::string> m_jetContainerName{this, "JetContainer", "", "Input jet container"}; + + SG::WriteDecorHandleKey<xAOD::JetContainer> m_assocFracKey{this, "AssociationFractionName", "AssociationFraction", "SG key for output AssociationFraction decoration"}; + SG::WriteDecorHandleKey<xAOD::JetContainer> m_assocLinkKey{this, "AssociationLinkName", "AssociationLink", "SG key for output AssociationLink decoration"}; }; #endif diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetTrackSumMomentsTool.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetTrackSumMomentsTool.h index 729c7b63ed9db0c942829c2a4ce599f614965de1..bd9847b459e7554711a6aecca07be5df1b51cd8f 100644 --- a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetTrackSumMomentsTool.h +++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetTrackSumMomentsTool.h @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration */ // JetTrackSumMomentsTool.h @@ -22,10 +22,15 @@ #include "AsgTools/ToolHandle.h" #include "AsgTools/AsgTool.h" +#include "JetInterface/IJetDecorator.h" #include "JetInterface/IJetTrackSelector.h" -#include "JetRec/JetModifierBase.h" #include "JetEDM/TrackVertexAssociation.h" +#include "StoreGate/ReadDecorHandleKey.h" +#include "StoreGate/ReadDecorHandle.h" +#include "StoreGate/WriteDecorHandleKey.h" +#include "StoreGate/WriteDecorHandle.h" +#include "xAODJet/JetContainer.h" #include "xAODTracking/TrackParticle.h" #include "xAODTracking/TrackParticleContainer.h" #include "xAODTracking/Vertex.h" @@ -35,8 +40,9 @@ #include <string> -class JetTrackSumMomentsTool : public JetModifierBase { - ASG_TOOL_CLASS(JetTrackSumMomentsTool,IJetModifier) +class JetTrackSumMomentsTool : public asg::AsgTool, + virtual public IJetDecorator { + ASG_TOOL_CLASS(JetTrackSumMomentsTool,IJetDecorator) public: @@ -44,11 +50,11 @@ public: JetTrackSumMomentsTool(const std::string& name); // Initialization. - StatusCode initialize(); + virtual StatusCode initialize() override; // Inherited methods to modify a jet // Calls moment and puts the results in the jet - virtual int modifyJet(xAOD::Jet& jet) const; + virtual StatusCode decorate(const xAOD::JetContainer& jets) const override; // Local method to return the vector track sums std::pair<float,float> @@ -61,13 +67,17 @@ public: private: - std::string m_assocTracksName; - bool m_requireTrackPV; - ToolHandle<IJetTrackSelector> m_htsel; + Gaudi::Property<std::string> m_assocTracksName{this, "AssociatedTracks", "", "SG key for associated tracks container"}; + Gaudi::Property<bool> m_requireTrackPV{this, "RequireTrackPV", true, "Require track to be from the primary vertex?"}; + Gaudi::Property<std::string> m_jetContainerName{this, "JetContainer", "", "SG key for the input jet container"}; + ToolHandle<IJetTrackSelector> m_htsel{this, "TrackSelector", "", "Track selector"}; - SG::ReadHandleKey< xAOD::VertexContainer> m_vertexContainer_key; - SG::ReadHandleKey<jet::TrackVertexAssociation> m_tva_key; + SG::ReadHandleKey< xAOD::VertexContainer> m_vertexContainer_key{this, "VertexContainer", "", "Vertex container key"}; + SG::ReadHandleKey<jet::TrackVertexAssociation> m_tva_key{this, "TrackVertexAssociation", "", "Track vertex association key"}; + + SG::WriteDecorHandleKey<xAOD::JetContainer> m_trackSumPtKey{this, "SumPtName", "SumPt", "SG key for output track SumPt decoration"}; + SG::WriteDecorHandleKey<xAOD::JetContainer> m_trackSumMassKey{this, "SumMassName", "SumMass", "SG key for output track SumMass decoration"}; }; #endif diff --git a/Reconstruction/Jet/JetMomentTools/Root/JetCaloEnergies.cxx b/Reconstruction/Jet/JetMomentTools/Root/JetCaloEnergies.cxx index 94d6a21512c9de5514bae928c3f5cbbcf1486fd6..4a6ae0627097623de345051add25e23983c7c95e 100644 --- a/Reconstruction/Jet/JetMomentTools/Root/JetCaloEnergies.cxx +++ b/Reconstruction/Jet/JetMomentTools/Root/JetCaloEnergies.cxx @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration */ #include "JetMomentTools/JetCaloEnergies.h" @@ -16,43 +16,66 @@ //********************************************************************** JetCaloEnergies::JetCaloEnergies(const std::string& name) -: JetModifierBase(name) { } +: AsgTool(name) { } //********************************************************************** -int JetCaloEnergies::modifyJet(xAOD::Jet& jet) const { - ATH_MSG_VERBOSE("Begin modifying jet."); - static const xAOD::JetAttributeAccessor::AccessorWrapper< std::vector<float> >& ePerSamplingAcc = - *xAOD::JetAttributeAccessor::accessor< std::vector<float> >(xAOD::JetAttribute::EnergyPerSampling); - size_t numConstit = jet.numConstituents(); - std::vector<float> & ePerSampling = ePerSamplingAcc(jet); - ePerSampling.resize(CaloSampling::Unknown, 0.); - for ( float& e : ePerSampling ) e = 0.0; // re-initialize - - if ( numConstit == 0 ) { - ATH_MSG_VERBOSE("Jet has no constituents."); - return 1; +StatusCode JetCaloEnergies::initialize() { + ATH_MSG_INFO("Initializing JetCaloEnergies " << name()); + + if(m_jetContainerName.empty()){ + ATH_MSG_ERROR("JetCaloEnergies needs to have its input jet container configured!"); + return StatusCode::FAILURE; } + + m_ePerSamplingKey = m_jetContainerName + "." + m_ePerSamplingKey.key(); + m_emFracKey = m_jetContainerName + "." + m_emFracKey.key(); + m_hecFracKey = m_jetContainerName + "." + m_hecFracKey.key(); + m_psFracKey = m_jetContainerName + "." + m_psFracKey.key(); + + ATH_CHECK(m_ePerSamplingKey.initialize()); + ATH_CHECK(m_emFracKey.initialize()); + ATH_CHECK(m_hecFracKey.initialize()); + ATH_CHECK(m_psFracKey.initialize()); + + return StatusCode::SUCCESS; +} + +//********************************************************************** - // should find a more robust solution than using 1st constit type. - xAOD::Type::ObjectType ctype = jet.rawConstituent( 0 )->type(); - if ( ctype == xAOD::Type::CaloCluster ) { - ATH_MSG_VERBOSE(" Constituents are calo clusters."); - fillEperSamplingCluster(jet, ePerSampling); +StatusCode JetCaloEnergies::decorate(const xAOD::JetContainer& jets) const { + ATH_MSG_VERBOSE("Begin decorating jets."); + for(const xAOD::Jet* jet : jets) { + SG::WriteDecorHandle<xAOD::JetContainer, std::vector<float> > ePerSamplingHandle(m_ePerSamplingKey); + size_t numConstit = jet->numConstituents(); + ePerSamplingHandle(*jet) = std::vector<float>(CaloSampling::Unknown, 0.); + std::vector<float>& ePerSampling = ePerSamplingHandle(*jet); + for ( float& e : ePerSampling ) e = 0.0; // re-initialize + + if ( numConstit == 0 ) { + ATH_MSG_VERBOSE("Jet has no constituents."); + continue; + } - } else if (ctype == xAOD::Type::ParticleFlow) { - ATH_MSG_VERBOSE(" Constituents are pflow objects."); - fillEperSamplingPFO(jet, ePerSampling); + // should find a more robust solution than using 1st constit type. + xAOD::Type::ObjectType ctype = jet->rawConstituent( 0 )->type(); + if ( ctype == xAOD::Type::CaloCluster ) { + ATH_MSG_VERBOSE(" Constituents are calo clusters."); + fillEperSamplingCluster(*jet, ePerSampling); + + } else if (ctype == xAOD::Type::ParticleFlow) { + ATH_MSG_VERBOSE(" Constituents are pflow objects."); + fillEperSamplingPFO(*jet, ePerSampling); - }else { - ATH_MSG_VERBOSE("Constituents are not calo clusters nor pflow objects."); - } + }else { + ATH_MSG_VERBOSE("Constituents are not calo clusters nor pflow objects."); + } - return 1; + } + return StatusCode::SUCCESS; } - -void JetCaloEnergies::fillEperSamplingCluster(xAOD::Jet& jet, std::vector<float> & ePerSampling ) const { +void JetCaloEnergies::fillEperSamplingCluster(const xAOD::Jet& jet, std::vector<float> & ePerSampling ) const { // loop over raw constituents size_t numConstit = jet.numConstituents(); for ( size_t i=0; i<numConstit; i++ ) { @@ -62,19 +85,16 @@ void JetCaloEnergies::fillEperSamplingCluster(xAOD::Jet& jet, std::vector<float> } const xAOD::CaloCluster* constit = static_cast<const xAOD::CaloCluster*>(jet.rawConstituent(i)); for ( size_t s= CaloSampling::PreSamplerB; s< CaloSampling::Unknown; s++ ) { - ePerSampling[s] += constit->eSample( (xAOD::CaloCluster::CaloSample) s ); - } + ePerSampling[s] += constit->eSample( (xAOD::CaloCluster::CaloSample) s ); } - static const xAOD::JetAttributeAccessor::AccessorWrapper<float>& emFracAcc = - *xAOD::JetAttributeAccessor::accessor< float >(xAOD::JetAttribute::EMFrac); - static const xAOD::JetAttributeAccessor::AccessorWrapper<float>& hecFracAcc = - *xAOD::JetAttributeAccessor::accessor< float >(xAOD::JetAttribute::HECFrac); - static const xAOD::JetAttributeAccessor::AccessorWrapper<float>& psFracAcc = - *xAOD::JetAttributeAccessor::accessor< float >(xAOD::JetAttribute::PSFrac); + } + SG::WriteDecorHandle<xAOD::JetContainer, float> emFracHandle(m_emFracKey); + SG::WriteDecorHandle<xAOD::JetContainer, float> hecFracHandle(m_hecFracKey); + SG::WriteDecorHandle<xAOD::JetContainer, float> psFracHandle(m_psFracKey); - emFracAcc(jet) = jet::JetCaloQualityUtils::emFraction( ePerSampling ); - hecFracAcc(jet) = jet::JetCaloQualityUtils::hecF( &jet ); - psFracAcc(jet) = jet::JetCaloQualityUtils::presamplerFraction( &jet ); + emFracHandle(jet) = jet::JetCaloQualityUtils::emFraction( ePerSampling ); + hecFracHandle(jet) = jet::JetCaloQualityUtils::hecF( &jet ); + psFracHandle(jet) = jet::JetCaloQualityUtils::presamplerFraction( &jet ); } #define FillESamplingPFO( LAYERNAME ) \ @@ -83,7 +103,7 @@ void JetCaloEnergies::fillEperSamplingCluster(xAOD::Jet& jet, std::vector<float> ePerSampling[CaloSampling::LAYERNAME] += E_##LAYERNAME; \ } -void JetCaloEnergies::fillEperSamplingPFO(xAOD::Jet & jet, std::vector<float> & ePerSampling ) const { +void JetCaloEnergies::fillEperSamplingPFO(const xAOD::Jet & jet, std::vector<float> & ePerSampling ) const { float emTot=0; float hecTot=0; @@ -144,26 +164,24 @@ void JetCaloEnergies::fillEperSamplingPFO(xAOD::Jet & jet, std::vector<float> & } } - static const xAOD::JetAttributeAccessor::AccessorWrapper<float>& emFracAcc = - *xAOD::JetAttributeAccessor::accessor< float >(xAOD::JetAttribute::EMFrac); + SG::WriteDecorHandle<xAOD::JetContainer, float> emFracHandle(m_emFracKey); if(eTot != 0.0){ - emFracAcc(jet) = emTot/eTot; + emFracHandle(jet) = emTot/eTot; /* * Ratio of EM layer calorimeter energy of neutrals to sum of all constituents * at EM scale (note charged PFO have an EM scale at track scale, and charged weights are ignored) * */ } else { - emFracAcc(jet) = 0.; + emFracHandle(jet) = 0.; } - static const xAOD::JetAttributeAccessor::AccessorWrapper<float>& hecFracAcc = - *xAOD::JetAttributeAccessor::accessor< float >(xAOD::JetAttribute::HECFrac); + SG::WriteDecorHandle<xAOD::JetContainer, float> hecFracHandle(m_hecFracKey); if (eTot != 0.0){ - hecFracAcc(jet) = hecTot/eTot; + hecFracHandle(jet) = hecTot/eTot; } else{ - hecFracAcc(jet) = 0.; + hecFracHandle(jet) = 0.; } } diff --git a/Reconstruction/Jet/JetMomentTools/Root/JetForwardJvtTool.cxx b/Reconstruction/Jet/JetMomentTools/Root/JetForwardJvtTool.cxx index d11c475a09f16d6cae55d5991b3f48583dd429f0..f4fac18e6572d12d47a9de738687e4e018319fb7 100644 --- a/Reconstruction/Jet/JetMomentTools/Root/JetForwardJvtTool.cxx +++ b/Reconstruction/Jet/JetMomentTools/Root/JetForwardJvtTool.cxx @@ -1,7 +1,7 @@ ///////////////////////// -*- C++ -*- ///////////////////////////// /* - Copyright (C) 2002-2017, 2019 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration */ // JetForwardJvtTool.cxx @@ -13,14 +13,6 @@ #include "JetMomentTools/JetForwardJvtTool.h" // Jet EDM -#include "xAODJet/JetAttributes.h" - -// Shallow copy -//#include "xAODCore/ShallowCopy.h" - - static SG::AuxElement::Decorator<char> isHS("isJvtHS"); - static SG::AuxElement::Decorator<char> isPU("isJvtPU"); - static SG::AuxElement::Decorator<float> fjvt_dec("fJvt"); /////////////////////////////////////////////////////////////////// // Public methods: @@ -29,27 +21,8 @@ // Constructors //////////////// JetForwardJvtTool::JetForwardJvtTool(const std::string& name) : - AsgTool(name), - m_fjvtThresh(15e3) - { - declareProperty("OverlapDec", m_orLabel = "" ); - declareProperty("OutputDec", m_outLabel = "passFJVT" ); - declareProperty("EtaThresh", m_etaThresh = 2.5 ); - declareProperty("ForwardMinPt", m_forwardMinPt = 20e3 ); - declareProperty("ForwardMaxPt", m_forwardMaxPt = 50e3 ); - declareProperty("CentralMinPt", m_centerMinPt = 20e3 ); - declareProperty("CentralMaxPt", m_centerMaxPt = -1 ); - declareProperty("CentralJvtThresh", m_centerJvtThresh = 0.14 ); - declareProperty("JvtMomentName", m_jvtMomentName = "Jvt" ); - declareProperty("CentralDrptThresh", m_centerDrptThresh = 0.2 ); - declareProperty("CentralMaxStochPt", m_maxStochPt = 35e3 ); - declareProperty("JetScaleFactor", m_jetScaleFactor = 0.4 ); - //declareProperty("FJVTThreshold", m_fjvtThresh = 15e3 );//15GeV->92%,11GeV->85% - declareProperty("UseTightOP", m_tightOP = false );//Tight or Loose - - declareProperty("VertexContainer", m_vertexContainer_key="PrimaryVertices"); - declareProperty("Met_Track", m_trkMET_key="Met_Track"); - + AsgTool(name) { + declareInterface<IJetDecorator>(this); } // Destructor @@ -61,15 +34,34 @@ //////////////////////////// StatusCode JetForwardJvtTool::initialize() { - ATH_MSG_DEBUG("initializing version with data handles"); ATH_MSG_INFO ("Initializing " << name() << "..."); if (m_tightOP) m_fjvtThresh = 0.4; else m_fjvtThresh = 0.5; - if (m_orLabel!="") m_Dec_OR = new SG::AuxElement::Decorator<char>(m_orLabel); - m_Dec_out = new SG::AuxElement::Decorator<char>(m_outLabel); - ATH_CHECK(m_vertexContainer_key.initialize()); - ATH_CHECK(m_trkMET_key.initialize()); + ATH_CHECK(m_vertexContainerName.initialize()); + ATH_CHECK(m_trkMETName.initialize()); + + if(m_jetContainerName.empty()) { + ATH_MSG_ERROR("JetForwardJvtTool needs to have its input jet container configured!"); + return StatusCode::FAILURE; + } + m_orKey = m_jetContainerName + "." + m_orKey.key(); + m_outKey = m_jetContainerName + "." + m_outKey.key(); + m_isHSKey = m_jetContainerName + "." + m_isHSKey.key(); + m_isPUKey = m_jetContainerName + "." + m_isPUKey.key(); + m_fjvtDecKey = m_jetContainerName + "." + m_fjvtDecKey.key(); + m_widthKey = m_jetContainerName + "." + m_widthKey.key(); + m_jvtMomentKey = m_jetContainerName + "." + m_jvtMomentKey.key(); + m_sumPtsKey = m_jetContainerName + "." + m_sumPtsKey.key(); + + ATH_CHECK(m_orKey.initialize()); + ATH_CHECK(m_outKey.initialize()); + ATH_CHECK(m_isHSKey.initialize()); + ATH_CHECK(m_isPUKey.initialize()); + ATH_CHECK(m_fjvtDecKey.initialize()); + ATH_CHECK(m_widthKey.initialize()); + ATH_CHECK(m_jvtMomentKey.initialize()); + ATH_CHECK(m_sumPtsKey.initialize()); return StatusCode::SUCCESS; } @@ -80,16 +72,19 @@ return StatusCode::SUCCESS; } - StatusCode JetForwardJvtTool::modify(xAOD::JetContainer& jetCont) const { + StatusCode JetForwardJvtTool::decorate(const xAOD::JetContainer& jetCont) const { + SG::WriteDecorHandle<xAOD::JetContainer, char> outHandle(m_outKey); + SG::WriteDecorHandle<xAOD::JetContainer, float> fjvtDecHandle(m_fjvtDecKey); + getPV(); if (jetCont.size() > 0) calculateVertexMomenta(&jetCont); for(const auto& jetF : jetCont) { - (*m_Dec_out)(*jetF) = 1; - fjvt_dec(*jetF) = 0; + outHandle(*jetF) = 1; + fjvtDecHandle(*jetF) = 0; if (!forwardJet(jetF)) continue; double fjvt = getFJVT(jetF)/jetF->pt(); - if (fjvt>m_fjvtThresh) (*m_Dec_out)(*jetF) = 0; - fjvt_dec(*jetF) = fjvt; + if (fjvt>m_fjvtThresh) outHandle(*jetF) = 0; + fjvtDecHandle(*jetF) = fjvt; } return StatusCode::SUCCESS; } @@ -109,28 +104,25 @@ void JetForwardJvtTool::calculateVertexMomenta(const xAOD::JetContainer *jets) const { m_pileupMomenta.clear(); - auto trkMETContainer = SG::makeHandle (m_trkMET_key); - if (!trkMETContainer.isValid()){ - ATH_MSG_WARNING("Invalid xAOD::MissingETContainer datahandle"); + SG::ReadHandle<xAOD::MissingETContainer> trkMetHandle(m_trkMETName); + SG::ReadHandle<xAOD::VertexContainer> vertexContainerHandle(m_vertexContainerName); + if( !trkMetHandle.isValid() ) { + ATH_MSG_WARNING("xAOD::MissingETContainer " << m_trkMETName.key() << " is invalid"); return; } - auto trkMet = trkMETContainer.cptr(); - - auto vertexContainer = SG::makeHandle (m_vertexContainer_key); - if (!vertexContainer.isValid()){ - ATH_MSG_WARNING("Invalid xAOD::VertexContainer datahandle"); + if( !vertexContainerHandle.isValid() ) { + ATH_MSG_WARNING("xAOD::VertexContainer " << m_vertexContainerName.key() << " is invalid"); return; } - auto vxCont = vertexContainer.cptr(); - for(const auto& vx : *vxCont) { + for(const auto& vx : *vertexContainerHandle) { if(vx->vertexType()!=xAOD::VxType::PriVtx && vx->vertexType()!=xAOD::VxType::PileUp) continue; TString vname = "PVTrack_vx"; vname += vx->index(); m_pileupMomenta.push_back(\ (vx->index()==m_pvind ? \ 0:\ - -(1./m_jetScaleFactor))*TVector2(0.5*(*trkMet)[vname.Data()]->mpx(),0.5*(*trkMet)[vname.Data()]->mpy())); + -(1./m_jetScaleFactor))*TVector2(0.5*(*trkMetHandle)[vname.Data()]->mpx(),0.5*(*trkMetHandle)[vname.Data()]->mpy())); } for (const auto& jet : *jets) { @@ -144,7 +136,8 @@ float Width = 0; float CWidth = 0; float ptsum = 0; - jet->getAttribute(xAOD::JetAttribute::Width,Width); + SG::ReadDecorHandle<xAOD::JetContainer, float> widthHandle(m_widthKey); + Width = widthHandle(*jet); xAOD::JetConstituentVector constvec = jet->getConstituents(); for (xAOD::JetConstituentVector::iterator it = constvec.begin(); it != constvec.end(); it++) { const xAOD::CaloCluster *cl = static_cast<const xAOD::CaloCluster*>((*it)->rawConstituent()); @@ -164,11 +157,13 @@ } bool JetForwardJvtTool::centralJet(const xAOD::Jet *jet) const { + SG::ReadDecorHandle<xAOD::JetContainer, char> orHandle(m_orKey); if (fabs(jet->eta())>m_etaThresh) return false; if (jet->pt()<m_centerMinPt || (m_centerMaxPt>0 && jet->pt()>m_centerMaxPt)) return false; - if (m_Dec_OR && !(*m_Dec_OR)(*jet)) return false; + if ( !orHandle(*jet) ) return false; float jvt = 0; - jet->getAttribute<float>(m_jvtMomentName,jvt); + SG::ReadDecorHandle<xAOD::JetContainer, float> jvtMomentHandle(m_jvtMomentKey); + jvt = jvtMomentHandle(*jet); if (jvt>m_centerJvtThresh) return false; if (jet->pt()<m_maxStochPt && getDrpt(jet)<m_centerDrptThresh) return false; return true; @@ -176,7 +171,8 @@ int JetForwardJvtTool::getJetVertex(const xAOD::Jet *jet) const { std::vector<float> sumpts; - jet->getAttribute<std::vector<float> >("SumPtTrkPt500",sumpts); + SG::ReadDecorHandle<xAOD::JetContainer, std::vector<float> > sumPtsHandle(m_sumPtsKey); + sumpts = sumPtsHandle(*jet); double firstVal = 0; int bestMatch = -1; for (size_t i = 0; i < sumpts.size(); i++) { @@ -190,7 +186,8 @@ float JetForwardJvtTool::getDrpt(const xAOD::Jet *jet) const { std::vector<float> sumpts; - jet->getAttribute<std::vector<float> >("SumPtTrkPt500",sumpts); + SG::ReadDecorHandle<xAOD::JetContainer, std::vector<float> > sumPtsHandle(m_sumPtsKey); + sumpts = sumPtsHandle(*jet); if (sumpts.size()<2) return 0; std::nth_element(sumpts.begin(),sumpts.begin()+sumpts.size()/2,sumpts.end(),std::greater<int>()); @@ -202,7 +199,7 @@ void JetForwardJvtTool::getPV() const { - auto vertexContainer = SG::makeHandle (m_vertexContainer_key); + auto vertexContainer = SG::makeHandle (m_vertexContainerName); if (!vertexContainer.isValid()){ ATH_MSG_WARNING("Invalid xAOD::VertexContainer datahandle"); return; @@ -222,6 +219,8 @@ } StatusCode JetForwardJvtTool::tagTruth(const xAOD::JetContainer *jets,const xAOD::JetContainer *truthJets) { + SG::WriteDecorHandle<xAOD::JetContainer, bool> isHSHandle(m_isHSKey); + SG::WriteDecorHandle<xAOD::JetContainer, bool> isPUHandle(m_isPUKey); for(const auto& jet : *jets) { bool ishs = false; bool ispu = true; @@ -229,8 +228,8 @@ if (tjet->p4().DeltaR(jet->p4())<0.3 && tjet->pt()>10e3) ishs = true; if (tjet->p4().DeltaR(jet->p4())<0.6) ispu = false; } - isHS(*jet)=ishs; - isPU(*jet)=ispu; + isHSHandle(*jet)=ishs; + isPUHandle(*jet)=ispu; } return StatusCode::SUCCESS; } diff --git a/Reconstruction/Jet/JetMomentTools/Root/JetForwardJvtToolBDT.cxx b/Reconstruction/Jet/JetMomentTools/Root/JetForwardJvtToolBDT.cxx new file mode 100644 index 0000000000000000000000000000000000000000..d56444edd38b1940c4b279a78f22a9315b0414c7 --- /dev/null +++ b/Reconstruction/Jet/JetMomentTools/Root/JetForwardJvtToolBDT.cxx @@ -0,0 +1,450 @@ +///////////////////////// -*- C++ -*- ///////////////////////////// +/* + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration +*/ +// JetForwardJvtToolBDT.cxx +// Implementation file for class JetForwardJvtToolBDT +// Author: Louis Portales <louis.portales@cern.ch> +/////////////////////////////////////////////////////////////////// + +// JetForwardJvtToolBDT includes +#include "JetMomentTools/JetForwardJvtToolBDT.h" +// Jet EDM +#include "xAODJet/JetAttributes.h" +#include "xAODMetaData/FileMetaData.h" +#include "PathResolver/PathResolver.h" + + +const double GeV = 1000.; +/////////////////////////////////////////////////////////////////// +// Public methods: +/////////////////////////////////////////////////////////////////// + +// Constructors +//////////////// +JetForwardJvtToolBDT::JetForwardJvtToolBDT(const std::string& name) : + asg::AsgTool(name) +{ + declareInterface<IJetDecorator>(this); +} + +// Destructor +/////////////// +JetForwardJvtToolBDT::~JetForwardJvtToolBDT(){} + +// Athena algtool's Hooks +//////////////////////////// +StatusCode JetForwardJvtToolBDT::initialize() +{ + ATH_MSG_INFO ("Initializing " << name() << "..."); + + if(m_isAna){ + // -- Retrieve MVfJVT WP configFile ONLY if tool used in 'Analysis mode' + std::string filename = PathResolverFindCalibFile(m_configDir+m_wpFile); + if (filename.empty()){ + ATH_MSG_ERROR ( "Could NOT resolve file name " << m_wpFile); + return StatusCode::FAILURE; + } else{ + ATH_MSG_INFO(" Config Files Path found = "<<filename); + } + + // -- Retrieve WP histograms + m_wpFileIn = std::make_unique<TFile> (filename.c_str(),"read"); + + if ( m_OP=="TIGHTER") { + m_mvfjvtThresh = std::unique_ptr< TH3D >( dynamic_cast<TH3D*>( m_wpFileIn->Get( "MVfJVT_tighter" ) ) ); + } else if ( m_OP=="TIGHT" ) { + m_mvfjvtThresh = std::unique_ptr< TH3D >( dynamic_cast<TH3D*>( m_wpFileIn->Get( "MVfJVT_tight" ) ) ); + } else if ( m_OP=="DEFAULT" || m_OP=="LOOSE" ) { + m_mvfjvtThresh = std::unique_ptr< TH3D >( dynamic_cast<TH3D*>( m_wpFileIn->Get( "MVfJVT_loose" ) ) ); + } else { + ATH_MSG_ERROR(m_OP << " working point doesn't exist." ); + return StatusCode::FAILURE; + } + m_mvfjvtThresh->SetDirectory(0); + m_wpFileIn->Close(); + } + + // -- Setup the tagger + m_MVreader = std::make_unique< TMVA::Reader > ( "Silent" ); + float fjvt,width,time,cllambda2,cletawidth,cle,cliso,clemprob; + m_MVreader->AddVariable( "fjvtdist", &fjvt ); + m_MVreader->AddVariable( "Width_jet", &width ); + m_MVreader->AddVariable( "timedist", &time ); + m_MVreader->AddVariable( "jet_LeadingClusterSecondLambda", &cllambda2 ); + m_MVreader->AddVariable( "cl_etaWidthLead", &cletawidth ); + m_MVreader->AddVariable( "clsum_e", &cle ); + m_MVreader->AddVariable( "cl_ISOLATIONsumE", &cliso ); + m_MVreader->AddVariable( "cl_EM_PROBABILITYsumE", &clemprob ); + for(unsigned int i = 0; i<m_MVconfig.size(); ++i) m_MVreader->BookMVA(TString::Format("BDT_%i",i+1),PathResolverFindCalibFile(m_configDir+m_MVconfig.value().at(i))); + + // "passMVfJVT" flag + m_outMVKey = m_jetContainerName + "." + m_outMVKey.key(); + ATH_CHECK(m_outMVKey.initialize()); + + // Moments values + m_mvfjvtKey = m_jetContainerName + "." + m_mvfjvtKey.key(); + ATH_CHECK(m_mvfjvtKey.initialize()); + + m_cllambda2Key = m_jetContainerName + "." + m_cllambda2Key.key(); + m_clwidthKey = m_jetContainerName + "." + m_clwidthKey.key(); + m_clisoKey = m_jetContainerName + "." + m_clisoKey.key(); + m_clemprobKey = m_jetContainerName + "." + m_clemprobKey.key(); + m_cleKey = m_jetContainerName + "." + m_cleKey.key(); + + ATH_CHECK(m_cllambda2Key.initialize()); + ATH_CHECK(m_clwidthKey.initialize()); + ATH_CHECK(m_clisoKey.initialize()); + ATH_CHECK(m_clemprobKey.initialize()); + ATH_CHECK(m_cleKey.initialize()); + + m_lcllambda2Key = m_jetContainerName + "." + m_lcllambda2Key.key(); + m_lcllambda2NTKey = m_jetContainerName + "." + m_lcllambda2NTKey.key(); + m_lclwidthKey = m_jetContainerName + "." + m_lclwidthKey.key(); + m_lclisoKey = m_jetContainerName + "." + m_lclisoKey.key(); + m_lclemprobKey = m_jetContainerName + "." + m_lclemprobKey.key(); + m_lcleKey = m_jetContainerName + "." + m_lcleKey.key(); + + ATH_CHECK(m_lcllambda2Key.initialize()); + ATH_CHECK(m_lcllambda2NTKey.initialize()); + ATH_CHECK(m_lclwidthKey.initialize()); + ATH_CHECK(m_lclisoKey.initialize()); + ATH_CHECK(m_lclemprobKey.initialize()); + ATH_CHECK(m_lcleKey.initialize()); + + ATH_CHECK(m_eventInfoKey.initialize()); + ATH_CHECK(m_vertexContainerKey.initialize()); + ATH_CHECK(m_caloClusterContainerKey.initialize()); + ATH_CHECK(m_trkMetKey.initialize()); + + // Truth information + m_isHSKey = m_jetContainerName + "." + m_isHSKey.key(); + m_isPUKey = m_jetContainerName + "." + m_isPUKey.key(); + + ATH_CHECK(m_isHSKey.initialize()); + ATH_CHECK(m_isPUKey.initialize()); + + return StatusCode::SUCCESS; +} + + +StatusCode JetForwardJvtToolBDT::decorate(const xAOD::JetContainer& jetCont) const { + + // -- Retrieve PV index if not provided by user + //pvind = (m_pvind.value()==-1) ? getPV() : m_pvind; + int pvind = m_pvind.value(); + if(pvind == -1) pvind = getPV(); + + ATH_MSG_DEBUG("In JetForwardJvtToolBDT::modify: PV index = " << pvind); + if( pvind == -1 ){ + ATH_MSG_WARNING( "Something went wrong with the HS primary vertex identification." ); + return StatusCode::FAILURE; + } + + SG::WriteDecorHandle<xAOD::JetContainer, char> outMVHandle(m_outMVKey); + SG::WriteDecorHandle<xAOD::JetContainer, float> mvfjvtHandle(m_mvfjvtKey); + + SG::WriteDecorHandle<xAOD::JetContainer, float> cllambda2Handle(m_cllambda2Key); + SG::WriteDecorHandle<xAOD::JetContainer, float> clwidthHandle(m_clwidthKey); + SG::WriteDecorHandle<xAOD::JetContainer, float> clisoHandle(m_clisoKey); + SG::WriteDecorHandle<xAOD::JetContainer, float> clemprobHandle(m_clemprobKey); + SG::WriteDecorHandle<xAOD::JetContainer, float> cleHandle(m_cleKey); + std::vector<TVector2> pileupMomenta; + for(const xAOD::Jet *jetF : jetCont) { + + float mvfjvt = -2; + outMVHandle(*jetF) = 1; + cllambda2Handle(*jetF) = 0; + clwidthHandle(*jetF) = 0; + cleHandle(*jetF) = 0; + clisoHandle(*jetF) = 0; + clemprobHandle(*jetF) = 0; + + // -- Get PU vertices momenta sums, then compute tagger value for forward jets + if ( forwardJet(jetF) ){ + if( pileupMomenta.size()==0 ) { + pileupMomenta = calculateVertexMomenta(&jetCont, pvind); + if( pileupMomenta.size()==0 ) { + ATH_MSG_DEBUG( "pileupMomenta is empty, this can happen for events with no PU vertices. fJVT won't be computed for this event and will be set to 0 instead." ); + mvfjvtHandle(*jetF) = 0; + continue; + } + } + mvfjvt = getMVfJVT(jetF, pvind, pileupMomenta); + if(m_isAna) outMVHandle(*jetF) = passMVfJVT( mvfjvt, jetF->pt()/(GeV), fabs(jetF->eta()) ); + mvfjvtHandle(*jetF) = mvfjvt; + } + } + return StatusCode::SUCCESS; + +} + + +float JetForwardJvtToolBDT::getFJVT(const xAOD::Jet *jet, int pvind, std::vector<TVector2> pileupMomenta) const { + + TVector2 fjet(-jet->pt()*cos(jet->phi()),-jet->pt()*sin(jet->phi())); + double fjvt = 0; + ATH_MSG_DEBUG("In JetForwardJvtToolBDT::getFJVT -----> Starting looping on vertices (pileupMomenta.size() = "<<pileupMomenta.size()); + for (size_t pui = 0; pui < pileupMomenta.size(); pui++) { + if (pui!=(size_t)pvind){ + double projection = pileupMomenta[pui]*fjet/fjet.Mod(); + if (projection>fjvt) fjvt = projection; + } + } + return fjvt; +} + + +float JetForwardJvtToolBDT::getMVfJVT(const xAOD::Jet *jet, int pvind, std::vector<TVector2> pileupMomenta) const { + + if(m_isAna && !m_getTagger) return jet->auxdata<float>("DFCommonJets_MVfJVT"); + + StatusCode sc = getInputs(jet); + if( sc.isFailure() ) { + ATH_MSG_WARNING(" Could not calculate BDT inputs"); + return -2; + } + + SG::ReadHandle<xAOD::EventInfo> eventInfoHandle(m_eventInfoKey); + if ( !eventInfoHandle.isValid() ) { + ATH_MSG_WARNING(" xAOD::EventInfo " << m_eventInfoKey.key() << "is invalid"); + return -2; + } + float mu = eventInfoHandle->actualInteractionsPerCrossing(); + + if (!forwardJet(jet)) return -2; + + SG::WriteDecorHandle<xAOD::JetContainer, float> cllambda2Handle(m_cllambda2Key); + SG::WriteDecorHandle<xAOD::JetContainer, float> clwidthHandle(m_clwidthKey); + SG::WriteDecorHandle<xAOD::JetContainer, float> clisoHandle(m_clisoKey); + SG::WriteDecorHandle<xAOD::JetContainer, float> clemprobHandle(m_clemprobKey); + SG::WriteDecorHandle<xAOD::JetContainer, float> cleHandle(m_cleKey); + + std::vector<float> MVinputs; + MVinputs.push_back( getFJVT(jet, pvind, pileupMomenta)/jet->pt() ); + MVinputs.push_back( jet->getAttribute<float>("Width") ); + MVinputs.push_back( jet->getAttribute<float>("Timing") ); + MVinputs.push_back( cllambda2Handle(*jet) ); + MVinputs.push_back( clwidthHandle(*jet) ); + MVinputs.push_back( cleHandle(*jet) ); + MVinputs.push_back( clisoHandle(*jet) ); + MVinputs.push_back( clemprobHandle(*jet) ); + + float pt = jet->pt()/(GeV); + float eta = fabs(jet->eta()); + + float score = -2.; + if ( pt < 30. && pt >= 20. && eta >= 3.2 && mu>=50. ) score = m_MVreader->EvaluateMVA( MVinputs, "BDT_1" ,1.); + else if ( pt < 30. && pt >= 20. && eta < 3.2 && mu>=50. ) score = m_MVreader->EvaluateMVA( MVinputs, "BDT_2" ,1.); + else if ( pt < 40. && pt >= 30. && eta >= 3.2 && mu>=50. ) score = m_MVreader->EvaluateMVA( MVinputs, "BDT_3" ,1.); + else if ( pt < 40. && pt >= 30. && eta < 3.2 && mu>=50. ) score = m_MVreader->EvaluateMVA( MVinputs, "BDT_4" ,1.); + else if ( pt < 50. && pt >= 40. && eta >= 3.2 && mu>=50. ) score = m_MVreader->EvaluateMVA( MVinputs, "BDT_5" ,1.); + else if ( pt < 50. && pt >= 40. && eta < 3.2 && mu>=50. ) score = m_MVreader->EvaluateMVA( MVinputs, "BDT_6" ,1.); + else if ( pt < 120. && pt >= 50. && eta >= 3.2 && mu>=50. ) score = m_MVreader->EvaluateMVA( MVinputs, "BDT_7" ,1.); + else if ( pt < 120. && pt >= 50. && eta < 3.2 && mu>=50. ) score = m_MVreader->EvaluateMVA( MVinputs, "BDT_8" ,1.); + else if ( pt < 30. && pt >= 20. && eta >= 3.2 && mu<50. ) score = m_MVreader->EvaluateMVA( MVinputs, "BDT_9" ,1.); + else if ( pt < 30. && pt >= 20. && eta < 3.2 && mu<50. ) score = m_MVreader->EvaluateMVA( MVinputs, "BDT_10" ,1.); + else if ( pt < 40. && pt >= 30. && eta >= 3.2 && mu<50. ) score = m_MVreader->EvaluateMVA( MVinputs, "BDT_11" ,1.); + else if ( pt < 40. && pt >= 30. && eta < 3.2 && mu<50. ) score = m_MVreader->EvaluateMVA( MVinputs, "BDT_12" ,1.); + else if ( pt < 50. && pt >= 40. && eta >= 3.2 && mu<50. ) score = m_MVreader->EvaluateMVA( MVinputs, "BDT_13" ,1.); + else if ( pt < 50. && pt >= 40. && eta < 3.2 && mu<50. ) score = m_MVreader->EvaluateMVA( MVinputs, "BDT_14" ,1.); + else if ( pt < 120. && pt >= 50. && eta >= 3.2 && mu<50. ) score = m_MVreader->EvaluateMVA( MVinputs, "BDT_15" ,1.); + else if ( pt < 120. && pt >= 50. && eta < 3.2 && mu<50. ) score = m_MVreader->EvaluateMVA( MVinputs, "BDT_16" ,1.); + + ATH_MSG_DEBUG("pt = " << pt << " | eta = " << eta << " | mu = " << mu << " || MVfJVT = " << score ); + + return score; +} + +bool JetForwardJvtToolBDT::passMVfJVT( float mvfjvt, float pt, float eta ) const { + + double mvfjvtThresh = -999.; + + SG::ReadHandle<xAOD::EventInfo> eventInfoHandle(m_eventInfoKey); + if ( !eventInfoHandle.isValid() ) { + ATH_MSG_WARNING(" xAOD::EventInfo " << m_eventInfoKey.key() << "is invalid"); + return true; + } + + float mu = eventInfoHandle->actualInteractionsPerCrossing(); + + // -- Grab WP from histogram + mvfjvtThresh = m_mvfjvtThresh->GetBinContent(m_mvfjvtThresh->GetXaxis()->FindBin(pt), + m_mvfjvtThresh->GetYaxis()->FindBin(eta), + m_mvfjvtThresh->GetZaxis()->FindBin(mu)); + + if(mvfjvt==-2 || mvfjvt>mvfjvtThresh) return true; + else return false; + +} + +StatusCode JetForwardJvtToolBDT::getInputs(const xAOD::Jet *jet) const { + SG::WriteDecorHandle<xAOD::JetContainer, float> cllambda2Handle(m_cllambda2Key); + SG::WriteDecorHandle<xAOD::JetContainer, float> clwidthHandle(m_clwidthKey); + SG::WriteDecorHandle<xAOD::JetContainer, float> clisoHandle(m_clisoKey); + SG::WriteDecorHandle<xAOD::JetContainer, float> clemprobHandle(m_clemprobKey); + SG::WriteDecorHandle<xAOD::JetContainer, float> cleHandle(m_cleKey); + + if(!m_getTagger){ + SG::ReadDecorHandle<xAOD::JetContainer, float> lcllambda2NTHandle(m_lcllambda2NTKey); + cllambda2Handle(*jet) = lcllambda2NTHandle(*jet); + + // -- Additional variables computed from cluster information + SG::ReadHandle<xAOD::CaloClusterContainer> clustersHandle(m_caloClusterContainerKey); + if( !clustersHandle.isValid() ) { + ATH_MSG_ERROR(" xAOD::CaloClusterContainer " << m_caloClusterContainerKey.key() << "is invalid"); + return StatusCode::FAILURE; + } + + int ind = 0; + float maxpt = 0; + float cle1 = 0; + float cliso1 = 0; + float clemprob1 = 0; + float cle2 = 0; + + // Loop over clusters within DeltaR<0.6 of the jet axis to compute the (energy-weighted) moment sums used in the BDT definitions + for (const xAOD::CaloCluster *cl: *clustersHandle) { + if(cl->p4().DeltaR(jet->p4())>0.6) continue; + cle1 += cl->e(); + cle2 += cl->e()*cl->e(); + cliso1 += cl->auxdata<float>("ISOLATION")*cl->e()*cl->e(); + clemprob1 += cl->auxdata<float>("EM_PROBABILITY")*cl->e()*cl->e(); + if(cl->rawE()/cosh(cl->rawEta()) > maxpt){ + maxpt = cl->rawE()/cosh(cl->rawEta()); + ind = cl->index(); + } + } + const xAOD::CaloCluster *cl = clustersHandle->at(ind); + clwidthHandle(*jet) = TMath::CosH(cl->rawEta()) * TMath::ATan2( TMath::Sqrt(cl->auxdata<float>("SECOND_R")), + cl->auxdata<float>("CENTER_MAG")); + + cleHandle(*jet) = cle1; + clisoHandle(*jet)= cliso1/cle2; + clemprobHandle(*jet) =clemprob1/cle2; + + } else { + SG::ReadDecorHandle<xAOD::JetContainer, float> lcllambda2Handle(m_lcllambda2Key); + SG::ReadDecorHandle<xAOD::JetContainer, float> lclwidthHandle(m_lclwidthKey); + SG::ReadDecorHandle<xAOD::JetContainer, float> lclisoHandle(m_lclisoKey); + SG::ReadDecorHandle<xAOD::JetContainer, float> lclemprobHandle(m_lclemprobKey); + SG::ReadDecorHandle<xAOD::JetContainer, float> lcleHandle(m_lcleKey); + + cllambda2Handle(*jet) = lcllambda2Handle(*jet); + clwidthHandle(*jet) = lclwidthHandle(*jet); + clisoHandle(*jet) = lclisoHandle(*jet); + clemprobHandle(*jet) = lclemprobHandle(*jet); + cleHandle(*jet) = lcleHandle(*jet); + } + return StatusCode::SUCCESS; +} + +std::vector<TVector2> JetForwardJvtToolBDT::calculateVertexMomenta(const xAOD::JetContainer *jets, int pvind) const { + + std::vector<TVector2> pileupMomenta; + + SG::ReadHandle<xAOD::MissingETContainer> trkMetHandle(m_trkMetKey); + if( !trkMetHandle.isValid() ) { + ATH_MSG_WARNING(" xAOD::MissingETContainer " << m_trkMetKey.key() << "is invalid"); + return pileupMomenta; + } + SG::ReadHandle<xAOD::VertexContainer> vxContHandle(m_vertexContainerKey); + if( !vxContHandle.isValid() ) { + ATH_MSG_WARNING(" xAOD::VertexContainer " << m_vertexContainerKey.key() << "is invalid"); + return pileupMomenta; + } + ATH_MSG_DEBUG("In JetForwardJvtToolBDT::calculateVertexMomenta : Starting vertex loop "); + for(const xAOD::Vertex *vx : *vxContHandle) { + ATH_MSG_DEBUG(" --> VertexType="<<vx->vertexType()); + if(vx->vertexType()!=xAOD::VxType::PriVtx && vx->vertexType()!=xAOD::VxType::PileUp) continue; + TString vname = "PVTrack_vx"; + vname += vx->index(); + pileupMomenta.push_back((vx->index()==(size_t)pvind?0:-(1./m_jetScaleFactor))*TVector2(0.5*(*trkMetHandle)[vname.Data()]->mpx(),0.5*(*trkMetHandle)[vname.Data()]->mpy())); + } + for (const xAOD::Jet *jet : *jets) { + if (!centralJet(jet)) continue; + int jetvert = getJetVertex(jet); + if (jetvert>=0) pileupMomenta[jetvert] += TVector2(0.5*jet->pt()*cos(jet->phi()),0.5*jet->pt()*sin(jet->phi())); + } + + return pileupMomenta; +} + +bool JetForwardJvtToolBDT::forwardJet(const xAOD::Jet *jet) const { + + if (fabs(jet->eta())<m_etaThresh) return false; + if (jet->pt()<m_forwardMinPt || jet->pt()>m_forwardMaxPt) return false; + return true; +} + +bool JetForwardJvtToolBDT::centralJet(const xAOD::Jet *jet) const { + + if (fabs(jet->eta())>m_etaThresh) return false; + if (jet->pt()<m_centerMinPt || (m_centerMaxPt>0 && jet->pt()>m_centerMaxPt)) return false; + float jvt = 0; + jet->getAttribute<float>(m_jvtMomentName,jvt); + if (jvt>m_centerJvtThresh) return false; + if (jet->pt()<m_maxStochPt && getDrpt(jet)<m_centerDrptThresh) return false; + return true; +} + +int JetForwardJvtToolBDT::getJetVertex(const xAOD::Jet *jet) const { + + std::vector<float> sumpts; + jet->getAttribute<std::vector<float> >("SumPtTrkPt500",sumpts); + double firstVal = 0; + int bestMatch = -1; + for (size_t i = 0; i < sumpts.size(); i++) { + if (sumpts[i]>firstVal) { + bestMatch = i; + firstVal = sumpts[i]; + } + } + return bestMatch; +} + +float JetForwardJvtToolBDT::getDrpt(const xAOD::Jet *jet) const { + + std::vector<float> sumpts; + jet->getAttribute<std::vector<float> >("SumPtTrkPt500",sumpts); + if (sumpts.size()<2) return 0; + + std::nth_element(sumpts.begin(),sumpts.begin()+sumpts.size()/2,sumpts.end(),std::greater<int>()); + double median = sumpts[sumpts.size()/2]; + std::nth_element(sumpts.begin(),sumpts.begin(),sumpts.end(),std::greater<int>()); + double max = sumpts[0]; + return (max-median)/jet->pt(); +} + +int JetForwardJvtToolBDT::getPV() const{ + + SG::ReadHandle<xAOD::VertexContainer> vxContHandle(m_vertexContainerKey); + if( !vxContHandle.isValid() ) { + ATH_MSG_WARNING(" xAOD::VertexContainer " << m_vertexContainerKey.key() << "is invalid"); + return 0; + } else { + ATH_MSG_DEBUG("Successfully retrieved primary vertex container"); + for(const xAOD::Vertex *vx : *vxContHandle) { + if(vx->vertexType()==xAOD::VxType::PriVtx) return vx->index(); + } + } + ATH_MSG_DEBUG("Couldn't identify the hard-scatter primary vertex (no vertex with \"vx->vertexType()==xAOD::VxType::PriVtx\" in the container)!"); + return 0; +} + +StatusCode JetForwardJvtToolBDT::tagTruth(const xAOD::JetContainer *jets,const xAOD::JetContainer *truthJets) { + SG::WriteDecorHandle<xAOD::JetContainer, char> isHSHandle(m_isHSKey); + SG::WriteDecorHandle<xAOD::JetContainer, char> isPUHandle(m_isPUKey); + + for(const xAOD::Jet *jet : *jets) { + bool ishs = false; + bool ispu = true; + for(const xAOD::Jet *tjet : *truthJets) { + if (tjet->p4().DeltaR(jet->p4())<0.3 && tjet->pt()>10e3) ishs = true; + if (tjet->p4().DeltaR(jet->p4())<0.6) ispu = false; + } + isHSHandle(*jet)=ishs; + isPUHandle(*jet)=ispu; + } + return StatusCode::SUCCESS; +} diff --git a/Reconstruction/Jet/JetMomentTools/Root/JetMomentToolsDict.h b/Reconstruction/Jet/JetMomentTools/Root/JetMomentToolsDict.h index 0904e1b560e8bf80e3e564fb819d148868ad44a1..56942c59ee9a7673168ea6a38ad99a16ea5d08ae 100644 --- a/Reconstruction/Jet/JetMomentTools/Root/JetMomentToolsDict.h +++ b/Reconstruction/Jet/JetMomentTools/Root/JetMomentToolsDict.h @@ -12,7 +12,6 @@ #include "JetMomentTools/JetECPSFractionTool.h" #include "JetMomentTools/JetForwardJvtTool.h" #include "JetMomentTools/JetLArHVTool.h" -#include "JetMomentTools/JetMuonSegmentMomentsTool.h" #include "JetMomentTools/JetOriginCorrectionTool.h" #include "JetMomentTools/JetPtAssociationTool.h" #include "JetMomentTools/JetTrackMomentsTool.h" diff --git a/Reconstruction/Jet/JetMomentTools/Root/JetMuonSegmentMomentsTool.cxx b/Reconstruction/Jet/JetMomentTools/Root/JetMuonSegmentMomentsTool.cxx deleted file mode 100644 index b781cb906e9e4f40581cd0172fa6f22f7c53496b..0000000000000000000000000000000000000000 --- a/Reconstruction/Jet/JetMomentTools/Root/JetMuonSegmentMomentsTool.cxx +++ /dev/null @@ -1,31 +0,0 @@ -/* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -*/ - - -#include "JetMomentTools/JetMuonSegmentMomentsTool.h" - -JetMuonSegmentMomentsTool::JetMuonSegmentMomentsTool(const std::string& name) - : JetModifierBase(name) - , m_assocMuonSegName("GhostMuonSegment") -{ - declareProperty("AssociatedMuonSegments", m_assocMuonSegName); -} - -int JetMuonSegmentMomentsTool::modifyJet(xAOD::Jet& jet) const { - // Get muon segments associated to the jet - // Note that there are often no segments - std::vector<const xAOD::IParticle*> segments; - bool havesegs = jet.getAssociatedObjects(m_assocMuonSegName, segments); - - if ( ! havesegs ) { - ATH_MSG_ERROR("Muon segments not found."); - return 1; - } else { - ATH_MSG_DEBUG("Jet has " << segments.size() << " muon segments."); - } - - return 0; -} - - diff --git a/Reconstruction/Jet/JetMomentTools/Root/JetOriginCorrectionTool.cxx b/Reconstruction/Jet/JetMomentTools/Root/JetOriginCorrectionTool.cxx index 038ff7295d15c3aa67b3563baf3f23f9991f8b79..705ce6a35908fc7779dfc3489cb45f58b97faed5 100644 --- a/Reconstruction/Jet/JetMomentTools/Root/JetOriginCorrectionTool.cxx +++ b/Reconstruction/Jet/JetMomentTools/Root/JetOriginCorrectionTool.cxx @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration */ // JetOriginCorrectionTool.cxx @@ -15,29 +15,44 @@ JetOriginCorrectionTool::JetOriginCorrectionTool(const std::string& myname) : asg::AsgTool(myname) { - declareProperty("OriginCorrectedName", m_correctionName="JetOriginConstitScaleMomentum"); - declareProperty("OnlyAssignPV", m_onlyAssignPV=false); - - declareProperty("VertexContainer", m_vertexContainer_key="PrimaryVertices"); - declareProperty("EventInfoName", m_eventInfo_key="EventInfo"); } //********************************************************************** StatusCode JetOriginCorrectionTool::initialize() { - ATH_MSG_DEBUG("initializing version with data handles"); - ATH_CHECK(m_vertexContainer_key.initialize()); ATH_CHECK(m_eventInfo_key.initialize()); + if(m_jetContainerName.empty()) { + ATH_MSG_ERROR("JetOriginCorrectionTool needs to have its input jet container configured!"); + return StatusCode::FAILURE; + } + + m_scaleMomentumPtKey = m_jetContainerName + "." + m_scaleMomentumName + "_" + m_scaleMomentumPtKey.key(); + m_scaleMomentumPhiKey = m_jetContainerName + "." + m_scaleMomentumName + "_" + m_scaleMomentumPhiKey.key(); + m_scaleMomentumEtaKey = m_jetContainerName + "." + m_scaleMomentumName + "_" + m_scaleMomentumEtaKey.key(); + m_scaleMomentumMKey = m_jetContainerName + "." + m_scaleMomentumName + "_" + m_scaleMomentumMKey.key(); + m_originVertexKey = m_jetContainerName + "." + m_originVertexKey.key(); + + ATH_CHECK(m_scaleMomentumPtKey.initialize()); + ATH_CHECK(m_scaleMomentumPhiKey.initialize()); + ATH_CHECK(m_scaleMomentumEtaKey.initialize()); + ATH_CHECK(m_scaleMomentumMKey.initialize()); + ATH_CHECK(m_originVertexKey.initialize()); return StatusCode::SUCCESS; } //********************************************************************** -StatusCode JetOriginCorrectionTool::modify(xAOD::JetContainer& jetCont) const { +StatusCode JetOriginCorrectionTool::decorate(const xAOD::JetContainer& jetCont) const { + SG::WriteDecorHandle<xAOD::JetContainer, float> scaleMomentumPtHandle(m_scaleMomentumPtKey); + SG::WriteDecorHandle<xAOD::JetContainer, float> scaleMomentumPhiHandle(m_scaleMomentumPhiKey); + SG::WriteDecorHandle<xAOD::JetContainer, float> scaleMomentumEtaHandle(m_scaleMomentumEtaKey); + SG::WriteDecorHandle<xAOD::JetContainer, float> scaleMomentumMHandle(m_scaleMomentumMKey); + SG::WriteDecorHandle<xAOD::JetContainer, ElementLink<xAOD::VertexContainer>> originVertexHandle(m_originVertexKey); + // static accessor for PV index access static SG::AuxElement::ConstAccessor<int> PVIndexAccessor("PVIndex"); @@ -48,9 +63,14 @@ StatusCode JetOriginCorrectionTool::modify(xAOD::JetContainer& jetCont) const { if (!handle.isValid()){ ATH_MSG_WARNING("Invalid VertexContainer datahandle: " << m_correctionName - << ": filling jet with null"); + << ": filling jet with -1"); xAOD::JetFourMom_t null; - for(xAOD::Jet * j : jetCont) j->setAttribute<xAOD::JetFourMom_t>(m_correctionName, null); + for(const xAOD::Jet * j : jetCont) { + scaleMomentumPtHandle(*j) = -1; + scaleMomentumPhiHandle(*j) = -1; + scaleMomentumEtaHandle(*j) = -1; + scaleMomentumMHandle(*j) = -1; + } return StatusCode::SUCCESS; } @@ -69,7 +89,7 @@ StatusCode JetOriginCorrectionTool::modify(xAOD::JetContainer& jetCont) const { auto eInfo = SG::makeHandle (m_eventInfo_key); if (!eInfo.isValid()){ - ATH_MSG_WARNING("Invalid eventInfo datahandle. Defaulting ti PV0 for " + ATH_MSG_WARNING("Invalid eventInfo datahandle. Defaulting to PV0 for " << m_correctionName); } else if (PVIndexAccessor.isAvailable(*(eInfo.cptr()))) { PVindex = PVIndexAccessor(*(eInfo.cptr())); @@ -77,10 +97,15 @@ StatusCode JetOriginCorrectionTool::modify(xAOD::JetContainer& jetCont) const { if (PVindex < 0 || static_cast<size_t>(PVindex) >= vxContainer->size()){ ATH_MSG_WARNING("Specified PV index of " - << PVindex << " is out of bounds. Filling jet with null " + << PVindex << " is out of bounds. Filling jet with -1" <<m_correctionName); xAOD::JetFourMom_t null; - for (xAOD::Jet* j : jetCont) j->setAttribute<xAOD::JetFourMom_t>(m_correctionName,null); + for (const xAOD::Jet* j : jetCont) { + scaleMomentumPtHandle(*j) = -1; + scaleMomentumPhiHandle(*j) = -1; + scaleMomentumEtaHandle(*j) = -1; + scaleMomentumMHandle(*j) = -1; + } return StatusCode::SUCCESS; } @@ -91,16 +116,19 @@ StatusCode JetOriginCorrectionTool::modify(xAOD::JetContainer& jetCont) const { const xAOD::Vertex *vx = vxContainer->at(PVindex); ATH_MSG_DEBUG(" correcting jets "); - for(xAOD::Jet * jet : jetCont){ + for(const xAOD::Jet * jet : jetCont){ ATH_MSG_DEBUG(" ----> jet "<< jet); ATH_MSG_DEBUG(" jet pT: "<< jet->pt()); if(!m_onlyAssignPV) { xAOD::JetFourMom_t fv = jet::clusterOriginCorrection(*jet,*vx); ATH_MSG_DEBUG(" " << m_correctionName << " pT: " << fv.pt()); - jet->setAttribute<xAOD::JetFourMom_t>(m_correctionName, fv); + scaleMomentumPtHandle(*jet) = fv.pt(); + scaleMomentumPhiHandle(*jet) = fv.phi(); + scaleMomentumEtaHandle(*jet) = fv.eta(); + scaleMomentumMHandle(*jet) = fv.M(); } - jet->setAssociatedObject("OriginVertex", vx); + originVertexHandle(*jet) = ElementLink<xAOD::VertexContainer>(*vxContainer, vxContainer->at(PVindex)->index()); } return StatusCode::SUCCESS; diff --git a/Reconstruction/Jet/JetMomentTools/Root/JetPtAssociationTool.cxx b/Reconstruction/Jet/JetMomentTools/Root/JetPtAssociationTool.cxx index 70531bb6fe00a24fb6ae2c33e440ceed30fe43e2..d44ab08d842db41d6c3605778b193cd123341bb1 100644 --- a/Reconstruction/Jet/JetMomentTools/Root/JetPtAssociationTool.cxx +++ b/Reconstruction/Jet/JetMomentTools/Root/JetPtAssociationTool.cxx @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration */ // JetPtAssociationTool.cxx @@ -16,9 +16,8 @@ using xAOD::JetConstituentVector; //********************************************************************** JetPtAssociationTool::JetPtAssociationTool(std::string myname) -: JetModifierBase(myname) { - declareProperty("AssociationName", m_aname); - declareProperty("InputContainer", m_jetContainer_key); + : asg::AsgTool(myname) { + declareInterface<IJetDecorator>(this); } @@ -26,62 +25,76 @@ JetPtAssociationTool::JetPtAssociationTool(std::string myname) StatusCode JetPtAssociationTool::initialize() { - ATH_MSG_DEBUG("initializing version with data handles"); - ATH_CHECK(m_jetContainer_key.initialize()); + if(m_jetContainerName.empty()) { + ATH_MSG_ERROR("JetPtAssociationTool needs to have its input jet container configured!"); + return StatusCode::FAILURE; + } + m_assocFracKey = m_jetContainerName + "." + m_aname + m_assocFracKey.key(); + m_assocLinkKey = m_jetContainerName + "." + m_aname + m_assocLinkKey.key(); + + ATH_CHECK(m_assocFracKey.initialize()); + ATH_CHECK(m_assocLinkKey.initialize()); + return StatusCode::SUCCESS; } //********************************************************************** -int JetPtAssociationTool::modifyJet(xAOD::Jet& jet) const { - ATH_MSG_DEBUG("Processing jet " << jet.index()); - APVector apins; +StatusCode JetPtAssociationTool::decorate(const xAOD::JetContainer& jets) const { - // Retrieve the association vector. - if ( ! jet.getAssociatedObjects(m_aname, apins) ) { - ATH_MSG_WARNING("Jet does not have association vector " << m_aname); - return 1; - } + SG::WriteDecorHandle<xAOD::JetContainer, float> assocFracHandle(m_assocFracKey); + SG::WriteDecorHandle<xAOD::JetContainer, ElementLink<JetContainer>> assocLinkHandle(m_assocLinkKey); + + for(const xAOD::Jet* jet : jets) { + ATH_MSG_DEBUG("Processing jet " << jet->index()); + APVector apins; - // Retrieve the container of jets to be matched. - auto handle = SG::makeHandle (m_jetContainer_key); - if (!handle.isValid()){ - ATH_MSG_WARNING("Matching jet container not found: " - << m_jetContainer_key.key()); - return 2; - } + // Retrieve the association vector. + if ( ! jet->getAssociatedObjects(m_aname, apins) ) { + ATH_MSG_WARNING("Jet does not have association vector " << m_aname); + return StatusCode::FAILURE; + } - auto pjets = handle.cptr(); + // Retrieve the container of jets to be matched. + SG::ReadHandle<xAOD::JetContainer> handle (m_jetContainerName); + if (!handle.isValid()){ + ATH_MSG_WARNING("Matching jet container not found: " + << m_jetContainerName); + return StatusCode::FAILURE; + } - // Match associated particle to jets. - FloatVector ptfs; - if ( ptfrac(apins, *pjets, ptfs) ) { - ATH_MSG_WARNING("Match of association vector to jets failed"); - return 3; - } - // Find the best match. - float frac = 0.0; - unsigned int ijet_matched = ptfs.size(); - for ( unsigned int ijet=0; ijet<ptfs.size(); ++ijet ) { - ATH_MSG_VERBOSE(" Pt fraction[" << ijet << "]: " << ptfs[ijet]); - if ( ptfs[ijet] > frac ) { - frac = ptfs[ijet]; - ijet_matched = ijet; + auto pjets = handle.cptr(); + + // Match associated particle to jets. + FloatVector ptfs; + if ( ptfrac(apins, *pjets, ptfs) ) { + ATH_MSG_WARNING("Match of association vector to jets failed"); + return StatusCode::FAILURE; } + // Find the best match. + float frac = 0.0; + unsigned int ijet_matched = ptfs.size(); + for ( unsigned int ijet=0; ijet<ptfs.size(); ++ijet ) { + ATH_MSG_VERBOSE(" Pt fraction[" << ijet << "]: " << ptfs[ijet]); + if ( ptfs[ijet] > frac ) { + frac = ptfs[ijet]; + ijet_matched = ijet; + } + } + ElementLink<JetContainer> el; + if ( ijet_matched < ptfs.size() ) { + el = ElementLink<JetContainer>(*pjets, ijet_matched); + } + string sfrac = m_aname + "AssociationFraction"; + string slink = m_aname + "AssociationLink"; + ATH_MSG_DEBUG(" Setting " << sfrac << " = " << frac); + assocFracHandle(*jet) = frac; + ATH_MSG_DEBUG(" Setting " << slink << " = " + << el.dataID() << "[" << el.index() << "]"); + assocLinkHandle(*jet) = el; } - ElementLink<JetContainer> el; - if ( ijet_matched < ptfs.size() ) { - el = ElementLink<JetContainer>(*pjets, ijet_matched); - } - string sfrac = m_aname + "AssociationFraction"; - string slink = m_aname + "AssociationLink"; - ATH_MSG_DEBUG(" Setting " << sfrac << " = " << frac); - jet.setAttribute(sfrac, frac); - ATH_MSG_DEBUG(" Setting " << slink << " = " - << el.dataID() << "[" << el.index() << "]"); - jet.setAttribute(slink, el); - return 0; + return StatusCode::SUCCESS; } //********************************************************************** diff --git a/Reconstruction/Jet/JetMomentTools/Root/JetTrackSumMomentsTool.cxx b/Reconstruction/Jet/JetMomentTools/Root/JetTrackSumMomentsTool.cxx index 5b90d1afd215b54b7592fd67ecb0c326b72e539f..b8c9ab45501acc130af86ccce10490e58165ca44 100644 --- a/Reconstruction/Jet/JetMomentTools/Root/JetTrackSumMomentsTool.cxx +++ b/Reconstruction/Jet/JetMomentTools/Root/JetTrackSumMomentsTool.cxx @@ -1,7 +1,7 @@ ///////////////////////// -*- C++ -*- //////////////////////////// /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration */ // JetTrackSumMomentsTool.cxx @@ -16,24 +16,13 @@ using xAOD::JetFourMom_t; //using xAOD::TrackParticle::FourMom_t; JetTrackSumMomentsTool::JetTrackSumMomentsTool(const std::string& name) - : JetModifierBase(name) - // , m_vertexContainer("") - , m_assocTracksName("") - , m_htsel("",this) - , m_vertexContainer_key("") - , m_tva_key("") -{ - declareProperty("VertexContainer",m_vertexContainer_key); - declareProperty("AssociatedTracks",m_assocTracksName); - declareProperty("TrackVertexAssociation",m_tva_key); - declareProperty("TrackSelector", m_htsel); - declareProperty("RequireTrackPV", m_requireTrackPV = true); + : asg::AsgTool(name) { + declareInterface<IJetDecorator>(this); } //********************************************************************** StatusCode JetTrackSumMomentsTool::initialize() { - ATH_MSG_DEBUG("initializing version with data handles"); ATH_MSG_INFO("Initializing JetTrackSumMomentsTool " << name()); if ( m_htsel.empty() ) { ATH_MSG_INFO(" No track selector."); @@ -44,20 +33,29 @@ StatusCode JetTrackSumMomentsTool::initialize() { ATH_CHECK(m_vertexContainer_key.initialize()); ATH_CHECK(m_tva_key.initialize()); + if(m_jetContainerName.empty()) { + ATH_MSG_ERROR("JetTrackSumMomentsTool needs to have its input jet container configured!"); + return StatusCode::FAILURE; + } + m_trackSumPtKey = m_jetContainerName + "." + m_trackSumPtKey.key(); + m_trackSumMassKey = m_jetContainerName + "." + m_trackSumMassKey.key(); + + ATH_CHECK(m_trackSumPtKey.initialize()); + ATH_CHECK(m_trackSumMassKey.initialize()); + return StatusCode::SUCCESS; } //********************************************************************** -int JetTrackSumMomentsTool::modifyJet(xAOD::Jet& jet) const { - +StatusCode JetTrackSumMomentsTool::decorate(const xAOD::JetContainer& jets) const { // Get input vertex collection auto handle_v = SG::makeHandle(m_vertexContainer_key); if (!handle_v.isValid()){ ATH_MSG_ERROR("Could not retrieve the VertexContainer: " << m_vertexContainer_key.key()); - return 1; + return StatusCode::FAILURE; } auto vertexContainer = handle_v.cptr(); @@ -67,34 +65,39 @@ int JetTrackSumMomentsTool::modifyJet(xAOD::Jet& jet) const { if (!handle_tva.isValid()){ ATH_MSG_ERROR("Could not retrieve the TrackVertexAssociation: " << m_tva_key.key()); - return 2; + return StatusCode::FAILURE; } auto tva = handle_tva.cptr(); + SG::WriteDecorHandle<xAOD::JetContainer, float> trackSumPtHandle(m_trackSumPtKey); + SG::WriteDecorHandle<xAOD::JetContainer, float> trackSumMassHandle(m_trackSumMassKey); + // Get the tracks associated to the jet // Note that there may be no tracks - this is both normal and an error case std::vector<const xAOD::TrackParticle*> tracks; - if ( ! jet.getAssociatedObjects(m_assocTracksName, tracks) ) { - ATH_MSG_DEBUG("Associated tracks not found."); - } + for(const xAOD::Jet* jet : jets) { + if ( ! jet->getAssociatedObjects(m_assocTracksName, tracks) ) { + ATH_MSG_DEBUG("Associated tracks not found."); + } - if (vertexContainer->size() == 0 ) { - ATH_MSG_WARNING("There are no vertices in the container. Exiting"); - return 4; - } + if (vertexContainer->size() == 0 ) { + ATH_MSG_WARNING("There are no vertices in the container. Exiting"); + return StatusCode::FAILURE; + } - const xAOD::Vertex* HSvertex = findHSVertex(vertexContainer); + const xAOD::Vertex* HSvertex = findHSVertex(vertexContainer); - const std::pair<float,float> tracksums = getJetTrackSums(HSvertex, tracks, tva); + const std::pair<float,float> tracksums = getJetTrackSums(HSvertex, tracks, tva); - jet.setAttribute("TrackSumPt",tracksums.first); - jet.setAttribute("TrackSumMass",tracksums.second); + trackSumPtHandle(*jet) = tracksums.first; + trackSumMassHandle(*jet) = tracksums.second; - ATH_MSG_VERBOSE("JetTrackSumMomentsTool " << name() - << ": TrackSumPt=" << tracksums.first - << ", TrackSumMass=" << tracksums.second ); - return 0; + ATH_MSG_VERBOSE("JetTrackSumMomentsTool " << name() + << ": TrackSumPt=" << tracksums.first + << ", TrackSumMass=" << tracksums.second ); + } + return StatusCode::SUCCESS; } //********************************************************************** diff --git a/Reconstruction/Jet/JetMomentTools/src/JetCaloCellQualityTool.cxx b/Reconstruction/Jet/JetMomentTools/src/JetCaloCellQualityTool.cxx index 2dd7e9957e6a676fc454bac9237fd8fbf5de59e6..710f2a886e5dcdcc2781b49aa41702deb0bfe615 100644 --- a/Reconstruction/Jet/JetMomentTools/src/JetCaloCellQualityTool.cxx +++ b/Reconstruction/Jet/JetMomentTools/src/JetCaloCellQualityTool.cxx @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration */ #ifndef XAOD_ANALYSIS @@ -9,15 +9,8 @@ #include <iostream> JetCaloCellQualityTool::JetCaloCellQualityTool(const std::string& name) - : JetCaloQualityTool(name) - , m_attSuffix("FromCells") - -{ - declareProperty("CellSuffix", m_attSuffix); - - declareProperty("LArQualityCut", m_LArQualityCut=4000); - declareProperty("TileQualityCut", m_TileQualityCut=254); - + : JetCaloQualityTool(name) { + declareInterface<IJetDecorator>(this); } StatusCode JetCaloCellQualityTool::initialize(){ diff --git a/Reconstruction/Jet/JetMomentTools/src/JetCaloCellQualityTool.h b/Reconstruction/Jet/JetMomentTools/src/JetCaloCellQualityTool.h index 33e3174ec1143054e82d65789bd3f50a8b143c92..bc1f52060a22e1b28dd7f3844cbf97b9000267a3 100644 --- a/Reconstruction/Jet/JetMomentTools/src/JetCaloCellQualityTool.h +++ b/Reconstruction/Jet/JetMomentTools/src/JetCaloCellQualityTool.h @@ -1,7 +1,7 @@ // this file is -*- C++ -*- /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration */ /** @@ -42,12 +42,11 @@ public: /// This objects holds a list of cell-based calculators jet::JetCaloCellCalculations m_cellCalculators; - /// Suffix added to the name of the jet attribute - std::string m_attSuffix; + Gaudi::Property<std::string> m_attSuffix{this, "CellSuffix", "FromCells", "Suffix added to the name of the jet attribute"}; // parameters for Quality cuts - int m_LArQualityCut; - int m_TileQualityCut; + Gaudi::Property<int> m_LArQualityCut{this, "LArQualityCut", 4000, "LAr quality cut parameter"}; + Gaudi::Property<int> m_TileQualityCut{this, "TileQualityCut", 254, "Tile quality cut parameter"}; }; #endif diff --git a/Reconstruction/Jet/JetMomentTools/src/JetIsolationTool.cxx b/Reconstruction/Jet/JetMomentTools/src/JetIsolationTool.cxx index 71f805d17f0599cb2c0c5d3e1f8a49de1fb8e15b..a6be32edc1ffaf6d67c315723c707c3316a346b4 100644 --- a/Reconstruction/Jet/JetMomentTools/src/JetIsolationTool.cxx +++ b/Reconstruction/Jet/JetMomentTools/src/JetIsolationTool.cxx @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration */ // JetIsolationTool.cxx @@ -60,18 +60,27 @@ namespace jet { virtual IsolationResult jetIsolation(const xAOD::Jet*, std::vector<jet::ParticlePosition> & ) const { return IsolationResult(); } - virtual void - setIsolationAttributes(xAOD::Jet* jet, std::vector<jet::ParticlePosition>& nearbyConstit) const { + virtual std::map<std::string, float> + calcIsolationAttributes(const xAOD::Jet* jet, std::vector<jet::ParticlePosition>& nearbyConstit) const { IsolationResult result = jetIsolation(jet, nearbyConstit); + std::map<std::string, float> result_map; for ( size_t i=0; i<m_kinematics.size(); ++i ) { switch( m_kinematics[i] ) { - case Perp : jet->setAttribute<float>( m_attNames[i], result.isoPerp(jet->p4().Vect()) ); break; - case SumPt : jet->setAttribute<float>( m_attNames[i], result.isoSumPt() ); break; - case Par : jet->setAttribute<float>( m_attNames[i], result.isoPar(jet->p4().Vect()) ); break; - case P : jet->setAttribute<float>( m_attNames[i], result.isoP().Vect().Mag() ); break; + case Perp: + result_map.insert(std::pair<std::string, float>(s_kname[i], result.isoPerp(jet->p4().Vect()) )); + break; + case SumPt: + result_map.insert(std::pair<std::string, float>(s_kname[i], result.isoSumPt() )); + break; + case Par: + result_map.insert(std::pair<std::string, float>(s_kname[i], result.isoPar(jet->p4().Vect()) )); + break; + case P: + result_map.insert(std::pair<std::string, float>(s_kname[i], result.isoP().Vect().Mag() )); + break; } } - + return result_map; } bool scheduleKinematicCalculation(string kname){ @@ -235,8 +244,8 @@ using namespace jet::JetIsolation; //********************************************************************** JetIsolationTool::JetIsolationTool(const string& name) -: JetModifierBase(name) { - declareProperty( "IsolationCalculations", m_isolationCodes); + : asg::AsgTool(name) { + declareInterface<IJetDecorator>(this); } //********************************************************************** @@ -272,11 +281,22 @@ StatusCode JetIsolationTool::initialize() { } } + if(m_jetContainerName.empty()) { + ATH_MSG_ERROR("JetIsolationTool needs to have its input jet container configured!"); + return StatusCode::FAILURE; + } + // Fill the iso calculator vector from the map + // Also fill DecorHandleKeyArrays at the same time for ( auto& pair : calcMap ){ m_isoCalculators.push_back(pair.second); ATH_MSG_DEBUG("Will use iso calculation : "<< pair.second->baseName() ); pair.second->dump(); + + m_perpKeys.emplace_back(m_jetContainerName + "." + pair.second->baseName() + "Perp"); + m_sumPtKeys.emplace_back(m_jetContainerName + "." + pair.second->baseName() + "SumPt"); + m_parKeys.emplace_back(m_jetContainerName + "." + pair.second->baseName() + "Par"); + m_pKeys.emplace_back(m_jetContainerName + "." + pair.second->baseName() + "P"); } ATH_MSG_INFO("Initialized JetIsolationTool " << name()); @@ -289,13 +309,24 @@ StatusCode JetIsolationTool::initialize() { } ATH_MSG_INFO(" Isolation calculations: " << m_isolationCodes); + m_inputTypeKey = m_jetContainerName + "." + m_inputTypeKey.key(); + ATH_CHECK(m_inputTypeKey.initialize()); + + ATH_CHECK(m_perpKeys.initialize()); + ATH_CHECK(m_sumPtKeys.initialize()); + ATH_CHECK(m_parKeys.initialize()); + ATH_CHECK(m_pKeys.initialize()); + + return StatusCode::SUCCESS; } //********************************************************************** -StatusCode JetIsolationTool::modify(xAOD::JetContainer& jets) const { +StatusCode JetIsolationTool::decorate(const xAOD::JetContainer& jets) const { + SG::ReadDecorHandle<xAOD::JetContainer, int> inputTypeHandle(m_inputTypeKey); + ATH_MSG_DEBUG("Modifying jets in container with size " << jets.size()); if ( jets.empty() ) return StatusCode::SUCCESS; @@ -314,10 +345,10 @@ StatusCode JetIsolationTool::modify(xAOD::JetContainer& jets) const { } // Loop over jets in this collection. - for ( xAOD::Jet* pjet : jets ) { + for (const xAOD::Jet* pjet : jets ) { // Check this jet has the same inputs. - // int jinp = pjet->getAttribute<int>("InputType"); + // int jinp = inputTypeHandle(*pjet); // This needs to be reimplemented when we decide how to better // encode this information -- right now this can't be matched // to the input PseudoJetContainer @@ -370,9 +401,19 @@ StatusCode JetIsolationTool::modify(xAOD::JetContainer& jets) const { << nearbyC.size() << ", in jet: "<< pjet->numConstituents() << ", total: "<< inputConstits->size() ); + std::map<std::string, const SG::WriteDecorHandleKeyArray<xAOD::JetContainer>* > decorKeyMap; + decorKeyMap.emplace(jet::JetIsolation::IsolationCalculator::s_kname[0], &m_sumPtKeys); + decorKeyMap.emplace(jet::JetIsolation::IsolationCalculator::s_kname[1], &m_parKeys); + decorKeyMap.emplace(jet::JetIsolation::IsolationCalculator::s_kname[2], &m_perpKeys); + decorKeyMap.emplace(jet::JetIsolation::IsolationCalculator::s_kname[3], &m_pKeys); + // loop over calculators, calculate isolation given the close-by particles not part of the jet. - for ( IsolationCalculator* calc : calculators ){ - calc->setIsolationAttributes(pjet, nearbyC); + for ( size_t iCalc = 0; iCalc < calculators.size(); iCalc++ ) { + std::map<std::string, float> results = calculators.at(iCalc)->calcIsolationAttributes(pjet, nearbyC); + for( const auto& [var_name, value]: results ) { + SG::WriteDecorHandle<xAOD::JetContainer, float> decorHandle(decorKeyMap.at(var_name)->at(iCalc)); + decorHandle(*pjet) = value; + } } } diff --git a/Reconstruction/Jet/JetMomentTools/src/components/JetMomentTools_entries.cxx b/Reconstruction/Jet/JetMomentTools/src/components/JetMomentTools_entries.cxx index 4b8bd610fc09d0912c42d18cf61c8221441a7bf7..bdc099dee4edcc37546ed51e81ac164db3d35105 100644 --- a/Reconstruction/Jet/JetMomentTools/src/components/JetMomentTools_entries.cxx +++ b/Reconstruction/Jet/JetMomentTools/src/components/JetMomentTools_entries.cxx @@ -9,7 +9,6 @@ #include "JetMomentTools/JetTrackSumMomentsTool.h" #include "JetMomentTools/JetClusterMomentsTool.h" #include "JetMomentTools/JetVoronoiMomentsTool.h" -#include "JetMomentTools/JetMuonSegmentMomentsTool.h" #include "JetMomentTools/JetPtAssociationTool.h" #include "JetMomentTools/JetIsolationTool.h" #include "JetMomentTools/JetLArHVTool.h" @@ -33,7 +32,6 @@ DECLARE_COMPONENT( JetTrackMomentsTool ) DECLARE_COMPONENT( JetTrackSumMomentsTool ) DECLARE_COMPONENT( JetClusterMomentsTool ) DECLARE_COMPONENT( JetVoronoiMomentsTool ) -DECLARE_COMPONENT( JetMuonSegmentMomentsTool ) DECLARE_COMPONENT( JetPtAssociationTool ) DECLARE_COMPONENT( JetIsolationTool ) DECLARE_COMPONENT( JetLArHVTool )