diff --git a/Reconstruction/MET/METUtilities/METUtilities/METMaker.h b/Reconstruction/MET/METUtilities/METUtilities/METMaker.h index c749d6a2d914654f32343a9f76a9649a20cce1e7..94844e962b35e64473c81fb59b94df3c0c2ffd80 100644 --- a/Reconstruction/MET/METUtilities/METUtilities/METMaker.h +++ b/Reconstruction/MET/METUtilities/METUtilities/METMaker.h @@ -4,10 +4,10 @@ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration */ -// METMaker.h +// METMaker.h // Header file for class METMaker // Author: T.J.Khoo<khoo@cern.ch> -/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// #ifndef METUTILITIES_MET_METMAKER_H #define METUTILITIES_MET_METMAKER_H 1 @@ -27,113 +27,165 @@ // Forward declaration namespace met { - + // typedefs typedef ElementLink<xAOD::IParticleContainer> obj_link_t; - + class METMaker - : virtual public asg::AsgTool, - virtual public IMETMaker - - { + : virtual public asg::AsgTool, + virtual public IMETMaker + + { // This macro defines the constructor with the interface declaration ASG_TOOL_CLASS(METMaker, IMETMaker) - - /////////////////////////////////////////////////////////////////// - // Public methods: - /////////////////////////////////////////////////////////////////// - public: - - // Copy constructor: - - /// Constructor with parameters: + + /////////////////////////////////////////////////////////////////// + // Public methods: + /////////////////////////////////////////////////////////////////// + public: + + // Copy constructor: + + /// Constructor with parameters: METMaker(const std::string& name); - - /// Destructor: - virtual ~METMaker(); - + + /// Destructor: + virtual ~METMaker(); + // Athena algtool's Hooks - virtual StatusCode initialize(); - virtual StatusCode finalize(); - // virtual StatusCode execute(); - - virtual StatusCode rebuildMET(xAOD::MissingET* met, - const xAOD::IParticleContainer* collection, - const xAOD::MissingETAssociationMap* map, - std::vector<const xAOD::IParticle*>& uniques); - virtual StatusCode rebuildMET(xAOD::MissingET* met, - const xAOD::IParticleContainer* collection, - const xAOD::MissingETAssociationMap* map, - std::vector<const xAOD::IParticle*>& uniques, - MissingETBase::UsageHandler::Policy p, - bool removeOverlap); - - inline virtual StatusCode rebuildMET(xAOD::MissingET* met, - const xAOD::IParticleContainer* collection, - const xAOD::MissingETAssociationMap* map) - { - std::vector<const xAOD::IParticle*> uniques; - return rebuildMET(met, collection, map, uniques); - } - - virtual StatusCode rebuildJetMET(xAOD::MissingET* metJet, - const xAOD::JetContainer* jets, - const xAOD::MissingETAssociationMap* map, - std::vector<const xAOD::IParticle*>& uniques, - xAOD::MissingET* metSoftClus, - const xAOD::MissingET* coreSoftClus, - xAOD::MissingET* metSoftTrk, - const xAOD::MissingET* coreSoftTrk); - - inline virtual StatusCode rebuildJetMET(xAOD::MissingET* metJet, - const xAOD::JetContainer* jets, - const xAOD::MissingETAssociationMap* map, - xAOD::MissingET* metSoftClus, - const xAOD::MissingET* coreSoftClus, - xAOD::MissingET* metSoftTrk, - const xAOD::MissingET* coreSoftTrk) - { - std::vector<const xAOD::IParticle*> uniques; - return rebuildJetMET(metJet, jets, map, uniques, - metSoftClus, coreSoftClus, - metSoftTrk, coreSoftTrk); - } - - - virtual StatusCode buildMETSum(const std::string& totalName, - xAOD::MissingETContainer* metCont, - MissingETBase::Types::bitmask_t softTermsSource=0); - - /////////////////////////////////////////////////////////////////// - // Const methods: + StatusCode initialize(); + StatusCode finalize(); + + StatusCode rebuildMET(const std::string& metKey, + xAOD::Type::ObjectType metType, + xAOD::MissingETContainer* metCont, + const xAOD::IParticleContainer* collection, + const xAOD::MissingETAssociationMap* map, + std::vector<const xAOD::IParticle*>& uniques); + // + StatusCode rebuildMET(const std::string& metKey, + xAOD::Type::ObjectType metType, + xAOD::MissingETContainer* metCont, + const xAOD::IParticleContainer* collection, + const xAOD::MissingETAssociationMap* map); + // + StatusCode rebuildMET(xAOD::MissingET* met, + const xAOD::IParticleContainer* collection, + const xAOD::MissingETAssociationMap* map); + // + StatusCode rebuildMET(xAOD::MissingET* met, + const xAOD::IParticleContainer* collection, + const xAOD::MissingETAssociationMap* map, + std::vector<const xAOD::IParticle*>& uniques); + // + StatusCode rebuildMET(xAOD::MissingET* met, + const xAOD::IParticleContainer* collection, + const xAOD::MissingETAssociationMap* map, + std::vector<const xAOD::IParticle*>& uniques, + MissingETBase::UsageHandler::Policy p, + bool removeOverlap); + + StatusCode rebuildJetMET(const std::string& metJetKey, + const std::string& softClusKey, + const std::string& softTrkKey, + xAOD::MissingETContainer* metCont, + const xAOD::JetContainer* jets, + const xAOD::MissingETContainer* metCoreCont, + const xAOD::MissingETAssociationMap* map, + bool doJetJVF); + StatusCode rebuildJetMET(const std::string& metJetKey, + const std::string& softClusKey, + const std::string& softTrkKey, + xAOD::MissingETContainer* metCont, + const xAOD::JetContainer* jets, + const xAOD::MissingETContainer* metCoreCont, + const xAOD::MissingETAssociationMap* map, + bool doJetJVF, + std::vector<const xAOD::IParticle*>& uniques); + // + StatusCode rebuildJetMET(const std::string& metJetKey, + const std::string& metSoftKey, + xAOD::MissingETContainer* metCont, + const xAOD::JetContainer* jets, + const xAOD::MissingETContainer* metCoreCont, + const xAOD::MissingETAssociationMap* map, + bool doJetJVF); + StatusCode rebuildJetMET(const std::string& metJetKey, + const std::string& metSoftKey, + xAOD::MissingETContainer* metCont, + const xAOD::JetContainer* jets, + const xAOD::MissingETContainer* metCoreCont, + const xAOD::MissingETAssociationMap* map, + bool doJetJVF, + std::vector<const xAOD::IParticle*>& uniques); + // + StatusCode rebuildJetMET(xAOD::MissingET* metJet, + const xAOD::JetContainer* jets, + const xAOD::MissingETAssociationMap* map, + xAOD::MissingET* metSoftClus, + const xAOD::MissingET* coreSoftClus, + xAOD::MissingET* metSoftTrk, + const xAOD::MissingET* coreSoftTrk, + bool doJetJVF, + std::vector<const xAOD::IParticle*>& uniques, + bool tracksForHardJets = 0); + + StatusCode rebuildTrackMET(const std::string& metJetKey, + const std::string& softTrkKey, + xAOD::MissingETContainer* metCont, + const xAOD::JetContainer* jets, + const xAOD::MissingETContainer* metCoreCont, + const xAOD::MissingETAssociationMap* map, + bool doJetJVF); + StatusCode rebuildTrackMET(const std::string& metJetKey, + const std::string& softTrkKey, + xAOD::MissingETContainer* metCont, + const xAOD::JetContainer* jets, + const xAOD::MissingETContainer* metCoreCont, + const xAOD::MissingETAssociationMap* map, + bool doJetJVF, + std::vector<const xAOD::IParticle*>& uniques); + StatusCode rebuildTrackMET(xAOD::MissingET* metJet, + const xAOD::JetContainer* jets, + const xAOD::MissingETAssociationMap* map, + xAOD::MissingET* metSoftTrk, + const xAOD::MissingET* coreSoftTrk, + bool doJetJVF, + std::vector<const xAOD::IParticle*>& uniques); + + StatusCode buildMETSum(const std::string& totalName, + xAOD::MissingETContainer* metCont, + MissingETBase::Types::bitmask_t softTermsSource=0); + + /////////////////////////////////////////////////////////////////// + // Const methods: + /////////////////////////////////////////////////////////////////// + + /////////////////////////////////////////////////////////////////// + // Non-const methods: + /////////////////////////////////////////////////////////////////// + + /////////////////////////////////////////////////////////////////// + // Private data: /////////////////////////////////////////////////////////////////// - - /////////////////////////////////////////////////////////////////// - // Non-const methods: - /////////////////////////////////////////////////////////////////// - - /////////////////////////////////////////////////////////////////// - // Private data: - /////////////////////////////////////////////////////////////////// private: - - bool m_doJetJvf; + bool m_jetCorrectPhi; - + double m_jetMinPt; double m_jetMinAbsJvf; double m_jetMinEfrac; double m_jetMinWeightedPt; - + bool m_doPFlow; - + // Set up accessors to original object links in case of corrected copy containers SG::AuxElement::ConstAccessor<obj_link_t> m_objLinkAcc; - - /// Default constructor: + + /// Default constructor: METMaker(); - - }; - + + }; + } //> end namespace met #endif //> !METUTILITIES_MET_METMAKER_H diff --git a/Reconstruction/MET/METUtilities/METUtilities/METRebuilder.h b/Reconstruction/MET/METUtilities/METUtilities/METRebuilder.h index e43d5f2aa5c7e05f1a1dda552844e57386ad938f..b93fa5640fd9dae713d81955c2026ec9759e1500 100644 --- a/Reconstruction/MET/METUtilities/METUtilities/METRebuilder.h +++ b/Reconstruction/MET/METUtilities/METUtilities/METRebuilder.h @@ -4,10 +4,10 @@ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration */ -// METRebuilder.h +// METRebuilder.h // Header file for class METRebuilder // Author: T.J.Khoo<khoo@cern.ch> -/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// #ifndef METUTILITIES_MET_METREBUILDER_H #define METUTILITIES_MET_METREBUILDER_H 1 @@ -15,95 +15,118 @@ #include <string> // FrameWork includes -#include "AsgTools/ToolHandle.h" #include "AsgTools/AsgTool.h" +#include "AsgTools/ToolHandle.h" // METInterface includes #include "METInterface/IMETRebuilder.h" +// Tracking Tool +#include "InDetTrackSelectionTool/IInDetTrackSelectionTool.h" + // EDM includes #include "xAODJet/JetContainer.h" - #include "xAODTracking/VertexFwd.h" #include "xAODTracking/TrackParticleFwd.h" // Forward declaration namespace met { - + // typedefs typedef ElementLink<xAOD::IParticleContainer> obj_link_t; - + class METRebuilder - : virtual public asg::AsgTool, - virtual public IMETRebuilder - - { + : virtual public asg::AsgTool, + virtual public IMETRebuilder + + { // This macro defines the constructor with the interface declaration ASG_TOOL_CLASS(METRebuilder, IMETRebuilder) - - /////////////////////////////////////////////////////////////////// - // Public methods: - /////////////////////////////////////////////////////////////////// - public: - - // Copy constructor: - - /// Constructor with parameters: + + /////////////////////////////////////////////////////////////////// + // Public methods: + /////////////////////////////////////////////////////////////////// + public: + + // Copy constructor: + + /// Constructor with parameters: METRebuilder(const std::string& name); - - /// Destructor: - virtual ~METRebuilder(); - + + /// Destructor: + virtual ~METRebuilder(); + // Athena algtool's Hooks - virtual StatusCode initialize(); - virtual StatusCode finalize(); - virtual StatusCode execute(); - - virtual StatusCode rebuildMET(const std::string& metKey, - xAOD::MissingETContainer* metCont, - const xAOD::IParticleContainer* collection, - const xAOD::MissingETComponentMap* metMap, - bool doTracks=true); - virtual StatusCode rebuildMET(xAOD::MissingET* met, - const xAOD::IParticleContainer* collection, - const xAOD::MissingETComponent* component, - bool doTracks); - - virtual StatusCode rebuildJetMET(const std::string& jetKey, - const std::string& softKey, - xAOD::MissingETContainer* metCont, - const xAOD::JetContainer* jets, - const xAOD::MissingETComponentMap* metMap, - bool doTracks=true); - virtual StatusCode rebuildJetMET(xAOD::MissingET* metJet, - xAOD::MissingET* metSoft, - const xAOD::MissingET* oldSoft, - const xAOD::JetContainer* jets, - const xAOD::MissingETComponent* component, - bool doTracks); - virtual StatusCode buildMETSum(const std::string& totalName, - xAOD::MissingETContainer* metCont); - - /////////////////////////////////////////////////////////////////// - // Const methods: + StatusCode initialize(); + StatusCode finalize(); + StatusCode execute(); + + StatusCode copyMET(const std::string& metKey, + xAOD::MissingETContainer* metCont, + const xAOD::MissingETComponentMap* metMap); + + StatusCode rebuildMET(const std::string& metKey, + xAOD::MissingETContainer* metCont, + const xAOD::IParticleContainer* collection, + const xAOD::MissingETComponentMap* metMap, + bool doTracks=true); + + StatusCode rebuildMET(xAOD::MissingET* met, + const xAOD::IParticleContainer* collection, + const xAOD::MissingETComponent* component, + bool doTracks=true); + + StatusCode rebuildJetMET(const std::string& jetKey, + const std::string& softKey, + xAOD::MissingETContainer* metCont, + const xAOD::JetContainer* jets, + const xAOD::MissingETComponentMap* metMap, + bool doTracks=true) { + return rebuildJetMET(jetKey,softKey,metCont,jets,metMap,doTracks, + m_jetDoJvf,m_pureTrkSoft,m_softJetScale); + } + + StatusCode rebuildJetMET(const std::string& jetKey, + const std::string& softKey, + xAOD::MissingETContainer* metCont, + const xAOD::JetContainer* jets, + const xAOD::MissingETComponentMap* metMap, + bool doTracks, + bool doJvfCut, + bool pureTrkSoft, + const std::string& jetScale); + + StatusCode rebuildJetMET(xAOD::MissingET* metJet, + xAOD::MissingET* metSoft, + const xAOD::JetContainer* jets, + const xAOD::MissingETComponent* component, + bool doTracks, + bool doJvfCut, + bool pureTrkSoft, + const std::string& jetScale); + StatusCode buildMETSum(const std::string& totalName, + xAOD::MissingETContainer* metCont); + /////////////////////////////////////////////////////////////////// - - /////////////////////////////////////////////////////////////////// - // Non-const methods: - /////////////////////////////////////////////////////////////////// - - /////////////////////////////////////////////////////////////////// - // Private data: - /////////////////////////////////////////////////////////////////// - private: - bool isPVTrack(const xAOD::TrackParticle* trk, - const xAOD::Vertex* pv) const; + // Const methods: + /////////////////////////////////////////////////////////////////// + + /////////////////////////////////////////////////////////////////// + // Non-const methods: + /////////////////////////////////////////////////////////////////// + + /////////////////////////////////////////////////////////////////// + // Private data: + /////////////////////////////////////////////////////////////////// + private: + bool acceptTrack(const xAOD::TrackParticle* trk, + const xAOD::Vertex* pv) const; void associateTracks(const xAOD::IParticle* obj); - - /// Default constructor: + + /// Default constructor: METRebuilder(); - + std::string m_eleColl; std::string m_gammaColl; std::string m_tauColl; @@ -116,26 +139,23 @@ namespace met { std::string m_jetTerm; std::string m_muonTerm; std::string m_softTerm; + std::string m_softTermType; // - std::string m_inputMap; + std::string m_inputMap; std::string m_outMETCont; std::string m_outMETTerm; bool m_warnOfDupes; - + bool m_doEle; bool m_doGamma; bool m_doTau; bool m_doMuon; - + bool m_rebuildEle; bool m_rebuildGamma; bool m_rebuildTau; bool m_rebuildMuon; - - // Set up accessors to original object links in case of corrected copy containers - SG::AuxElement::Accessor<obj_link_t> m_objLinkAcc; - SG::AuxElement::Decorator<char> m_trkUsedDec; - + // For jet/soft term -- eventually break off into a separate tool double m_jetPtCut; bool m_jetDoJvf; @@ -143,14 +163,15 @@ namespace met { std::string m_softJetScale; bool m_doTracks; bool m_pureTrkSoft; - - // Should be in a track selector + + // Decorate tracks to state that they have been used for a MET calc + SG::AuxElement::Decorator<char> m_trkUsedDec; + bool m_trk_doPVsel; + ToolHandle<InDet::IInDetTrackSelectionTool> m_trkseltool; std::string m_vtxColl; std::string m_clusColl; - double m_trk_d0Max; - double m_trk_z0Max; }; - + } //> end namespace met #endif //> !METUTILITIES_MET_METREBUILDER_H diff --git a/Reconstruction/MET/METUtilities/METUtilities/METSystematicsTool.h b/Reconstruction/MET/METUtilities/METUtilities/METSystematicsTool.h new file mode 100644 index 0000000000000000000000000000000000000000..7d2cc1b0e0337e573c55df89122216b945ce4bc0 --- /dev/null +++ b/Reconstruction/MET/METUtilities/METUtilities/METSystematicsTool.h @@ -0,0 +1,167 @@ +///////////////////////// -*- C++ -*- ///////////////////////////// + +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +// METSystematicsTool.h +// Header file for class METSystematicsTool +// Author: Russell Smith <rsmith@cern.ch> +/////////////////////////////////////////////////////////////////// +#ifndef METUTILITIES_METSYSTEMATICSTOOL_H +#define METUTILITIES_METSYSTEMATICSTOOL_H + +class TH3D; +class TH2D; +class TH1D; + +#include "METInterface/IMETSystematicsTool.h" +#include "AsgTools/AsgTool.h" + +#include "PATInterfaces/SystematicsTool.h" + +#include "TRandom3.h" + +#include "xAODMissingET/MissingETContainer.h" +#include "xAODMissingET/MissingETAssociationMap.h" + +#include "xAODJet/Jet.h" + +namespace met { + + namespace softCaloAffSyst { + const static CP::SystematicVariation MET_SoftCalo_ScaleUp ("MET_SoftCalo_ScaleUp" ); + const static CP::SystematicVariation MET_SoftCalo_ScaleDown("MET_SoftCalo_ScaleDown"); + const static CP::SystematicVariation MET_SoftCalo_Reso ("MET_SoftCalo_Reso"); + } + + namespace softTrkAffSyst { + const static CP::SystematicVariation MET_SoftTrk_ScaleUp ("MET_SoftTrk_ScaleUp" ); + const static CP::SystematicVariation MET_SoftTrk_ScaleDown("MET_SoftTrk_ScaleDown"); + const static CP::SystematicVariation MET_SoftTrk_ResoPara ("MET_SoftTrk_ResoPara" ); + const static CP::SystematicVariation MET_SoftTrk_ResoPerp ("MET_SoftTrk_ResoPerp" ); + const static CP::SystematicVariation MET_SoftTrk_ResoCorr ("MET_SoftTrk_ResoCorr" ); + } + + namespace jetTrkAffSyst { + const static CP::SystematicVariation MET_JetTrk_ScaleUp ("MET_JetTrk_ScaleUp" ); + const static CP::SystematicVariation MET_JetTrk_ScaleDown ("MET_JetTrk_ScaleDown"); + } + + enum SystType { + SOFTCALO, + SOFTTRK, + JETTRK + }; + + //this is to speed up applyCorrection lookup. + //we still store the full list in appliedSystematics()/m_appliedSystematics + enum SystApplied { + NONE , + MET_SOFTTRK_SCALEUP , + MET_SOFTTRK_SCALEDOWN , + MET_SOFTTRK_RESOPARA , + MET_SOFTTRK_RESOPERP , + MET_SOFTTRK_RESOCORR , + MET_SOFTCALO_SCALEUP , + MET_SOFTCALO_SCALEDOWN, + MET_SOFTCALO_RESO , + MET_JETTRK_SCALEUP , + MET_JETTRK_SCALEDOWN + }; + + class METSystematicsTool : public virtual IMETSystematicsTool, + public virtual CP::SystematicsTool, + public virtual asg::AsgTool + { + // This macro defines the constructor with the interface declaration + ASG_TOOL_CLASS(METSystematicsTool, IMETSystematicsTool) + + public: + + xAOD::EventInfo const * getDefaultEventInfo() const; + + //Constructor + METSystematicsTool(const std::string& name); + + //Destructor + virtual ~METSystematicsTool(){} + + + StatusCode softTrkSystInitialize(); //initialize softTrk scale/reo histos from config file + StatusCode softCaloSystInitialize();//initialize softCalo scale/reso histos from config file + StatusCode jetTrkSystInitialize(); //initialize jetTrk scale/reso histos from config file + + //required asg tool interface functions + StatusCode initialize(); + StatusCode finalize(); + + //required correction tool functions + //we don't inherit from CorrectionTool directly, since we don't want to implement applyContainerCorrection + CP::CorrectionCode applyCorrection(xAOD::MissingET& inputMet, xAOD::MissingETAssociationMap * map = nullptr ) const; + CP::CorrectionCode correctedCopy(const xAOD::MissingET& met, xAOD::MissingET*& outputmet,xAOD::MissingETAssociationMap * map = nullptr) const; + //We don't want these for MET since we only apply systematics to the soft term, and this may be unclear + //virtual CP::CorrectionCode applyContainerCorrection(xAOD::MissingETContainer& inputs, const CP::xAOD::EventInfo& eInfo) const; + //virtual CP::CorrectionCode applyContainerCorrection(xAOD::MissingETContainer& inputs, const CP::xAOD::EventInfo& eInfo) const; + + + //required systematic tool functions + bool isAffectedBySystematic (const CP::SystematicVariation& var) const{return CP::SystematicsTool::isAffectedBySystematic(var) ;} + CP::SystematicSet affectingSystematics () const{ return CP::SystematicsTool::affectingSystematics () ;} + CP::SystematicSet recommendedSystematics () const{ return CP::SystematicsTool::recommendedSystematics() ;} + CP::SystematicCode applySystematicVariation(const CP::SystematicSet& set){ return CP::SystematicsTool::applySystematicVariation(set) ;} + CP::SystematicCode sysApplySystematicVariation(const CP::SystematicSet&); //when inheriting from SystematicsTool, we should only have to implement this one + + void setRandomSeed(int seed); + + private: + SystApplied m_appliedSystEnum; + + CP::CorrectionCode internalSoftTermApplyCorrection(xAOD::MissingET& softMet, xAOD::MissingETContainer const * METcont , xAOD::EventInfo const& eInfo +) const; + CP::CorrectionCode calcJetTrackMETWithSyst(xAOD::MissingET& jettrkmet, const xAOD::MissingETAssociationMap* map, const xAOD::Jet* jet) const; + CP::CorrectionCode getCorrectedJetTrackMET(xAOD::MissingET& jettrkmet, const xAOD::MissingETAssociationMap* map, const xAOD::JetContainer* vecJets) const; + + + + //declared properties + std::string m_jetColl; + std::string m_configPrefix; + std::string m_configSoftTrkFile; + std::string m_configJetTrkFile; + std::string m_configSoftCaloFile; + std::string m_truthCont; + std::string m_truthObj; + std::string m_vertexCont; + std::string m_eventInfo; + int m_randSeed ; + + TH3D* shiftpara_pthard_njet_mu; + TH3D* resopara_pthard_njet_mu; + TH3D* resoperp_pthard_njet_mu; + TH2D* jet_systRpt_pt_eta; + TH1D* h_calosyst_scale; + TH1D* h_calosyst_reso; + + mutable TRandom3 m_rand;//so that we can call this from applyCorrection + + int getNPV() const; + + StatusCode addMETAffectingSystematics(); + StatusCode extractHistoPath(std::string & histfile, std::string & systpath, std::string & configdir, std::string & suffix, SystType const & type); + + xAOD::MissingET calcPtHard(xAOD::MissingETContainer const * const cont) const; + xAOD::MissingET caloSyst_scale(xAOD::MissingET const &softTerms , double const scale) const; + xAOD::MissingET caloSyst_reso (xAOD::MissingET const &softTerms) const; + xAOD::MissingET softTrkSyst_scale (xAOD::MissingET const & softTerms, xAOD::MissingET const & ptHard, double const shift) const; + xAOD::MissingET softTrkSyst_reso (xAOD::MissingET const & softTerms, xAOD::MissingET const & ptHard, double const shift, double const smearpara, + double const smearperp) const; + xAOD::MissingET projectST (xAOD::MissingET const & softTerms, xAOD::MissingET const & ptHard) const; + + };//METSystematicsTool + +} + +#endif //METUTILIES_METSYSTEMATICSTOOL_H + +// LocalWords: METSystematicsTool diff --git a/Reconstruction/MET/METUtilities/METUtilities/METUtilitiesDict.h b/Reconstruction/MET/METUtilities/METUtilities/METUtilitiesDict.h new file mode 100644 index 0000000000000000000000000000000000000000..ac5ee294a5b1e216db142100fb381d92351cb4c4 --- /dev/null +++ b/Reconstruction/MET/METUtilities/METUtilities/METUtilitiesDict.h @@ -0,0 +1,16 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef METUTILITIESDICT_H +#define METUTILITIESDICT_H + +#if defined(__GCCXML__) and not defined(EIGEN_DONT_VECTORIZE) +#define EIGEN_DONT_VECTORIZE +#endif // __GCCXML__ + +#include "METUtilities/METRebuilder.h" +#include "METUtilities/METMaker.h" +#include "METUtilities/METSystematicsTool.h" + +#endif //METUTILITIESDICT_H diff --git a/Reconstruction/MET/METUtilities/METUtilities/selection.xml b/Reconstruction/MET/METUtilities/METUtilities/selection.xml new file mode 100644 index 0000000000000000000000000000000000000000..61000da12c4c79b939728d7a3d5faadf5346b00c --- /dev/null +++ b/Reconstruction/MET/METUtilities/METUtilities/selection.xml @@ -0,0 +1,12 @@ +<lcgdict> + <class name="met::METRebuilder" /> + <class name="met::METMaker" /> + <class name="met::METSystematicsTool" /> + + <!-- Suppress the unwanted classes found by ROOT 6. --> + <!-- Hopefully we can remove these extra lines at one point... --> + <exclusion> + <class name="SG::IConstAuxStore" /> + <class name="DataLink<SG::IConstAuxStore>" /> + </exclusion> +</lcgdict> diff --git a/Reconstruction/MET/METUtilities/Root/METMaker.cxx b/Reconstruction/MET/METUtilities/Root/METMaker.cxx index 33d2f71f6ef66872b8818bf09aa6744e4be4bd5b..6849c7d85c0884abf6b3a770859954f19d41e2d7 100644 --- a/Reconstruction/MET/METUtilities/Root/METMaker.cxx +++ b/Reconstruction/MET/METUtilities/Root/METMaker.cxx @@ -4,10 +4,10 @@ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration */ -// METMaker.cxx +// METMaker.cxx // Implementation file for class METMaker // Author: T.J.Khoo<khoo@cern.ch> -/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// // METUtilities includes #include "METUtilities/METMaker.h" @@ -25,9 +25,9 @@ #include "xAODTracking/TrackParticle.h" namespace met { - + using std::vector; - + using xAOD::MissingET; using xAOD::MissingETContainer; using xAOD::MissingETAssociation; @@ -42,64 +42,115 @@ namespace met { using xAOD::JetConstituentVector; // using xAOD::TrackParticle; - - /////////////////////////////////////////////////////////////////// - // Public methods: - /////////////////////////////////////////////////////////////////// - + + /////////////////////////////////////////////////////////////////// + // Public methods: + /////////////////////////////////////////////////////////////////// + // Constructors //////////////// - METMaker::METMaker(const std::string& name) : - AsgTool(name), - m_objLinkAcc("originalObjectLink") + METMaker::METMaker(const std::string& name) : + AsgTool(name), + m_objLinkAcc("originalObjectLink") { // // Property declaration - // - + // + declareProperty("JetPtCut", m_jetMinPt = 20e3 ); - declareProperty("DoJetJvfCut", m_doJetJvf = true ); declareProperty("JetJvfCut", m_jetMinAbsJvf = 0.25 ); declareProperty("JetMinEFrac", m_jetMinEfrac = 0.5 ); declareProperty("JetMinWeightedPt", m_jetMinWeightedPt = 0. ); declareProperty("CorrectJetPhi", m_jetCorrectPhi = false ); declareProperty("DoPFlow", m_doPFlow = false ); } - + // Destructor /////////////// METMaker::~METMaker() {} - + // Athena algtool's Hooks //////////////////////////// StatusCode METMaker::initialize() { ATH_MSG_INFO ("Initializing " << name() << "..."); - + return StatusCode::SUCCESS; } - + StatusCode METMaker::finalize() { ATH_MSG_INFO ("Finalizing " << name() << "..."); - + return StatusCode::SUCCESS; } - + // StatusCode METMaker::execute() // { // ATH_MSG_DEBUG ( name() << " in execute..."); - + // return StatusCode::SUCCESS; // } - + // **** Rebuild generic MET term **** + + inline StatusCode METMaker::rebuildMET(xAOD::MissingET* met, + const xAOD::IParticleContainer* collection, + const xAOD::MissingETAssociationMap* map) { + + std::vector<const xAOD::IParticle*> uniques; + return rebuildMET(met,collection,map,uniques); + + } + inline StatusCode METMaker::rebuildMET(const std::string& metKey, + xAOD::Type::ObjectType metType, + xAOD::MissingETContainer* metCont, + const xAOD::IParticleContainer* collection, + const xAOD::MissingETAssociationMap* map) + { + std::vector<const xAOD::IParticle*> uniques; + return rebuildMET(metKey,metType,metCont,collection,map,uniques); + } + + StatusCode METMaker::rebuildMET(const std::string& metKey, + xAOD::Type::ObjectType metType, + xAOD::MissingETContainer* metCont, + const xAOD::IParticleContainer* collection, + const xAOD::MissingETAssociationMap* map, + std::vector<const xAOD::IParticle*>& uniques) + { + MissingETBase::Types::bitmask_t metSource; + switch(metType) { + case xAOD::Type::Electron: + metSource = MissingETBase::Source::electron(); + break; + case xAOD::Type::Photon: + metSource = MissingETBase::Source::photon(); + break; + case xAOD::Type::Tau: + metSource = MissingETBase::Source::tau(); + break; + case xAOD::Type::Muon: + metSource = MissingETBase::Source::muon(); + break; + case xAOD::Type::Jet: + ATH_MSG_WARNING("Incorrect use of rebuildMET -- use rebuildJetMET for RefJet term"); + return StatusCode::FAILURE; + default: + ATH_MSG_WARNING("Invalid object type provided: " << metType); + return StatusCode::FAILURE; + } + MissingET* met = new MissingET(metKey,metSource); + metCont->push_back(met); + return rebuildMET(met,collection,map,uniques); + } + StatusCode METMaker::rebuildMET(xAOD::MissingET* met, - const xAOD::IParticleContainer* collection, - const xAOD::MissingETAssociationMap* map, - std::vector<const xAOD::IParticle*>& uniques) + const xAOD::IParticleContainer* collection, + const xAOD::MissingETAssociationMap* map, + std::vector<const xAOD::IParticle*>& uniques) { MissingETBase::UsageHandler::Policy p = MissingETBase::UsageHandler::OnlyCluster; bool removeOverlap = true; @@ -109,100 +160,244 @@ namespace met { if(!collection->empty()) { const IParticle* obj = collection->front(); if(obj->type()==xAOD::Type::Muon) { - if(!m_doPFlow) { - p = MissingETBase::UsageHandler::OnlyTrack; - } - removeOverlap = false; + if(!m_doPFlow) { + p = MissingETBase::UsageHandler::OnlyTrack; + } + removeOverlap = false; } } return rebuildMET(met,collection,map,uniques,p,removeOverlap); } - + StatusCode METMaker::rebuildMET(xAOD::MissingET* met, - const xAOD::IParticleContainer* collection, - const xAOD::MissingETAssociationMap* map, - std::vector<const xAOD::IParticle*>& uniques, - MissingETBase::UsageHandler::Policy p, - bool removeOverlap) { + const xAOD::IParticleContainer* collection, + const xAOD::MissingETAssociationMap* map, + std::vector<const xAOD::IParticle*>& uniques, + MissingETBase::UsageHandler::Policy p, + bool removeOverlap) { if(!met || !collection || !map) { ATH_MSG_WARNING("Invalid pointer supplied for " - << "MET (" << met << "), " - << "collection (" << collection << "), " - << " or map (" << map << ")."); + << "MET (" << met << "), " + << "collection (" << collection << "), " + << " or map (" << map << ")."); return StatusCode::SUCCESS; } - ATH_MSG_VERBOSE("Rebuilding MET term " << met->name()); + ATH_MSG_VERBOSE("Building MET term " << met->name()); uniques.clear(); if(!collection->empty()) { bool originalInputs = !collection->isAvailable<ElementLink<IParticleContainer> >("originalObjectLink"); for(const auto& obj : *collection) { - const MissingETAssociation* assoc = 0; - if(originalInputs) { - assoc = MissingETComposition::getAssociation(map,obj); - } else { - const IParticle* orig = *m_objLinkAcc(*obj); - assoc = MissingETComposition::getAssociation(map,orig); - } - if(assoc) { - bool selected(false); - // Don't overlap remove muons, but flag the non-overlapping muons to take out their tracks from jets - selected = !assoc->hasOverlaps(obj,p); - if(selected || !removeOverlap) {*met += obj;} - ATH_MSG_VERBOSE(obj->type() << " with pt " << obj->pt() - << " is " << ( selected ? "non-" : "") << "overlapping"); - assoc->setObjSelectionFlag(obj,selected); - if(selected) { - uniques.push_back(obj); - } - } else { - // special case for muons, which might not be associated to a jet (isolated) - if(obj->type()==xAOD::Type::Muon) { - *met += obj; - } - } + std::vector<const MissingETAssociation*> assoc; + if(originalInputs) { + assoc = MissingETComposition::getAssociations(map,obj); + } else { + const IParticle* orig = *m_objLinkAcc(*obj); + assoc = MissingETComposition::getAssociations(map,orig); + } + bool selected(true); + // Don't overlap remove muons, but flag the non-overlapping muons to take out their tracks from jets + for(size_t i = 0; i < assoc.size(); i++) selected = selected && !assoc[i]->hasOverlaps(obj,p); + if (map->getMiscAssociation()->containsPhysics(obj)) selected = selected && !map->getMiscAssociation()->hasOverlaps(obj,p); + if(selected || !removeOverlap) {*met += obj;} + ATH_MSG_VERBOSE(obj->type() << " with pt " << obj->pt() + << " is " << ( selected ? "non-" : "") << "overlapping"); + for(size_t i = 0; i < assoc.size(); i++) assoc[i]->setObjSelectionFlag(obj,selected); + if (map->getMiscAssociation()->containsPhysics(obj)) map->getMiscAssociation()->setObjSelectionFlag(obj,selected); + if(selected) { + uniques.push_back(obj); + } } } return StatusCode::SUCCESS; } + + // **** Rebuild jet & soft MET terms **** + inline StatusCode METMaker::rebuildJetMET(const std::string& metJetKey, + const std::string& softKey, + xAOD::MissingETContainer* metCont, + const xAOD::JetContainer* jets, + const xAOD::MissingETContainer* metCoreCont, + const xAOD::MissingETAssociationMap* map, + bool doJetJVF) + { + std::vector<const xAOD::IParticle*> uniques; + + return rebuildJetMET(metJetKey,softKey,metCont, + jets,metCoreCont,map,doJetJVF,uniques); + } // **** Rebuild jet & soft MET terms **** + inline StatusCode METMaker::rebuildTrackMET(const std::string& metJetKey, + const std::string& softKey, + xAOD::MissingETContainer* metCont, + const xAOD::JetContainer* jets, + const xAOD::MissingETContainer* metCoreCont, + const xAOD::MissingETAssociationMap* map, + bool doJetJVF) + { + std::vector<const xAOD::IParticle*> uniques; + + return rebuildTrackMET(metJetKey,softKey,metCont, + jets,metCoreCont,map,doJetJVF,uniques); + } + + StatusCode METMaker::rebuildJetMET(const std::string& metJetKey, + const std::string& softKey, + xAOD::MissingETContainer* metCont, + const xAOD::JetContainer* jets, + const xAOD::MissingETContainer* metCoreCont, + const xAOD::MissingETAssociationMap* map, + bool doJetJVF, + std::vector<const xAOD::IParticle*>& uniques) + { + ATH_MSG_VERBOSE("Rebuild jet term: " << metJetKey << " and soft term: " << softKey); + + MissingET* metJet = new MissingET(metJetKey,MissingETBase::Source::jet()); + metCont->push_back(metJet); + + const MissingET *coreSoftClus(0), *coreSoftTrk(0); + MissingET *metSoftClus(0), *metSoftTrk(0); + + const MissingET* coreSoft = (*metCoreCont)[softKey+"Core"]; + if(!coreSoft) { + ATH_MSG_WARNING("Invalid soft term key supplied: " << softKey); + return StatusCode::FAILURE; + } + if(MissingETBase::Source::isTrackTerm(coreSoft->source())) { + coreSoftTrk = coreSoft; + metSoftTrk = new MissingET(softKey,coreSoftTrk->source()); + metCont->push_back(metSoftTrk); + } else { + coreSoftClus = coreSoft; + metSoftClus = new MissingET(softKey,coreSoftClus->source()); + metCont->push_back(metSoftClus); + } + + return rebuildJetMET(metJet, jets, map, + metSoftClus, coreSoftClus, + metSoftTrk, coreSoftTrk, + doJetJVF,uniques); + } + + StatusCode METMaker::rebuildTrackMET(const std::string& metJetKey, + const std::string& softKey, + xAOD::MissingETContainer* metCont, + const xAOD::JetContainer* jets, + const xAOD::MissingETContainer* metCoreCont, + const xAOD::MissingETAssociationMap* map, + bool doJetJVF, + std::vector<const xAOD::IParticle*>& uniques) + { + ATH_MSG_VERBOSE("Rebuild jet term: " << metJetKey << " and soft term: " << softKey); + + MissingET* metJet = new MissingET(metJetKey,MissingETBase::Source::jet()); + metCont->push_back(metJet); + + const MissingET *coreSoftTrk(0); + MissingET *metSoftTrk(0); + + const MissingET* coreSoft = (*metCoreCont)[softKey+"Core"]; + if(!coreSoft) { + ATH_MSG_WARNING("Invalid soft term key supplied: " << softKey); + return StatusCode::FAILURE; + } + coreSoftTrk = coreSoft; + metSoftTrk = new MissingET(softKey,coreSoftTrk->source()); + metCont->push_back(metSoftTrk); + + return rebuildTrackMET(metJet, jets, map, + metSoftTrk, coreSoftTrk, + doJetJVF,uniques); + } + + inline StatusCode METMaker::rebuildJetMET(const std::string& metJetKey, + const std::string& softClusKey, + const std::string& softTrkKey, + xAOD::MissingETContainer* metCont, + const xAOD::JetContainer* jets, + const xAOD::MissingETContainer* metCoreCont, + const xAOD::MissingETAssociationMap* map, + bool doJetJVF) + { + std::vector<const xAOD::IParticle*> uniques; + return rebuildJetMET(metJetKey,softClusKey,softTrkKey,metCont, + jets,metCoreCont,map,doJetJVF); + } + + StatusCode METMaker::rebuildJetMET(const std::string& metJetKey, + const std::string& softClusKey, + const std::string& softTrkKey, + xAOD::MissingETContainer* metCont, + const xAOD::JetContainer* jets, + const xAOD::MissingETContainer* metCoreCont, + const xAOD::MissingETAssociationMap* map, + bool doJetJVF, + std::vector<const xAOD::IParticle*>& uniques) + { + ATH_MSG_VERBOSE("Create Jet MET " << metJetKey); + MissingET* metJet = new MissingET(metJetKey,MissingETBase::Source::jet()); + metCont->push_back(metJet); + ATH_MSG_VERBOSE("Create SoftClus MET " << softClusKey); + const MissingET* coreSoftClus = (*metCoreCont)[softClusKey+"Core"]; + ATH_MSG_VERBOSE("Create SoftTrk MET " << softTrkKey); + const MissingET* coreSoftTrk = (*metCoreCont)[softTrkKey+"Core"]; + if(!coreSoftClus) { + ATH_MSG_WARNING("Invalid cluster soft term key supplied: " << softClusKey); + return StatusCode::FAILURE; + } + if(!coreSoftTrk) { + ATH_MSG_WARNING("Invalid track soft term key supplied: " << softTrkKey); + return StatusCode::FAILURE; + } + MissingET* metSoftClus = new MissingET(softClusKey,coreSoftClus->source()); + metCont->push_back(metSoftClus); + MissingET* metSoftTrk = new MissingET(softTrkKey,coreSoftTrk->source()); + metCont->push_back(metSoftTrk); + + return rebuildJetMET(metJet, jets, map, + metSoftClus, coreSoftClus, + metSoftTrk, coreSoftTrk, + doJetJVF,uniques); + } StatusCode METMaker::rebuildJetMET(xAOD::MissingET* metJet, - const xAOD::JetContainer* jets, - const xAOD::MissingETAssociationMap* map, - std::vector<const xAOD::IParticle*>& uniques, - xAOD::MissingET* metSoftClus, - const xAOD::MissingET* coreSoftClus, - xAOD::MissingET* metSoftTrk, - const xAOD::MissingET* coreSoftTrk) { + const xAOD::JetContainer* jets, + const xAOD::MissingETAssociationMap* map, + xAOD::MissingET* metSoftClus, + const xAOD::MissingET* coreSoftClus, + xAOD::MissingET* metSoftTrk, + const xAOD::MissingET* coreSoftTrk, + bool doJetJVF, + std::vector<const xAOD::IParticle*>& uniques, + bool tracksForHardJets) { if(!metJet || !jets || !map) { ATH_MSG_WARNING( "Invalid pointer supplied for " - << "RefJet term (" << metJet << "), " - << "jet collection (" << jets << "), " - << " or map (" << map << ")." ); + << "RefJet term (" << metJet << "), " + << "jet collection (" << jets << "), " + << " or map (" << map << ")." ); return StatusCode::SUCCESS; } - ATH_MSG_VERBOSE("Rebuilding MET jet term " << metJet->name()); + ATH_MSG_VERBOSE("Building MET jet term " << metJet->name()); if(!metSoftClus && !metSoftTrk) { ATH_MSG_WARNING("Neither soft cluster nor soft track term has been supplied!"); return StatusCode::SUCCESS; } if(metSoftClus) { if(coreSoftClus) { - ATH_MSG_VERBOSE("Rebuilding MET soft cluster term " << metSoftClus->name()); - *metSoftClus += *coreSoftClus; + ATH_MSG_VERBOSE("Building MET soft cluster term " << metSoftClus->name()); + *metSoftClus += *coreSoftClus; } else { - ATH_MSG_WARNING("Soft cluster term provided without a core term!"); - return StatusCode::SUCCESS; + ATH_MSG_WARNING("Soft cluster term provided without a core term!"); + return StatusCode::SUCCESS; } } if(metSoftTrk) { if(coreSoftTrk) { - ATH_MSG_VERBOSE("Rebuilding MET soft track term " << metSoftTrk->name()); - *metSoftTrk += *coreSoftTrk; + ATH_MSG_VERBOSE("Building MET soft track term " << metSoftTrk->name()); + *metSoftTrk += *coreSoftTrk; } else { - ATH_MSG_WARNING("Soft track term provided without a core term!"); - return StatusCode::SUCCESS; + ATH_MSG_WARNING("Soft track term provided without a core term!"); + return StatusCode::SUCCESS; } } uniques.clear(); @@ -210,172 +405,183 @@ namespace met { for(const auto& jet : *jets) { const MissingETAssociation* assoc = 0; if(originalInputs) { - assoc = MissingETComposition::getAssociation(map,jet); + assoc = MissingETComposition::getAssociation(map,jet); } else { - const IParticle* orig = *m_objLinkAcc(*jet); - assoc = MissingETComposition::getAssociation(map,orig); + const IParticle* orig = *m_objLinkAcc(*jet); + assoc = MissingETComposition::getAssociation(map,static_cast<const xAOD::Jet*>(orig)); } if(assoc) { - ATH_MSG_VERBOSE( "Jet calib pt = " << jet->pt()); - bool selected = jet->pt()>m_jetMinPt; - if(m_doJetJvf) { - vector<float> jvf; - jet->getAttribute<vector<float> >(xAOD::JetAttribute::JVF,jvf); - selected = (jet->pt()>50e3 || fabs(jet->eta())>2.4 || fabs(jvf[0])>m_jetMinAbsJvf); - ATH_MSG_VERBOSE("Jet " << (selected ? "passes" : "fails") <<" JVF selection"); - } - bool hardJet(false); - MissingETBase::Types::constvec_t calvec = assoc->overlapCalVec(); - bool caloverlap = false; - caloverlap = calvec.ce()>0; - ATH_MSG_DEBUG("Jet " << jet->index() << " is " << ( caloverlap ? "" : "non-") << "overlapping"); - - xAOD::JetFourMom_t constjet = jet->jetP4("JetConstitScaleMomentum"); - double jpx = constjet.Px(); - double jpy = constjet.Py(); - double jpt = constjet.Pt(); - double opx = jpx - calvec.cpx(); - double opy = jpy - calvec.cpy(); - double opt = sqrt( opx*opx+opy*opy ); - double uniquefrac = 1. - (calvec.ce() / constjet.E()); - ATH_MSG_VERBOSE( "Jet constscale px, py, pt = " << jpx << ", " << jpy << ", " << jpt ); - ATH_MSG_VERBOSE( "Jet overlap E = " << calvec.ce() ); - ATH_MSG_VERBOSE( "Jet unique const E = " << constjet.E() - calvec.ce() ); - ATH_MSG_VERBOSE( "Jet OR px, py, pt = " << opx << ", " << opy << ", " << opt ); - if(selected) { - if(!caloverlap) { - // add jet full four-vector - hardJet = true; - *metJet += jet; - } else { - // check unique fraction - if(uniquefrac>m_jetMinEfrac && opt>m_jetMinWeightedPt) { - // add jet corrected for overlaps - hardJet = true; - if(m_jetCorrectPhi) { - metJet->add(opx/jpx*jet->px(),opy/jpy*jet->py(),opt/jpt*jet->pt()); - } else { - metJet->add(uniquefrac*jet->px(),uniquefrac*jet->py(),uniquefrac*jet->pt()); - } - } - } - } // hard jet selection - - if(hardJet){ - ATH_MSG_VERBOSE("Jet added at full scale"); - uniques.push_back(jet); - } else { - if(metSoftClus) { - // add fractional contribution - // FIXME: extend to allow pileup-suppressed versions - ATH_MSG_VERBOSE("Jet added at const scale"); - metSoftClus->add(opx,opy,opt); - } - if(metSoftTrk){ - // use jet tracks - // remove any tracks already used by other objects - MissingETBase::Types::constvec_t trkvec = assoc->overlapTrkVec(); - MissingETBase::Types::constvec_t jettrkvec = assoc->jetTrkVec(); - if(jettrkvec.ce()>1e-9) { - jpx = jettrkvec.cpx(); - jpy = jettrkvec.cpy(); - jpt = jettrkvec.sumpt(); - jettrkvec -= trkvec; - opx = jettrkvec.cpx(); - opy = jettrkvec.cpy(); - opt = jettrkvec.sumpt(); - ATH_MSG_VERBOSE( "Jet track px, py, sumpt = " << jpx << ", " << jpy << ", " << jpt ); - ATH_MSG_VERBOSE( "Jet OR px, py, sumpt = " << opx << ", " << opy << ", " << opt ); - } else { - opx = opy = opt = 0; - ATH_MSG_VERBOSE( "This jet has no associated tracks" ); - } - metSoftTrk->add(opx,opy,opt); - } // soft track - } // soft jet + ATH_MSG_VERBOSE( "Jet calib pt = " << jet->pt()); + bool selected = jet->pt()>m_jetMinPt; + if(doJetJVF) { + vector<float> jvf; + jet->getAttribute<vector<float> >(xAOD::JetAttribute::JVF,jvf); + selected = (jet->pt()>50e3 || fabs(jet->eta())>2.4 || fabs(jvf[0])>m_jetMinAbsJvf); + ATH_MSG_VERBOSE("Jet " << (selected ? "passes" : "fails") <<" JVF selection"); + } + bool hardJet(false); + MissingETBase::Types::constvec_t calvec = assoc->overlapCalVec(); + bool caloverlap = false; + caloverlap = calvec.ce()>0; + ATH_MSG_DEBUG("Jet " << jet->index() << " is " << ( caloverlap ? "" : "non-") << "overlapping"); + + xAOD::JetFourMom_t constjet = jet->jetP4("JetConstitScaleMomentum"); + double jpx = constjet.Px(); + double jpy = constjet.Py(); + double jpt = constjet.Pt(); + double opx = jpx - calvec.cpx(); + double opy = jpy - calvec.cpy(); + double opt = sqrt( opx*opx+opy*opy ); + double uniquefrac = 1. - (calvec.ce() / constjet.E()); + ATH_MSG_VERBOSE( "Jet constscale px, py, pt = " << jpx << ", " << jpy << ", " << jpt ); + ATH_MSG_VERBOSE( "Jet overlap E = " << calvec.ce() ); + ATH_MSG_VERBOSE( "Jet unique const E = " << constjet.E() - calvec.ce() ); + ATH_MSG_VERBOSE( "Jet OR px, py, pt = " << opx << ", " << opy << ", " << opt ); + if(selected) { + if(!caloverlap) { + // add jet full four-vector + hardJet = true; + if (!tracksForHardJets) *metJet += jet; + } else { + // check unique fraction + if(uniquefrac>m_jetMinEfrac && opt>m_jetMinWeightedPt) { + // add jet corrected for overlaps + hardJet = true; + if(m_jetCorrectPhi && !tracksForHardJets) { + metJet->add(opx/jpx*jet->px(),opy/jpy*jet->py(),opt/jpt*jet->pt()); + } else if (!tracksForHardJets) { + metJet->add(uniquefrac*jet->px(),uniquefrac*jet->py(),uniquefrac*jet->pt()); + } + } + } + } // hard jet selection + + if(hardJet){ + ATH_MSG_VERBOSE("Jet added at full scale"); + uniques.push_back(jet); + } else { + if(metSoftClus) { + // add fractional contribution + // FIXME: extend to allow pileup-suppressed versions + ATH_MSG_VERBOSE("Jet added at const scale"); + metSoftClus->add(opx,opy,opt); + } + } + if(metSoftTrk && (!hardJet || tracksForHardJets)){ + // use jet tracks + // remove any tracks already used by other objects + MissingETBase::Types::constvec_t trkvec = assoc->overlapTrkVec(); + MissingETBase::Types::constvec_t jettrkvec = assoc->jetTrkVec(); + if(jettrkvec.ce()>1e-9) { + jpx = jettrkvec.cpx(); + jpy = jettrkvec.cpy(); + jpt = jettrkvec.sumpt(); + jettrkvec -= trkvec; + opx = jettrkvec.cpx(); + opy = jettrkvec.cpy(); + opt = jettrkvec.sumpt(); + ATH_MSG_VERBOSE( "Jet track px, py, sumpt = " << jpx << ", " << jpy << ", " << jpt ); + ATH_MSG_VERBOSE( "Jet OR px, py, sumpt = " << opx << ", " << opy << ", " << opt ); + } else { + opx = opy = opt = 0; + ATH_MSG_VERBOSE( "This jet has no associated tracks" ); + } + if (hardJet) metJet->add(opx,opy,opt); + else metSoftTrk->add(opx,opy,opt); + } // soft track } // association exists else { - ATH_MSG_WARNING( "WARNING: Jet without association found!" ); + ATH_MSG_WARNING( "WARNING: Jet without association found!" ); } } // jet loop - + if(metSoftTrk) { // supplement track term with any tracks associated to isolated muons // these are recorded in the misc association const MissingETAssociation* assoc = map->getMiscAssociation(); if(assoc) { - MissingETBase::Types::constvec_t trkvec = assoc->overlapTrkVec(); - double opx = trkvec.cpx(); - double opy = trkvec.cpy(); - double osumpt = trkvec.sumpt();; - ATH_MSG_VERBOSE( "Misc track px, py, sumpt = " << opx << ", " << opy << ", " << osumpt ); - metSoftTrk->add(opx,opy,osumpt); + MissingETBase::Types::constvec_t trkvec = assoc->overlapTrkVec(); + double opx = trkvec.cpx(); + double opy = trkvec.cpy(); + double osumpt = trkvec.sumpt();; + ATH_MSG_VERBOSE( "Misc track px, py, sumpt = " << opx << ", " << opy << ", " << osumpt ); + metSoftTrk->add(opx,opy,osumpt); } } - + if(metSoftClus) { // supplement cluster term with any clusters associated to isolated e/gamma // these are recorded in the misc association const MissingETAssociation* assoc = map->getMiscAssociation(); if(assoc) { - MissingETBase::Types::constvec_t calvec = assoc->overlapCalVec(); - double opx = calvec.cpx(); - double opy = calvec.cpy(); - double osumpt = calvec.sumpt(); - ATH_MSG_VERBOSE( "Misc cluster px, py, sumpt = " << opx << ", " << opy << ", " << osumpt ); - metSoftClus->add(opx,opy,osumpt); + MissingETBase::Types::constvec_t calvec = assoc->overlapCalVec(); + double opx = calvec.cpx(); + double opy = calvec.cpy(); + double osumpt = calvec.sumpt(); + ATH_MSG_VERBOSE( "Misc cluster px, py, sumpt = " << opx << ", " << opy << ", " << osumpt ); + metSoftClus->add(opx,opy,osumpt); } } - + return StatusCode::SUCCESS; } - // **** Sum up MET terms **** + StatusCode METMaker::rebuildTrackMET(xAOD::MissingET* metJet, + const xAOD::JetContainer* jets, + const xAOD::MissingETAssociationMap* map, + xAOD::MissingET* metSoftTrk, + const xAOD::MissingET* coreSoftTrk, + bool doJetJVF, + std::vector<const xAOD::IParticle*>& uniques) { + return rebuildJetMET(metJet,jets,map,NULL,NULL,metSoftTrk,coreSoftTrk,doJetJVF,uniques,true); + } + // **** Sum up MET terms **** + StatusCode METMaker::buildMETSum(const std::string& totalName, - xAOD::MissingETContainer* metCont, - MissingETBase::Types::bitmask_t softTermsSource) + xAOD::MissingETContainer* metCont, + MissingETBase::Types::bitmask_t softTermsSource) { ATH_MSG_DEBUG("Build MET total: " << totalName); - + MissingET* metFinal = new MissingET(0.,0.,0.,totalName,MissingETBase::Source::total()); metCont->push_back(metFinal); - + for(const auto& met : *metCont) { - if(met->source()==MissingETBase::Source::total()) continue; + if(MissingETBase::Source::isTotalTerm(met->source())) continue; if(softTermsSource && MissingETBase::Source::isSoftTerm(met->source())) { - if(!MissingETBase::Source::hasPattern(met->source(),softTermsSource)) continue; + if(!MissingETBase::Source::hasPattern(met->source(),softTermsSource)) continue; } ATH_MSG_VERBOSE("Add MET term " << met->name() ); *metFinal += *met; } - + ATH_MSG_DEBUG( "Rebuilt MET Final --" - << " mpx: " << metFinal->mpx() - << " mpy: " << metFinal->mpy() - ); + << " mpx: " << metFinal->mpx() + << " mpy: " << metFinal->mpy() + ); return StatusCode::SUCCESS; } - + /////////////////////////////////////////////////////////////////// // Const methods: /////////////////////////////////////////////////////////////////// - + /////////////////////////////////////////////////////////////////// // Non-const methods: /////////////////////////////////////////////////////////////////// - + /////////////////////////////////////////////////////////////////// // Protected methods: /////////////////////////////////////////////////////////////////// - + /////////////////////////////////////////////////////////////////// // Const methods: /////////////////////////////////////////////////////////////////// - + /////////////////////////////////////////////////////////////////// // Non-const methods: /////////////////////////////////////////////////////////////////// - + } //> end namespace met diff --git a/Reconstruction/MET/METUtilities/Root/METRebuilder.cxx b/Reconstruction/MET/METUtilities/Root/METRebuilder.cxx index 6bbcbc362e400ca2187461f148cc171ab7c91894..8aa87c8eda64ef9da7e8e5a9b898868f026b7cd9 100644 --- a/Reconstruction/MET/METUtilities/Root/METRebuilder.cxx +++ b/Reconstruction/MET/METUtilities/Root/METRebuilder.cxx @@ -4,10 +4,10 @@ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration */ -// METRebuilder.cxx +// METRebuilder.cxx // Implementation file for class METRebuilder // Author: T.J.Khoo<khoo@cern.ch> -/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// // METUtilities includes #include "METUtilities/METRebuilder.h" @@ -30,52 +30,50 @@ // Tracking EDM #include "xAODTracking/TrackParticle.h" +#ifdef ROOTCORE +#include "InDetTrackSelectionTool/InDetTrackSelectionTool.h" +#endif // Track errors #include "EventPrimitives/EventPrimitivesHelpers.h" -namespace met { +#ifndef ROOTCORE +#include "GaudiKernel/IJobOptionsSvc.h" +#endif +namespace met { + + // Set up accessors to original object links in case of corrected copy containers + const static SG::AuxElement::Accessor<obj_link_t> objLinkAcc("originalObjectLink"); + using std::vector; - - using xAOD::MissingET; - using xAOD::MissingETContainer; - using xAOD::MissingETComponent; - using xAOD::MissingETComponentMap; - using xAOD::MissingETAuxContainer; - using xAOD::MissingETComposition; - // - using xAOD::IParticle; - using xAOD::IParticleContainer; - // - using xAOD::JetContainer; - using xAOD::JetConstituentVector; - // - using xAOD::TrackParticle; - - /////////////////////////////////////////////////////////////////// - // Public methods: - /////////////////////////////////////////////////////////////////// - + + using namespace xAOD; + + /////////////////////////////////////////////////////////////////// + // Public methods: + /////////////////////////////////////////////////////////////////// + // Constructors //////////////// - METRebuilder::METRebuilder(const std::string& name) : - AsgTool(name), - m_doEle(false), - m_doGamma(false), - m_doTau(false), - m_doMuon(false), - m_rebuildEle(false), - m_rebuildGamma(false), - m_rebuildTau(false), - m_rebuildMuon(false), - m_objLinkAcc("originalObjectLink"), - m_trkUsedDec("usedByMET"), - m_pureTrkSoft(true) + METRebuilder::METRebuilder(const std::string& name) : + AsgTool(name), + m_doEle(false), + m_doGamma(false), + m_doTau(false), + m_doMuon(false), + m_rebuildEle(false), + m_rebuildGamma(false), + m_rebuildTau(false), + m_rebuildMuon(false), + m_doTracks(true), + m_pureTrkSoft(true), + m_trkUsedDec("usedByMET"), + m_trkseltool("") { // // Property declaration - // + // declareProperty( "EleColl", m_eleColl = "ElectronCollection" ); declareProperty( "GammaColl", m_gammaColl = "PhotonCollection" ); declareProperty( "TauColl", m_tauColl = "TauRecContainer" ); @@ -87,125 +85,163 @@ namespace met { declareProperty( "TauTerm", m_tauTerm = "RefTau" ); declareProperty( "JetTerm", m_jetTerm = "RefJet" ); declareProperty( "MuonTerm", m_muonTerm = "Muons" ); - declareProperty( "SoftTerm", m_softTerm = "PVSoftTrk" ); + declareProperty( "SoftTerm", m_softTerm = "" ); // declareProperty( "InputMap", m_inputMap = "METMap_RefFinal" ); declareProperty( "OutputContainer", m_outMETCont = "MET_MyRefFinal" ); declareProperty( "OutputTotal", m_outMETTerm = "Final" ); declareProperty( "WarnIfDuplicate", m_warnOfDupes = true ); - - // migrate to new tool at some point + // + // for partial autoconfiguration + declareProperty( "SoftTermType", m_softTermType = "TrackSoftTerm" ); + // or "ClusterSoftTerm" or "Reference" + // declareProperty( "CalibJetPtCut", m_jetPtCut = 20e3 ); declareProperty( "DoJetJVFCut", m_jetDoJvf = true ); declareProperty( "CalibJetJvfCut", m_jetJvfCut = 0.25 ); - declareProperty( "SoftJetScale", m_softJetScale = "JetConstitScaleMomentum" ); - - // should be in a track tool -- also should use track isolation decorations + declareProperty( "SoftJetScale", m_softJetScale = "" ); + // declareProperty( "DoTrackPVSel", m_trk_doPVsel = true ); declareProperty( "VertexColl", m_vtxColl = "PrimaryVertices" ); declareProperty( "ClusterColl", m_clusColl = "CaloCalTopoCluster" ); - declareProperty( "TrackD0Max", m_trk_d0Max = 1.5 ); - declareProperty( "TrackZ0Max", m_trk_z0Max = 1.5 ); - + declareProperty( "TrkSelTool", m_trkseltool ); + } - + // Destructor /////////////// METRebuilder::~METRebuilder() {} - + // Athena algtool's Hooks //////////////////////////// StatusCode METRebuilder::initialize() { ATH_MSG_INFO ("Initializing " << name() << "..."); - + if( m_inputMap.size()==0 ) { ATH_MSG_FATAL("Input MissingETComponentMap name must be provided."); return StatusCode::FAILURE; } ATH_MSG_INFO ("Input MET Map: " << m_inputMap); - + if( m_outMETCont.size()==0 ) { ATH_MSG_FATAL("Output MissingETContainer name must be provided."); return StatusCode::FAILURE; } ATH_MSG_INFO ("Output MET Container: " << m_outMETCont); - + ATH_MSG_INFO ("Configured to rebuild following MET Terms:"); if( m_eleTerm!="" ) { m_doEle = true; if( m_eleColl!="" ) { - m_rebuildEle = true; - ATH_MSG_INFO(" Electrons: " << m_eleColl - << " > " << m_eleTerm ); + m_rebuildEle = true; + ATH_MSG_INFO(" Electrons: " << m_eleColl + << " > " << m_eleTerm ); } } if( m_gammaTerm!="" ) { m_doGamma = true; if( m_gammaColl!="" ) { - m_rebuildGamma = true; - ATH_MSG_INFO(" Photons: " << m_gammaColl - << " > " << m_gammaTerm ); + m_rebuildGamma = true; + ATH_MSG_INFO(" Photons: " << m_gammaColl + << " > " << m_gammaTerm ); } } if( m_tauTerm!="" ) { m_doTau = true; if( m_tauColl!="" ) { - m_rebuildTau = true; - ATH_MSG_INFO(" Taus: " << m_tauColl - << " > " << m_tauTerm ); + m_rebuildTau = true; + ATH_MSG_INFO(" Taus: " << m_tauColl + << " > " << m_tauTerm ); } } if( m_muonTerm!="" ) { m_doMuon = true; if( m_muonColl!="" ) { - m_rebuildMuon = true; - ATH_MSG_INFO(" Muons: " << m_muonColl - << " > " << m_muonTerm ); + m_rebuildMuon = true; + ATH_MSG_INFO(" Muons: " << m_muonColl + << " > " << m_muonTerm ); } } if( m_jetColl!="" && m_jetTerm!="") { ATH_MSG_INFO(" Jets: " << m_jetColl - << " > " << m_jetTerm ); + << " > " << m_jetTerm ); } else { ATH_MSG_FATAL("Error in configuration -- jet input and term keys must both be specified."); return StatusCode::FAILURE; } + if(m_softTerm.empty()) { + if(m_softTermType=="TrackSoftTerm") { + m_softTerm = "PVSoftTrk"; + m_jetDoJvf = true; + m_softJetScale = "JetTrackScale"; + } else if(m_softTermType=="ClusterSoftTerm") { + m_softTerm = "SoftClus"; + m_jetDoJvf = false; + m_softJetScale = "JetConstitScaleMomentum"; + } + if(m_softTermType=="Reference") { + m_softTerm = "SoftClus"; + m_jetDoJvf = false; + m_softJetScale = "JetConstitScaleMomentum"; + m_jetPtCut = 0.; + } + } ATH_MSG_INFO (" Soft: " << m_softTerm); - + m_trkUsedDec = SG::AuxElement::Decorator<char>("usedBy"+m_outMETCont); - - m_doTracks = (m_softTerm == "PVSoftTrk"); - m_pureTrkSoft = (m_softJetScale == "JetTrackScale"); - + + if(m_doTracks) { + if(m_trkseltool.empty()) { +#ifdef ROOTCORE + InDet::InDetTrackSelectionTool* trkSelTool = new InDet::InDetTrackSelectionTool("IDTrkSel_MET"); + ATH_CHECK( trkSelTool->setProperty("maxZ0SinTheta",1.5) ); + ATH_CHECK( trkSelTool->setProperty("maxD0overSigmaD0",3.) ); + trkSelTool->setCutLevel(InDet::CutLevel::Loose); + ATH_CHECK( trkSelTool->initialize() ); + m_trkseltool = ToolHandle<InDet::IInDetTrackSelectionTool>(trkSelTool); +#else + ServiceHandle<IJobOptionsSvc> josvc("JobOptionsSvc",name()); + std::string toolName = "IDTrkSel_METUtil"; + ATH_MSG_INFO("METRebuilder: Autoconfiguring " << toolName); + m_trkseltool.setTypeAndName("InDet::InDetTrackSelectionTool/"+toolName); + std::string fullToolName = "ToolSvc."+toolName; + ATH_CHECK( josvc->addPropertyToCatalogue(fullToolName,FloatProperty("maxZ0SinTheta",1.5)) ); + ATH_CHECK( josvc->addPropertyToCatalogue(fullToolName,FloatProperty("maxD0overSigmaD0",1.5)) ); + ATH_CHECK( josvc->addPropertyToCatalogue(fullToolName,StringProperty("CutLevel","Loose")) ); +#endif + } + ATH_CHECK( m_trkseltool.retrieve() ); + } + return StatusCode::SUCCESS; } - + StatusCode METRebuilder::finalize() { ATH_MSG_INFO ("Finalizing " << name() << "..."); - + return StatusCode::SUCCESS; } - + StatusCode METRebuilder::execute() { ATH_MSG_DEBUG ( name() << " in execute..."); - + const MissingETComponentMap* metMap = 0; if( evtStore()->retrieve(metMap, m_inputMap).isFailure() ) { ATH_MSG_WARNING("Unable to retrieve MissingETComponentMap: " << m_inputMap); return StatusCode::SUCCESS; } - + if( evtStore()->contains<MissingETContainer>(m_outMETCont) ) { if(m_warnOfDupes) - { ATH_MSG_WARNING("MET container " << m_outMETCont << " already in StoreGate"); } + { ATH_MSG_WARNING("MET container " << m_outMETCont << " already in StoreGate"); } return StatusCode::SUCCESS; } - + // Create a MissingETContainer with its aux store MissingETContainer* outCont = new MissingETContainer(); if( evtStore()->record(outCont, m_outMETCont).isFailure() ) { @@ -218,125 +254,116 @@ namespace met { return StatusCode::SUCCESS; } outCont->setStore(metAuxCont); - + if(m_doEle) { if(m_rebuildEle) { - const xAOD::ElectronContainer* elec = 0; - if( evtStore()->retrieve(elec, m_eleColl).isFailure() ) { - ATH_MSG_WARNING("Unable to retrieve ElectronContainer: " << m_eleColl); - return StatusCode::SUCCESS; - } - ATH_CHECK( rebuildMET(m_eleTerm, outCont, elec, metMap, m_doTracks) ); + const xAOD::ElectronContainer* elec = 0; + if( evtStore()->retrieve(elec, m_eleColl).isFailure() ) { + ATH_MSG_WARNING("Unable to retrieve ElectronContainer: " << m_eleColl); + return StatusCode::SUCCESS; + } + ATH_CHECK( rebuildMET(m_eleTerm, outCont, elec, metMap, m_doTracks) ); } else { - MissingET* metele = new MissingET(); - const MissingET* metele_ref = MissingETComposition::getMissingET(metMap,m_eleTerm); - outCont->push_back(metele); - metele->setName(metele_ref->name()); - metele->setSource(metele_ref->source()); - *metele += *metele_ref; + ATH_CHECK( copyMET(m_eleTerm,outCont,metMap) ); } } - + if(m_doGamma) { - if(m_rebuildGamma) { - const xAOD::PhotonContainer* gamma = 0; - if( evtStore()->retrieve(gamma, m_gammaColl).isFailure() ) { - ATH_MSG_WARNING("Unable to retrieve GammaContainer: " << m_gammaColl); - return StatusCode::SUCCESS; - } - ATH_CHECK( rebuildMET(m_gammaTerm, outCont, gamma, metMap, m_doTracks) ); - } else { - MissingET* metgamma = new MissingET(); - const MissingET* metgamma_ref = MissingETComposition::getMissingET(metMap,m_gammaTerm); - outCont->push_back(metgamma); - metgamma->setName(metgamma_ref->name()); - metgamma->setSource(metgamma_ref->source()); - *metgamma += *metgamma_ref; + if(m_rebuildGamma) { + const xAOD::PhotonContainer* gamma = 0; + if( evtStore()->retrieve(gamma, m_gammaColl).isFailure() ) { + ATH_MSG_WARNING("Unable to retrieve GammaContainer: " << m_gammaColl); + return StatusCode::FAILURE; + } + ATH_CHECK( rebuildMET(m_gammaTerm, outCont, gamma, metMap, m_doTracks) ); + } else { + ATH_CHECK( copyMET(m_gammaTerm,outCont,metMap) ); } } - - + if(m_doTau) { - if(m_rebuildTau) { - const xAOD::TauJetContainer* taujet = 0; - if( evtStore()->retrieve(taujet, m_tauColl).isFailure() ) { - ATH_MSG_WARNING("Unable to retrieve TauJetContainer: " << m_tauColl); - return StatusCode::SUCCESS; - } - ATH_CHECK( rebuildMET(m_tauTerm, outCont, taujet, metMap, m_doTracks) ); - } else { - MissingET* mettau = new MissingET(); - const MissingET* mettau_ref = MissingETComposition::getMissingET(metMap,m_tauTerm); - outCont->push_back(mettau); - mettau->setName(mettau_ref->name()); - mettau->setSource(mettau_ref->source()); - *mettau += *mettau_ref; + if(m_rebuildTau) { + const xAOD::TauJetContainer* taujet = 0; + if( evtStore()->retrieve(taujet, m_tauColl).isFailure() ) { + ATH_MSG_WARNING("Unable to retrieve TauJetContainer: " << m_tauColl); + return StatusCode::FAILURE; + } + ATH_CHECK( rebuildMET(m_tauTerm, outCont, taujet, metMap, m_doTracks) ); + } else { + ATH_CHECK( copyMET(m_tauTerm,outCont,metMap) ); } } - - + + if(m_doMuon) { - if(m_rebuildMuon) { - // May need implementation of Eloss correction - // Place in separate tool (?) - const xAOD::MuonContainer* muon = 0; - if( evtStore()->retrieve(muon, m_muonColl).isFailure() ) { - ATH_MSG_WARNING("Unable to retrieve MuonContainer: " << m_muonColl); - return StatusCode::SUCCESS; - } - ATH_CHECK( rebuildMET(m_muonTerm, outCont, muon, metMap, m_doTracks) ); - } else { - MissingET* metmuon = new MissingET(); - const MissingET* metmuon_ref = MissingETComposition::getMissingET(metMap,m_muonTerm); - outCont->push_back(metmuon); - metmuon->setName(metmuon_ref->name()); - metmuon->setSource(metmuon_ref->source()); - *metmuon += *metmuon_ref; + if(m_rebuildMuon) { + // May need implementation of Eloss correction + // Place in separate tool (?) + const xAOD::MuonContainer* muon = 0; + if( evtStore()->retrieve(muon, m_muonColl).isFailure() ) { + ATH_MSG_WARNING("Unable to retrieve MuonContainer: " << m_muonColl); + return StatusCode::FAILURE; + } + ATH_CHECK( rebuildMET(m_muonTerm, outCont, muon, metMap, m_doTracks) ); + } else { + ATH_CHECK( copyMET(m_muonTerm,outCont,metMap) ); } } - - + + // Implementation of the jet/soft term rebuilding // Place in separate tool (?) const xAOD::JetContainer* jet = 0; if( evtStore()->retrieve(jet, m_jetColl).isFailure() ) { ATH_MSG_WARNING("Unable to retrieve JetContainer: " << m_jetColl); - return StatusCode::SUCCESS; + return StatusCode::FAILURE; } ATH_CHECK( rebuildJetMET(m_jetTerm, m_softTerm, outCont, jet, metMap, m_doTracks) ); - + ATH_CHECK( buildMETSum(m_outMETTerm, outCont) ); - + return StatusCode::SUCCESS; } - + + StatusCode METRebuilder::copyMET(const std::string& metKey, + xAOD::MissingETContainer* metCont, + const xAOD::MissingETComponentMap* metMap) { + + MissingET* metterm = new MissingET(); + const MissingET* metterm_ref = MissingETComposition::getMissingET(metMap,metKey); + metCont->push_back(metterm); + *metterm = *metterm_ref; + + return StatusCode::SUCCESS; + } + // **** Rebuild generic MET term **** - + StatusCode METRebuilder::rebuildMET(const std::string& metKey, - xAOD::MissingETContainer* metCont, - const xAOD::IParticleContainer* collection, - const xAOD::MissingETComponentMap* metMap, - bool doTracks) { + xAOD::MissingETContainer* metCont, + const xAOD::IParticleContainer* collection, + const xAOD::MissingETComponentMap* metMap, + bool doTracks) { ATH_MSG_DEBUG("Rebuild MET term: " << metKey); const MissingETComponent* component = MissingETComposition::getComponent(metMap,metKey); if(!component) { ATH_MSG_WARNING("Could not find METComponent for " << metKey << " in MET Map!"); - return StatusCode::SUCCESS; + return StatusCode::FAILURE; } MissingET* met = new MissingET(0.,0.,0.,metKey,component->metObject()->source()); metCont->push_back(met); return rebuildMET(met,collection,component,doTracks); } - + StatusCode METRebuilder::rebuildMET(xAOD::MissingET* met, - const xAOD::IParticleContainer* collection, - const xAOD::MissingETComponent* component, - bool doTracks) { - + const xAOD::IParticleContainer* collection, + const xAOD::MissingETComponent* component, + bool doTracks) { + if(component->size()==0) return StatusCode::SUCCESS; - + ATH_MSG_VERBOSE("Rebuilding MET term " << component->metObject()->name()); - + const IParticleContainer* testCollection = dynamic_cast<const IParticleContainer*>(component->objects().front()->container()); bool originalInputs = (testCollection == collection); bool matchCollection = true; @@ -344,131 +371,133 @@ namespace met { // Consistency test: check that the collection supplied is the original one // used for MET reconstruction, or is derived from this collection if(!originalInputs) { - const IParticle* pObj = collection->front(); - if(!m_objLinkAcc.isAvailable(*pObj)) { - ATH_MSG_WARNING("Modified container provided without originalObjectLink -- cannot proceed."); - matchCollection = false; - } else { - const IParticleContainer* sourceCollection = dynamic_cast<const IParticleContainer*>((*m_objLinkAcc(*pObj))->container()); - matchCollection = (sourceCollection == testCollection); - } + const IParticle* pObj = collection->front(); + if(!objLinkAcc.isAvailable(*pObj)) { + ATH_MSG_WARNING("Modified container provided without originalObjectLink -- cannot proceed."); + matchCollection = false; + } else { + const IParticleContainer* sourceCollection = dynamic_cast<const IParticleContainer*>((*objLinkAcc(*pObj))->container()); + matchCollection = (sourceCollection == testCollection); + } } if(!matchCollection) { - ATH_MSG_WARNING("Problem with input object container -- skipping this term."); - return StatusCode::SUCCESS; + ATH_MSG_WARNING("Problem with input object container -- skipping this term."); + return StatusCode::SUCCESS; } } - + // Method flow: // 1. Loop over the objects in the collection // 2. Find them or their originals in the METComponent // 3. Add to the MET term with appropriate weights for( IParticleContainer::const_iterator iObj=collection->begin(); - iObj!=collection->end(); ++iObj ) { - + iObj!=collection->end(); ++iObj ) { + const IParticle* pObj = *iObj; // check if this is a copy - if so, get the original object pointer - if(!originalInputs) pObj = *m_objLinkAcc(*pObj); - + if(!originalInputs) pObj = *objLinkAcc(*pObj); + if(component->findIndex(pObj) != MissingETBase::Numerical::invalidIndex()) { - MissingETBase::Types::weight_t objWeight = component->weight(pObj); - ATH_MSG_VERBOSE( "Object with pt " << (*iObj)->pt() << " has weight " << objWeight.wet() ); - - if(doTracks) { - associateTracks(*iObj); - } - - met->add((*iObj)->pt()*cos((*iObj)->phi())*objWeight.wpx(), - (*iObj)->pt()*sin((*iObj)->phi())*objWeight.wpy(), - (*iObj)->pt()*objWeight.wet()); + MissingETBase::Types::weight_t objWeight = component->weight(pObj); + ATH_MSG_VERBOSE( "Object with pt " << (*iObj)->pt() << " has weight " << objWeight.wet() ); + + if(doTracks) { + associateTracks(*iObj); + } + + met->add((*iObj)->pt()*cos((*iObj)->phi())*objWeight.wpx(), + (*iObj)->pt()*sin((*iObj)->phi())*objWeight.wpy(), + (*iObj)->pt()*objWeight.wet()); } // used object in MET else { - ATH_MSG_VERBOSE( "Object with pt " << (*iObj)->pt() << " not found." ); + ATH_MSG_VERBOSE( "Object with pt " << (*iObj)->pt() << " not found." ); } } - + ATH_MSG_DEBUG( "Original " << component->metObject()->name() << " MET --" - << " mpx: " << component->metObject()->mpx() - << " mpy: " << component->metObject()->mpy() - ); + << " mpx: " << component->metObject()->mpx() + << " mpy: " << component->metObject()->mpy() + ); ATH_MSG_DEBUG( "Rebuilt " << component->metObject()->name() << " MET --" - << " mpx: " << met->mpx() - << " mpy: " << met->mpy() - ); - + << " mpx: " << met->mpx() + << " mpy: " << met->mpy() + ); + return StatusCode::SUCCESS; } - + // **** Rebuild jet & soft MET terms **** - + StatusCode METRebuilder::rebuildJetMET(const std::string& jetKey, - const std::string& softKey, - xAOD::MissingETContainer* metCont, - const xAOD::JetContainer* jets, - const xAOD::MissingETComponentMap* metMap, - bool doTracks) { + const std::string& softKey, + xAOD::MissingETContainer* metCont, + const xAOD::JetContainer* jets, + const xAOD::MissingETComponentMap* metMap, + bool doTracks, + bool doJvfCut, + bool pureTrkSoft, + const std::string& softJetScale) { const MissingETComponent* component = MissingETComposition::getComponent(metMap,jetKey); if(!component) { ATH_MSG_WARNING("Could not find METComponent for " << jetKey << " in MET Map!"); - return StatusCode::SUCCESS; + return StatusCode::FAILURE; } MissingET* metJet = new MissingET(0.,0.,0.,jetKey,component->metObject()->source()); metCont->push_back(metJet); - const MissingET* oldSoft = xAOD::MissingETComposition::getMissingET(metMap,softKey); - if(!component) { - ATH_MSG_WARNING("Could not find METComponent for " << softKey << " in MET Map!"); - return StatusCode::SUCCESS; - } - // need -1 because this assumes adding a particle, not a MET + ATH_CHECK( copyMET(softKey,metCont,metMap) ); // copy constructor needs correcting. - MissingET* metSoft = new MissingET(); - metCont->push_back(metSoft); - metSoft->setName(oldSoft->name()); - metSoft->setSource(oldSoft->source()); - *metSoft += *oldSoft; - return rebuildJetMET(metJet,metSoft,oldSoft,jets,component,doTracks); + MissingET* metSoft = (*metCont)[softKey]; + return rebuildJetMET(metJet,metSoft,jets,component,doTracks, + doJvfCut,pureTrkSoft,softJetScale); } - + StatusCode METRebuilder::rebuildJetMET(xAOD::MissingET* metJet, - xAOD::MissingET* metSoft, - const xAOD::MissingET* oldSoft, - const xAOD::JetContainer* jets, - const xAOD::MissingETComponent* component, - bool doTracks) { - + xAOD::MissingET* metSoft, + const xAOD::JetContainer* jets, + const xAOD::MissingETComponent* component, + bool doTracks, + bool doJvfCut, + bool pureTrkSoft, + const std::string& jetScale) { + if(component->size()==0) return StatusCode::SUCCESS; - - const xAOD::VertexContainer* vtx = 0; - if(m_trk_doPVsel && doTracks) { - if( evtStore()->retrieve( vtx, m_vtxColl).isFailure() ) { + + const VertexContainer* vtxCont = 0; + const Vertex* pv = 0; + if(doJvfCut || (m_trk_doPVsel && doTracks)) { + if( evtStore()->retrieve( vtxCont, m_vtxColl).isFailure() ) { ATH_MSG_WARNING("Unable to retrieve input primary vertex container"); - return StatusCode::SUCCESS; + return StatusCode::FAILURE; + } + for(const auto& vx : *vtxCont) { + if(vx->vertexType()==VxType::PriVtx) + {pv = vx; break;} } - if(vtx->size()>0) { - ATH_MSG_DEBUG("Main primary vertex has z = " << (*vtx)[0]->z()); + if(!pv) { + ATH_MSG_WARNING("Event has no primary vertex"); + return StatusCode::SUCCESS; } else{ - ATH_MSG_WARNING("Event has no primary vertices"); - return StatusCode::SUCCESS; + ATH_MSG_DEBUG("Main primary vertex has z = " << pv->z()); } } - - const IParticleContainer* testCollection = dynamic_cast<const IParticleContainer*>(component->objects().front()->container()); - const IParticleContainer* collcast = dynamic_cast<const IParticleContainer*>(jets); + + const IParticleContainer* testCollection = static_cast<const IParticleContainer*>(component->objects().front()->container()); + const IParticleContainer* collcast = static_cast<const IParticleContainer*>(jets); bool originalInputs = (testCollection == collcast); bool matchCollection = true; if(jets->size()>0) { // Consistency test: check that the collection supplied is the original one // used for MET reconstruction, or is derived from this collection if(!originalInputs) { - const IParticle* pJet = jets->front(); - if(!m_objLinkAcc.isAvailable(*pJet)) { - ATH_MSG_WARNING("Modified container provided without originalObjectLink -- cannot proceed."); - matchCollection = false; - } else { - const IParticleContainer* sourceCollection = dynamic_cast<const IParticleContainer*>((*m_objLinkAcc(*pJet))->container()); - matchCollection = (sourceCollection == testCollection); - } + const IParticle* pJet = jets->front(); + if(!objLinkAcc.isAvailable(*pJet)) { + ATH_MSG_WARNING("Modified container provided without originalObjectLink -- cannot proceed."); + matchCollection = false; + } else { + const IParticleContainer* sourceCollection = static_cast<const IParticleContainer*>((*objLinkAcc(*pJet))->container()); + matchCollection = (sourceCollection == testCollection); + } } } if(!matchCollection) { @@ -479,208 +508,197 @@ namespace met { // 2. Find them or their originals in the METComponent // 3. Add to the MET term with appropriate weights for( JetContainer::const_iterator iJet=jets->begin(); - iJet!=jets->end(); ++iJet ) { - + iJet!=jets->end(); ++iJet ) { + const xAOD::IParticle* pJet = *iJet; - if(!originalInputs) pJet = *m_objLinkAcc(*pJet); - + if(!originalInputs) pJet = *objLinkAcc(*pJet); + if(component->findIndex(pJet) != MissingETBase::Numerical::invalidIndex()) { - - MissingETBase::Types::weight_t jetWeight = component->weight(pJet); - bool passJVF = true; - if(m_jetDoJvf) { - vector<float> jvf; - (*iJet)->getAttribute<vector<float> >(xAOD::JetAttribute::JVF,jvf); - passJVF = (*iJet)->pt()>50e3 || fabs((*iJet)->eta())>2.4 || fabs(jvf[0])>m_jetJvfCut; - ATH_MSG_VERBOSE("Jet with pt " << (*iJet)->pt() << " has jvf = " << jvf[0]); - } - - if((*iJet)->pt()>m_jetPtCut && passJVF) { - - ATH_MSG_VERBOSE("Retain jet with pt " << (*iJet)->pt() << " at full scale."); - - metJet->add((*iJet)->px()*jetWeight.wpx(), - (*iJet)->py()*jetWeight.wpy(), - (*iJet)->pt()*jetWeight.wet()); - } // minimum pt cut for jet calibration - else { - double trkjetpx(0), trkjetpy(0), trkjetpt(0); - if(doTracks) { - ATH_MSG_VERBOSE("Add tracks from jet with pt " << (*iJet)->pt()); - vector<const TrackParticle*> jettracks = (*iJet)->getAssociatedObjects<TrackParticle>(xAOD::JetAttribute::GhostTrack); - ATH_MSG_VERBOSE("Got jet tracks"); - for(vector<const TrackParticle*>::const_iterator iTrk = jettracks.begin(); - iTrk!=jettracks.end(); ++iTrk) { - if(!*iTrk) continue; - // duplicate ST track selection -- should be in a tool - if(fabs((*iTrk)->pt())>500/*MeV*/ && fabs((*iTrk)->eta())<2.5) { - uint8_t nPixHits(0), nSctHits(0); - (*iTrk)->summaryValue(nPixHits,xAOD::numberOfPixelHits); - (*iTrk)->summaryValue(nSctHits,xAOD::numberOfSCTHits); - if(nPixHits>=1 && nSctHits>=6) { - bool badTrack = false; - if( (fabs((*iTrk)->eta())<1.5 && (*iTrk)->pt()>200e3) || - (fabs((*iTrk)->eta())>=1.5 && (*iTrk)->pt()>120e3) ) { - // Get relative error on qoverp - float Rerr = Amg::error((*iTrk)->definingParametersCovMatrix(),4)/fabs((*iTrk)->qOverP()); - // Simplified cut -- remove tracks that are more energetic than the jet - if(Rerr>0.4 || (*iTrk)->pt()>2*(*iJet)->pt()) badTrack = true; - } // additional cuts against high pt mismeasured tracks - bool uniqueTrack = true; - uniqueTrack = !m_trkUsedDec(**iTrk); - if(!badTrack && uniqueTrack && - (!m_trk_doPVsel||isPVTrack(*iTrk,(*vtx)[0]))) { - ATH_MSG_VERBOSE(" + track with pt " << (*iTrk)->pt()); - trkjetpx += (*iTrk)->pt()*cos((*iTrk)->phi()); - trkjetpy += (*iTrk)->pt()*sin((*iTrk)->phi()); - trkjetpt += (*iTrk)->pt(); - } else { ATH_MSG_VERBOSE(" - track failed badtrack/uniqueness/PV"); } - } // PIX & SCT hits - else { ATH_MSG_VERBOSE(" - track failed hit requirement"); } - } // pt, eta - else { ATH_MSG_VERBOSE(" - track failed pt/eta"); } - } // track loop - } // track-based soft term - if(m_pureTrkSoft){ - metSoft->add(trkjetpx,trkjetpy,trkjetpt); - } else { - // just add the weighted constituent-scale jet - xAOD::JetFourMom_t jetP = (*iJet)->jetP4(m_softJetScale); - ATH_MSG_VERBOSE("Soft jet pt = " << jetP.Pt() << ", track pt = " << trkjetpt); - if(trkjetpt>jetP.Pt() || !passJVF) { // use tracks if higher scale than calo jet or fail JVF cut in central region - ATH_MSG_VERBOSE("Add jet with pt " << (*iJet)->pt() - << " at track scale (pt = " << trkjetpt << ")."); - metSoft->add(trkjetpx,trkjetpy,trkjetpt); - } else { - ATH_MSG_VERBOSE("Add jet with pt " << (*iJet)->pt() - << " at constituent scale (pt = " << jetP.Pt() << ")."); - metSoft->add(jetP.Px()*jetWeight.wpx(), - jetP.Py()*jetWeight.wpy(), - jetP.Pt()*jetWeight.wet()); - } // otherwise, use calo jet at chosen scale - } // cluster-based soft term - } // jets below threshold should be added to the soft terms + + MissingETBase::Types::weight_t jetWeight = component->weight(pJet); + bool passJVF = true; + if(doJvfCut) { + vector<float> jvf; + (*iJet)->getAttribute<vector<float> >(xAOD::JetAttribute::JVF,jvf); + passJVF = (*iJet)->pt()>50e3 || fabs((*iJet)->eta())>2.4 || fabs(jvf[pv->index()])>m_jetJvfCut; + ATH_MSG_VERBOSE("Jet with pt " << (*iJet)->pt() << " has jvf = " << jvf[pv->index()]); + } + + if((*iJet)->pt()>m_jetPtCut && passJVF) { + + ATH_MSG_VERBOSE("Retain jet with pt " << (*iJet)->pt() << " at full scale."); + + metJet->add((*iJet)->px()*jetWeight.wpx(), + (*iJet)->py()*jetWeight.wpy(), + (*iJet)->pt()*jetWeight.wet()); + } // minimum pt cut for jet calibration + else { + double trkjetpx(0), trkjetpy(0), trkjetpt(0); + if(doTracks) { + ATH_MSG_VERBOSE("Add tracks from jet with pt " << (*iJet)->pt()); + vector<const TrackParticle*> jettracks = (*iJet)->getAssociatedObjects<TrackParticle>(xAOD::JetAttribute::GhostTrack); + ATH_MSG_VERBOSE("Got jet tracks"); + for(vector<const TrackParticle*>::const_iterator iTrk = jettracks.begin(); + iTrk!=jettracks.end(); ++iTrk) { + if(!*iTrk) continue; + bool badTrack = false; + if( (fabs((*iTrk)->eta())<1.5 && (*iTrk)->pt()>200e3) || + (fabs((*iTrk)->eta())>=1.5 && (*iTrk)->pt()>120e3) ) { + // Get relative error on qoverp + float Rerr = Amg::error((*iTrk)->definingParametersCovMatrix(),4)/fabs((*iTrk)->qOverP()); + // Simplified cut -- remove tracks that are more energetic than the jet + if(Rerr>0.4 || (*iTrk)->pt()>2*(*iJet)->pt()) badTrack = true; + } // additional cuts against high pt mismeasured tracks + bool uniqueTrack = true; + uniqueTrack = !m_trkUsedDec(**iTrk); + if(!badTrack && uniqueTrack && (acceptTrack(*iTrk,pv))) { + ATH_MSG_VERBOSE(" + track with pt " << (*iTrk)->pt()); + trkjetpx += (*iTrk)->pt()*cos((*iTrk)->phi()); + trkjetpy += (*iTrk)->pt()*sin((*iTrk)->phi()); + trkjetpt += (*iTrk)->pt(); + } else { ATH_MSG_VERBOSE(" - track failed badtrack/uniqueness/PV"); } + } // track loop + } // track-based soft term + if(pureTrkSoft){ + metSoft->add(trkjetpx,trkjetpy,trkjetpt); + } else { + // just add the weighted constituent-scale jet + xAOD::JetFourMom_t jetP = (*iJet)->jetP4(jetScale); + ATH_MSG_VERBOSE("Soft jet pt = " << jetP.Pt() << ", track pt = " << trkjetpt); + if(trkjetpt>jetP.Pt() || !passJVF) { // use tracks if higher scale than calo jet or fail JVF cut in central region + ATH_MSG_VERBOSE("Add jet with pt " << (*iJet)->pt() + << " at track scale (pt = " << trkjetpt << ")."); + metSoft->add(trkjetpx,trkjetpy,trkjetpt); + } else { + ATH_MSG_VERBOSE("Add jet with pt " << (*iJet)->pt() + << " at constituent scale (pt = " << jetP.Pt() << ")."); + metSoft->add(jetP.Px()*jetWeight.wpx(), + jetP.Py()*jetWeight.wpy(), + jetP.Pt()*jetWeight.wet()); + } // otherwise, use calo jet at chosen scale + } // cluster-based soft term + } // jets below threshold should be added to the soft terms } // used jet in MET } - + ATH_MSG_DEBUG( "Original jet MET --" - << " mpx: " << component->metObject()->mpx() - << " mpy: " << component->metObject()->mpy() - ); + << " mpx: " << component->metObject()->mpx() + << " mpy: " << component->metObject()->mpy() + ); ATH_MSG_DEBUG( "Rebuilt jet MET --" - << " mpx: " << metJet->mpx() - << " mpy: " << metJet->mpy() - ); - - ATH_MSG_DEBUG("Original MET Soft --" - << " mpx: " << oldSoft->mpx() - << " mpy: " << oldSoft->mpy()); - + << " mpx: " << metJet->mpx() + << " mpy: " << metJet->mpy() + ); + ATH_MSG_DEBUG( "Rebuilt MET soft --" - << " mpx: " << metSoft->mpx() - << " mpy: " << metSoft->mpy() - ); - + << " mpx: " << metSoft->mpx() + << " mpy: " << metSoft->mpy() + ); + return StatusCode::SUCCESS; } - + // **** Sum up MET terms **** - + StatusCode METRebuilder::buildMETSum(const std::string& totalName, - xAOD::MissingETContainer* metCont) + xAOD::MissingETContainer* metCont) { ATH_MSG_DEBUG("Build MET total: " << totalName); - + MissingET* metFinal = new MissingET(0.,0.,0.,"Final",MissingETBase::Source::total()); metCont->push_back(metFinal); - + for(MissingETContainer::const_iterator iMET=metCont->begin(); - iMET!=metCont->end(); ++iMET) { + iMET!=metCont->end(); ++iMET) { if(*iMET==metFinal) continue; *metFinal += **iMET; } - + ATH_MSG_DEBUG( "Rebuilt MET Final --" - << " mpx: " << metFinal->mpx() - << " mpy: " << metFinal->mpy() - ); + << " mpx: " << metFinal->mpx() + << " mpy: " << metFinal->mpy() + ); return StatusCode::SUCCESS; } - - /////////////////////////////////////////////////////////////////// - // Const methods: + + /////////////////////////////////////////////////////////////////// + // Const methods: + /////////////////////////////////////////////////////////////////// + + /////////////////////////////////////////////////////////////////// + // Non-const methods: + /////////////////////////////////////////////////////////////////// + /////////////////////////////////////////////////////////////////// - - /////////////////////////////////////////////////////////////////// - // Non-const methods: - /////////////////////////////////////////////////////////////////// - - /////////////////////////////////////////////////////////////////// // Protected methods: /////////////////////////////////////////////////////////////////// - + // Implement for now, but should move to common tools when possible - bool METRebuilder::isPVTrack(const xAOD::TrackParticle* trk, - const xAOD::Vertex* pv) const + bool METRebuilder::acceptTrack(const xAOD::TrackParticle* trk, + const xAOD::Vertex* pv) const { - - if(trk->d0()>m_trk_d0Max) return false; - if(fabs(trk->z0()+trk->vz() - pv->z()) > m_trk_z0Max) return false; - - return true; + + // if(trk->d0()>m_trk_d0Max) return false; + // if(fabs(trk->z0()+trk->vz() - pv->z()) > m_trk_z0Max) return false; + if(m_trk_doPVsel) {return m_trkseltool->accept( *trk, pv );} + else {return m_trkseltool->accept( trk );} } - + void METRebuilder::associateTracks(const xAOD::IParticle* obj) { if(obj->type()==xAOD::Type::Electron) { - const xAOD::Electron* el = dynamic_cast<const xAOD::Electron*>(obj); + const xAOD::Electron* el = static_cast<const xAOD::Electron*>(obj); for(size_t iTrk=0; iTrk<el->nTrackParticles(); ++iTrk) { - const TrackParticle* eltrk = xAOD::EgammaHelpers::getOriginalTrackParticleFromGSF(el->trackParticle(iTrk)); - if(eltrk) m_trkUsedDec(*eltrk) = true; + const TrackParticle* eltrk = xAOD::EgammaHelpers::getOriginalTrackParticleFromGSF(el->trackParticle(iTrk)); + if(eltrk) m_trkUsedDec(*eltrk) = true; } } if(obj->type()==xAOD::Type::Photon) { - const xAOD::Photon* ph = dynamic_cast<const xAOD::Photon*>(obj); + const xAOD::Photon* ph = static_cast<const xAOD::Photon*>(obj); for(size_t iVtx=0; iVtx<ph->nVertices(); ++iVtx) { - const xAOD::Vertex* phvx = ph->vertex(iVtx); - - if(phvx) { - for(size_t iTrk=0; iTrk<phvx->nTrackParticles(); ++iTrk) { - const TrackParticle* phtrk = xAOD::EgammaHelpers::getOriginalTrackParticleFromGSF(phvx->trackParticle(iTrk)); - if(phtrk) m_trkUsedDec(*phtrk) = true; - } - } - } - } - if(obj->type()==xAOD::Type::Tau) { - const xAOD::TauJet* tau = dynamic_cast<const xAOD::TauJet*>(obj); - // now find associated tracks - for(size_t iTrk=0; iTrk<tau->nTracks(); ++iTrk) { - m_trkUsedDec(*tau->track(iTrk)) = true; - } - for(size_t iTrk=0; iTrk<tau->nOtherTracks(); ++iTrk) { - const xAOD::TrackParticle* trk = tau->otherTrack(iTrk); - double dR = (*tau->jetLink())->p4().DeltaR(trk->p4()); - if(dR<0.2) { - m_trkUsedDec(*trk) = true; - } + const xAOD::Vertex* phvx = ph->vertex(iVtx); + + if(phvx) { + for(size_t iTrk=0; iTrk<phvx->nTrackParticles(); ++iTrk) { + const TrackParticle* phtrk = xAOD::EgammaHelpers::getOriginalTrackParticleFromGSF(phvx->trackParticle(iTrk)); + if(phtrk) m_trkUsedDec(*phtrk) = true; + } + } } } + if(obj->type()==xAOD::Type::Tau) { + const xAOD::TauJet* tau = static_cast<const xAOD::TauJet*>(obj); + // now find associated tracks + const std::vector< ElementLink< xAOD::TrackParticleContainer > >& trackLinks = tau->trackLinks(); + for(const auto& trklink : trackLinks) { + if(trklink.isValid()) {m_trkUsedDec(**trklink) = true;} + else { ATH_MSG_WARNING("Invalid tau track link"); } + } + const std::vector< ElementLink< xAOD::TrackParticleContainer > >& otherTrackLinks = tau->otherTrackLinks(); + for(const auto& trklink : otherTrackLinks) { + if(trklink.isValid()) { + double dR = (*tau->jetLink())->p4().DeltaR((*trklink)->p4()); + if(dR<0.2) { + m_trkUsedDec(**trklink) = true; + } + } else { ATH_MSG_WARNING("Invalid tau track link"); } + } + } if(obj->type()==xAOD::Type::Muon) { - const xAOD::Muon* mu = dynamic_cast<const xAOD::Muon*>(obj); + const xAOD::Muon* mu = static_cast<const xAOD::Muon*>(obj); if(mu->inDetTrackParticleLink().isValid()) { - m_trkUsedDec(**mu->inDetTrackParticleLink()) = true; + m_trkUsedDec(**mu->inDetTrackParticleLink()) = true; } } } - + /////////////////////////////////////////////////////////////////// // Const methods: /////////////////////////////////////////////////////////////////// - + /////////////////////////////////////////////////////////////////// // Non-const methods: /////////////////////////////////////////////////////////////////// - + } //> end namespace met diff --git a/Reconstruction/MET/METUtilities/Root/METSystematicsTool.cxx b/Reconstruction/MET/METUtilities/Root/METSystematicsTool.cxx new file mode 100644 index 0000000000000000000000000000000000000000..aeaeff2dc517e121d581f2db40c823fb672805ea --- /dev/null +++ b/Reconstruction/MET/METUtilities/Root/METSystematicsTool.cxx @@ -0,0 +1,709 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#include "METUtilities/METSystematicsTool.h" +#include "TFile.h" +#include "TEnv.h" +#include "TSystem.h" +#include "TH3.h" +#include "TH2.h" + +#include <iostream> +// #include <boost/filesystem.hpp> +#include "PathResolver/PathResolver.h" + +// xAOD includes +#include "xAODCore/ShallowCopy.h" +#include "xAODTracking/VertexContainer.h" +#include "xAODMissingET/MissingETAuxContainer.h" +#include "xAODMissingET/MissingETComposition.h" +#include "xAODJet/JetContainer.h" + + +namespace met { + + using namespace xAOD; + + METSystematicsTool::METSystematicsTool(const std::string& name) : asg::AsgTool::AsgTool(name), + m_appliedSystEnum(NONE), + shiftpara_pthard_njet_mu(nullptr), + resopara_pthard_njet_mu(nullptr), + resoperp_pthard_njet_mu(nullptr), + jet_systRpt_pt_eta(nullptr), + h_calosyst_scale(nullptr), + h_calosyst_reso (nullptr), + m_rand(0) + { + ATH_MSG_DEBUG (__PRETTY_FUNCTION__ ); + + declareProperty( "JetColl", m_jetColl = "AntiKt4LCTopoJets" ); + declareProperty( "ConfigPrefix", m_configPrefix = "METUtilities/data12_8TeV/final2/" ); + declareProperty( "ConfigSoftTrkFile", m_configSoftTrkFile = "TrackSoftTerms.config" ); + declareProperty( "ConfigJetTrkFile", m_configJetTrkFile = "" ); + declareProperty( "ConfigSoftCaloFile",m_configSoftCaloFile= "METRefFinal.config" ); + declareProperty( "TruthContainer", m_truthCont = "MET_Truth" ); + declareProperty( "TruthObj", m_truthObj = "NonInt" ); + declareProperty( "VertexContainer", m_vertexCont = "PrimaryVertices" ); + declareProperty( "EventInfo", m_eventInfo = "EventInfo" ); + + applySystematicVariation(CP::SystematicSet()).ignore(); + } + + StatusCode METSystematicsTool::addMETAffectingSystematics(){ + ATH_MSG_VERBOSE (__PRETTY_FUNCTION__ ); + + if(!m_configSoftTrkFile.empty()){ + if( !(addAffectingSystematic( softTrkAffSyst::MET_SoftTrk_ScaleUp , true /*recommended */ ) && + addAffectingSystematic( softTrkAffSyst::MET_SoftTrk_ScaleDown, true /*recommended */ ) && + addAffectingSystematic( softTrkAffSyst::MET_SoftTrk_ResoPara , true /*recommended */ ) && + addAffectingSystematic( softTrkAffSyst::MET_SoftTrk_ResoPerp , true /*recommended */ ) && + addAffectingSystematic( softTrkAffSyst::MET_SoftTrk_ResoCorr , false /*not recommended */) ) ) { + ATH_MSG_ERROR("failed to properly add softTrk affecting systematics " ); + return StatusCode::FAILURE; + } + } + if(!m_configSoftCaloFile.empty()){ + if( !(addAffectingSystematic( softCaloAffSyst::MET_SoftCalo_ScaleUp , true /*recommended */ ) && + addAffectingSystematic( softCaloAffSyst::MET_SoftCalo_ScaleDown, true /*recommended */ ) && + addAffectingSystematic( softCaloAffSyst::MET_SoftCalo_Reso , true /*recommended */ ) ) ) { + ATH_MSG_ERROR("failed to properly add softCalo affecting systematics " ); + return StatusCode::FAILURE; + } + } + if(!m_configJetTrkFile.empty()){ + if( !(addAffectingSystematic( jetTrkAffSyst::MET_JetTrk_ScaleUp , true /*recommended */ ) && + addAffectingSystematic( jetTrkAffSyst::MET_JetTrk_ScaleDown, true /*recommended */ ) ) ){ + ATH_MSG_ERROR("failed to properly add jetTrk affecting systematics " ); + return StatusCode::FAILURE; + } + } + + ATH_MSG_INFO( "AffectingSystematics are:" ); + for(const auto& syst : m_affectingSystematics) { + ATH_MSG_INFO( syst.name() ); + } + + ATH_MSG_DEBUG("These systematics are set based on your config files: " ); + ATH_MSG_DEBUG("ConfigSoftTrkFile: " << m_configSoftTrkFile ); + ATH_MSG_DEBUG("ConfigSoftCaloFile: "<< m_configSoftCaloFile); + ATH_MSG_DEBUG("ConfigSoftJetFile: " << m_configJetTrkFile ); + + return StatusCode::SUCCESS; + } + + StatusCode METSystematicsTool::initialize() + { + ATH_MSG_VERBOSE (__PRETTY_FUNCTION__ ); + + const char lastchar = m_configPrefix.back(); + if(std::strncmp(&lastchar,"/",1)!=0) { + m_configPrefix.append("/"); + } + + if(m_configSoftTrkFile.empty() && + m_configSoftCaloFile.empty() && + m_configJetTrkFile.empty() ) ATH_MSG_WARNING("you have initialized the METSystematicsTool with no configuration file. The tool will do nothing. Please set the configuration file properties."); + + if(!m_configSoftTrkFile.empty()) ATH_CHECK( softTrkSystInitialize() ) ; + if(!m_configSoftCaloFile.empty()) ATH_CHECK( softCaloSystInitialize()); + if(!m_configJetTrkFile.empty()) ATH_CHECK( jetTrkSystInitialize() ); + + return addMETAffectingSystematics(); + } + + StatusCode METSystematicsTool::finalize() + { + ATH_MSG_DEBUG (__PRETTY_FUNCTION__); + + delete shiftpara_pthard_njet_mu; + delete resopara_pthard_njet_mu; + delete resoperp_pthard_njet_mu; + delete jet_systRpt_pt_eta; + delete h_calosyst_scale; + delete h_calosyst_reso; + + return StatusCode::SUCCESS; + } + + StatusCode METSystematicsTool::softCaloSystInitialize() + { + ATH_MSG_VERBOSE (__PRETTY_FUNCTION__); + + std::string histfile = ""; + std::string gsystpath = ""; + std::string histpath = ""; + std::string blank = ""; + + ATH_CHECK( extractHistoPath(histfile,gsystpath,histpath,blank ,SOFTCALO) );//properly sets the paths + + TFile infile(histpath.c_str()); + + ATH_MSG_INFO( "METSystematics: Read calo uncertainties" ); + h_calosyst_scale = dynamic_cast<TH1D*>( infile.Get((gsystpath+"/globsyst_scale").c_str())); + h_calosyst_reso = dynamic_cast<TH1D*> (infile.Get((gsystpath+"/globsyst_reso").c_str())); + + if( !(h_calosyst_scale && + h_calosyst_reso + ) + ) + { + ATH_MSG_ERROR("Could not get all calo histos from the config file:" << histfile); + return StatusCode::FAILURE; + } + + h_calosyst_scale->SetDirectory(0); + h_calosyst_reso ->SetDirectory(0); + + ATH_MSG_VERBOSE( __PRETTY_FUNCTION__ << " DONE!!" ); + return StatusCode::SUCCESS; + } + + + StatusCode METSystematicsTool::jetTrkSystInitialize() + { + ATH_MSG_VERBOSE (__PRETTY_FUNCTION__); + + std::string histfile = ""; + std::string gsystpath = ""; + std::string configdir = ""; + std::string blank = ""; + + ATH_CHECK( extractHistoPath(histfile,gsystpath, configdir,blank , JETTRK) );//properly sets the paths + + TFile infile((configdir).c_str()); + + jet_systRpt_pt_eta = dynamic_cast<TH2D*>( infile.Get("jet_systRpt_pt_eta") ) ; + + if( !jet_systRpt_pt_eta) + { + ATH_MSG_ERROR("Could not get jet track histo from the config file:" << histfile); + return StatusCode::FAILURE; + } + + jet_systRpt_pt_eta->SetDirectory(0); + + ATH_MSG_VERBOSE( __PRETTY_FUNCTION__ << " DONE!!"); + return StatusCode::SUCCESS; + } + + + StatusCode METSystematicsTool::softTrkSystInitialize() + { + ATH_MSG_VERBOSE (__PRETTY_FUNCTION__ ); + + std::string histfile = ""; + std::string psystpath = ""; + std::string histpath = ""; + std::string suffix = ""; + + ATH_CHECK( extractHistoPath(histfile,psystpath,histpath, suffix, SOFTTRK) );//properly sets the paths + + TFile infile(histpath.c_str()); + + resoperp_pthard_njet_mu = dynamic_cast<TH3D*>( infile.Get((psystpath+"/resoperp_"+suffix).c_str()) ); + resopara_pthard_njet_mu = dynamic_cast<TH3D*>( infile.Get((psystpath+"/resopara_"+suffix).c_str()) ); + shiftpara_pthard_njet_mu = dynamic_cast<TH3D*>( infile.Get((psystpath+"/shiftpara_"+suffix).c_str())); + + if( !(resoperp_pthard_njet_mu && + resopara_pthard_njet_mu && + shiftpara_pthard_njet_mu + ) + ) + { + ATH_MSG_ERROR("Could not get all track histos from the config file:" << histfile ); + return StatusCode::FAILURE; + } + + resoperp_pthard_njet_mu ->SetDirectory(0); + resopara_pthard_njet_mu ->SetDirectory(0); + shiftpara_pthard_njet_mu->SetDirectory(0); + + ATH_MSG_VERBOSE( __PRETTY_FUNCTION__ <<" DONE!!"); + return StatusCode::SUCCESS; + } + + CP::SystematicCode METSystematicsTool::sysApplySystematicVariation(const CP::SystematicSet& systSet){//this should already be filtered for MET systematics + ATH_MSG_VERBOSE (__PRETTY_FUNCTION__ ); + // Only a single systematic can be applied at a time: + // If at some point we can deal with multiple systematics, we will check here that the combination we are given will work + + if( systSet.size()==0 ) { + ATH_MSG_DEBUG("No affecting systematics received."); + return CP::SystematicCode::Ok; + } else if( systSet.size() > 1 ) { + ATH_MSG_WARNING("Tool does not support multiple systematics, returning unsupported" ); + return CP::SystematicCode::Unsupported; + } + CP::SystematicVariation systVar = *systSet.begin(); + if ( systVar == CP::SystematicVariation("") ) m_appliedSystEnum = NONE ; + else if( systVar == softTrkAffSyst::MET_SoftTrk_ScaleUp) m_appliedSystEnum = MET_SOFTTRK_SCALEUP ; + else if( systVar == softTrkAffSyst::MET_SoftTrk_ScaleDown) m_appliedSystEnum = MET_SOFTTRK_SCALEDOWN ; + else if( systVar == softTrkAffSyst::MET_SoftTrk_ResoPara) m_appliedSystEnum = MET_SOFTTRK_RESOPARA ; + else if( systVar == softTrkAffSyst::MET_SoftTrk_ResoPerp) m_appliedSystEnum = MET_SOFTTRK_RESOPERP ; + else if( systVar == softTrkAffSyst::MET_SoftTrk_ResoCorr) m_appliedSystEnum = MET_SOFTTRK_RESOCORR ; + else if( systVar == softCaloAffSyst::MET_SoftCalo_ScaleUp) m_appliedSystEnum = MET_SOFTCALO_SCALEUP ; + else if( systVar == softCaloAffSyst::MET_SoftCalo_ScaleDown) m_appliedSystEnum = MET_SOFTCALO_SCALEDOWN; + else if( systVar == softCaloAffSyst::MET_SoftCalo_Reso) m_appliedSystEnum = MET_SOFTCALO_RESO ; + else if( systVar == jetTrkAffSyst::MET_JetTrk_ScaleUp) m_appliedSystEnum = MET_JETTRK_SCALEUP ; + else if( systVar == jetTrkAffSyst::MET_JetTrk_ScaleDown) m_appliedSystEnum = MET_JETTRK_SCALEDOWN ; + else{ + ATH_MSG_WARNING("unsupported systematic applied " ); + return CP::SystematicCode::Unsupported; + } + + ATH_MSG_DEBUG("applied systematic is " << m_appliedSystEnum ); + + return CP::SystematicCode::Ok; + } + + CP::CorrectionCode METSystematicsTool::applyCorrection(xAOD::MissingET& inputMet, xAOD::MissingETAssociationMap * map) const{//if asking for jet track systematics, the user needs to give a met association map as well + ATH_MSG_VERBOSE (__PRETTY_FUNCTION__ ); + + if( MissingETBase::Source::isSoftTerm(inputMet.source())){ + if(getDefaultEventInfo() == nullptr){ + ATH_MSG_WARNING("event info is empty, returning without applying correction"); + return CP::CorrectionCode::Error; + } + xAOD::MissingETContainer const * METcont = static_cast<xAOD::MissingETContainer const*>(inputMet.container()); + if(METcont == nullptr){ + ATH_MSG_WARNING("MissingET object not owned by a container. Returning without applying correction" ); + return CP::CorrectionCode::Error; + } + return internalSoftTermApplyCorrection(inputMet, METcont, *getDefaultEventInfo()); + } + + if( MissingETBase::Source::isTrackTerm(inputMet.source()) ){//todo other source check? //todo + if( map == nullptr) { + ATH_MSG_WARNING("To calculate jet track systematics, you must give applyCorrection a MissingETAssociationMap, as applyCorrection(inputMet, map). Returning without applying correction "); + return CP::CorrectionCode::Error; + } + + xAOD::JetContainer const * jetCont = nullptr; + if(evtStore()->retrieve(jetCont, m_jetColl).isFailure()){ + ATH_MSG_WARNING( "retrieval of jet container " << m_jetColl << " failed" ); + return CP::CorrectionCode::Error; + } + + return getCorrectedJetTrackMET(inputMet, map, jetCont); + } + + ATH_MSG_ERROR("METSystematicsTool received a MissingET object it can't correct"); + return CP::CorrectionCode::Error;//todo should this be an error? + } + + + CP::CorrectionCode METSystematicsTool::correctedCopy(const xAOD::MissingET& met, xAOD::MissingET*& outputmet, xAOD::MissingETAssociationMap * map) const + { ATH_MSG_VERBOSE (__PRETTY_FUNCTION__ ); + xAOD::MissingET * copy = new xAOD::MissingET(met); + + if( MissingETBase::Source::isSoftTerm(met.source())){ + if(getDefaultEventInfo() == nullptr){ + ATH_MSG_WARNING("event info is empty, returning empty with unknown source"); + outputmet = nullptr; delete copy; + return CP::CorrectionCode::Error; + } + + xAOD::MissingETContainer const * METcont = dynamic_cast<xAOD::MissingETContainer const*>(met.container()); + if(METcont == nullptr){ + ATH_MSG_WARNING("MissingET object not owned by a container. Unable to apply correction, returning output MET object as null" ); + outputmet = nullptr; delete copy; + return CP::CorrectionCode::Error; + } + + if(internalSoftTermApplyCorrection(*copy, METcont, *getDefaultEventInfo()) != CP::CorrectionCode::Ok ){ + outputmet = nullptr; delete copy; + return CP::CorrectionCode::Error; + } + }//soft term source + //todo check this should be if not else if + if(MissingETBase::Source::isTrackTerm(met.source()) ){//todo additional source check? + if( map == nullptr) { + ATH_MSG_WARNING("To calculate jet track systematics, you must give correctedCopy a MissingETAssociationMap, as correctedCopy(inputMet, map). "); + outputmet = nullptr; delete copy; + return CP::CorrectionCode::Error; + } + + + xAOD::JetContainer const * jetCont = nullptr; + + if(evtStore()->retrieve(jetCont, m_jetColl).isFailure()){ + ATH_MSG_WARNING( "retrieval of jet container " << m_jetColl << " failed" ); + outputmet = nullptr; delete copy; + return CP::CorrectionCode::Error; + } + + if(getCorrectedJetTrackMET(*copy, map, jetCont) != CP::CorrectionCode::Ok){ + outputmet = nullptr; delete copy; + return CP::CorrectionCode::Error; + } + }//track term source + + outputmet = copy; + return CP::CorrectionCode::Ok; + } + + CP::CorrectionCode METSystematicsTool::internalSoftTermApplyCorrection(xAOD::MissingET & softMet, + xAOD::MissingETContainer const * METcont, + xAOD::EventInfo const & eInfo) const{ //this is equivalent of "getSoftTerms" + + ATH_MSG_VERBOSE( __PRETTY_FUNCTION__ ); + + + if( ! MissingETBase::Source::isSoftTerm(softMet.source()) ){ + ATH_MSG_WARNING("not soft met, cannot apply soft term correction to this MET"); //todo should this be a warning? + return CP::CorrectionCode::Error; + } + + bool doSyst = false; + if( MissingETBase::Source::isTrackTerm(softMet.source()) ) { + doSyst = m_appliedSystEnum>=MET_SOFTTRK_SCALEUP && m_appliedSystEnum<=MET_SOFTTRK_RESOCORR; + } else { + doSyst = m_appliedSystEnum>=MET_SOFTCALO_SCALEUP && m_appliedSystEnum<=MET_SOFTCALO_RESO; + } + + if(doSyst) { + + xAOD::JetContainer const * jetCont = nullptr; + + if(METcont == nullptr){ + ATH_MSG_WARNING("failed to retrieve MET container from passed object"); + return CP::CorrectionCode::Error; + } + if(evtStore()->retrieve(jetCont, m_jetColl).isFailure()){ + ATH_MSG_WARNING( "retrieval of jet container " << m_jetColl << " failed" ); + return CP::CorrectionCode::Error; + } + + xAOD::MissingET const ptHard = calcPtHard(METcont);//GeV? + + int phbin = shiftpara_pthard_njet_mu->GetXaxis()->FindBin(ptHard.met()); + if(phbin>shiftpara_pthard_njet_mu->GetNbinsX()) phbin = shiftpara_pthard_njet_mu->GetNbinsX(); + int const jetbin = shiftpara_pthard_njet_mu->GetYaxis()->FindBin(jetCont->size()); + int const mubin = shiftpara_pthard_njet_mu->GetZaxis()->FindBin(eInfo.actualInteractionsPerCrossing() ); + double const ptHardShift = shiftpara_pthard_njet_mu->GetBinContent(phbin,jetbin,mubin); + + double const randGaus = m_rand.Gaus(0.,1.); + + ATH_MSG_DEBUG("About to apply systematic " << appliedSystematicsString() ); + + //now we need to know what soft term systematics we are doing + //m_appliedSystEnum was cached by the applySystematicVariation method + switch( m_appliedSystEnum ){ + case MET_SOFTTRK_SCALEUP : { + softMet = softTrkSyst_scale(softMet, ptHard, ptHardShift); + break; + } + case MET_SOFTTRK_SCALEDOWN: { + softMet = softTrkSyst_scale(softMet, ptHard,-1*ptHardShift); + break; + } + case MET_SOFTTRK_RESOPARA : { + double const smearpara = resopara_pthard_njet_mu->GetBinContent(phbin,jetbin,mubin)*randGaus; + softMet = softTrkSyst_reso(softMet, ptHard, ptHardShift, smearpara, 0.); + break; + } + case MET_SOFTTRK_RESOPERP : { + double const smearperp = resoperp_pthard_njet_mu->GetBinContent(phbin,jetbin,mubin)*randGaus; + softMet = softTrkSyst_reso(softMet, ptHard, ptHardShift, 0., smearperp ); + break; + } + case MET_SOFTTRK_RESOCORR : { + double const smearpara = resopara_pthard_njet_mu->GetBinContent(phbin,jetbin,mubin)*randGaus; + double const smearperp = resoperp_pthard_njet_mu->GetBinContent(phbin,jetbin,mubin)*randGaus; + softMet = softTrkSyst_reso(softMet, ptHard, ptHardShift, smearpara , smearperp); + break; + } + case MET_SOFTCALO_SCALEUP : { + double const caloscale = 1. + h_calosyst_scale->GetBinContent(1); + softMet = caloSyst_scale(softMet,caloscale); + break; + } + case MET_SOFTCALO_SCALEDOWN : { + double const caloscale = 1. - h_calosyst_scale->GetBinContent(1); + softMet = caloSyst_scale(softMet,caloscale); + break; + } + case MET_SOFTCALO_RESO : { + softMet = caloSyst_reso(softMet) ; + break; + } + default:{ + ATH_MSG_DEBUG("No systematic applied, returning nominal MET term"); + } + } + } else { + ATH_MSG_DEBUG("Ignore irrelevant systematic."); + } + + return CP::CorrectionCode::Ok; + } + + CP::CorrectionCode METSystematicsTool::calcJetTrackMETWithSyst(xAOD::MissingET& jettrkmet, + const xAOD::MissingETAssociationMap* map, + const xAOD::Jet* jet) const + + { + ATH_MSG_VERBOSE(__PRETTY_FUNCTION__); + + if( jet_systRpt_pt_eta == nullptr ) { + ATH_MSG_ERROR("jet track systematics histogram not initialized properly.") ; + return CP::CorrectionCode::Error; + } + xAOD::MissingETAssociation const * const assoc = MissingETComposition::getAssociation(map,jet); + MissingETBase::Types::constvec_t trkvec = assoc->jetTrkVec(); + + int phbin = jet_systRpt_pt_eta->GetXaxis()->FindBin(jet->pt()/1e3); + if(phbin>jet_systRpt_pt_eta->GetNbinsX()) phbin = jet_systRpt_pt_eta->GetNbinsX(); + + int etabin = jet_systRpt_pt_eta->GetYaxis()->FindBin(fabs( jet->eta() )); + if(etabin>jet_systRpt_pt_eta->GetNbinsY()) etabin = jet_systRpt_pt_eta->GetNbinsY(); + + double uncert = 0.; + switch( m_appliedSystEnum ){ + + case MET_JETTRK_SCALEUP : { + uncert = jet_systRpt_pt_eta->GetBinContent(phbin,etabin); + break; + } + case MET_JETTRK_SCALEDOWN : { + uncert = -1.*jet_systRpt_pt_eta->GetBinContent(phbin,etabin); + break; + } + default:{ + ATH_MSG_DEBUG("No jet track systematic applied"); + break; + } + } + + jettrkmet.setMpx ( jettrkmet.mpx() - trkvec.cpx() * uncert); + jettrkmet.setMpy ( jettrkmet.mpy() - trkvec.cpy() * uncert); + + return CP::CorrectionCode::Ok; + } + + CP::CorrectionCode METSystematicsTool::getCorrectedJetTrackMET(xAOD::MissingET& jettrkmet, + const xAOD::MissingETAssociationMap* map, + const xAOD::JetContainer * vecJets + ) const + { + ATH_MSG_VERBOSE( __PRETTY_FUNCTION__ ); + if(!map) { + ATH_MSG_ERROR("MissingETAssociationMap null, error."); + return CP::CorrectionCode::Error; + } + //todo cleanup statuses + + for(xAOD::JetContainer::const_iterator iJet=vecJets->begin(); iJet!=vecJets->end(); ++iJet ){ + if(calcJetTrackMETWithSyst(jettrkmet, map , *iJet) != CP::CorrectionCode::Ok) return CP::CorrectionCode::Error; + } + return CP::CorrectionCode::Ok; + } + + xAOD::MissingET METSystematicsTool::caloSyst_scale(xAOD::MissingET const &softTerms, double const scale) const{ + ATH_MSG_VERBOSE( __PRETTY_FUNCTION__ ); + return xAOD::MissingET(softTerms.mpx()*scale, softTerms.mpy()*scale, softTerms.sumet(), + softTerms.name(),softTerms.source()); + } + + xAOD::MissingET METSystematicsTool::caloSyst_reso(xAOD::MissingET const &softTerms) const { + ATH_MSG_VERBOSE( __PRETTY_FUNCTION__ ); + ATH_MSG_VERBOSE("caloSyst_reso: input MET: " << softTerms.met()); + + double const metSigma = .7 * sqrt(softTerms.sumet()); + double const resUnc = h_calosyst_reso->GetBinContent(1); + double const smearedSigma = sqrt( (metSigma* (1. + resUnc))*(metSigma* (1. + resUnc)) - + metSigma * metSigma ); + + ATH_MSG_VERBOSE("caloSyst_reso: metSigma: " << metSigma << ", resUnc: " << resUnc << ", smearedSigma = " << smearedSigma); + + double const rand = gRandom->Gaus(0.,1.); + double const shift = softTerms.met()<1e-9 ? 0. : rand*smearedSigma / softTerms.met(); + + ATH_MSG_VERBOSE("caloSyst_reso: shift = " << shift); + + return xAOD::MissingET(softTerms.mpx()*(1.+shift),softTerms.mpy()*(1.+shift),softTerms.sumet(), + softTerms.name(),softTerms.source()); + } + + xAOD::MissingET METSystematicsTool::softTrkSyst_scale(xAOD::MissingET const &softTerms, xAOD::MissingET const &ptHard, double const shift) const + { ATH_MSG_VERBOSE(__PRETTY_FUNCTION__); + return softTrkSyst_reso(softTerms, ptHard, shift, 0. , 0.); + } + + xAOD::MissingET METSystematicsTool::softTrkSyst_reso(xAOD::MissingET const &softTerms, + xAOD::MissingET const &ptHard, + double const shift, + double const smearpara, + double const smearperp) const{ + ATH_MSG_VERBOSE( __PRETTY_FUNCTION__ ); + + xAOD::MissingET projection = projectST(softTerms,ptHard); + projection.setMpx(projection.mpx() + shift + smearpara ); + projection.setMpy(projection.mpy() + + smearperp ); + + return projectST(projection, ptHard); + } + + xAOD::MissingET METSystematicsTool::projectST(xAOD::MissingET const &softTerms, xAOD::MissingET const &ptHard) const + { + ATH_MSG_VERBOSE( __PRETTY_FUNCTION__ ); + + double const ux = ptHard.mpx()/ptHard.met(); + double const uy = ptHard.mpy()/ptHard.met(); + double const projL = ux*softTerms.mpx() + uy*softTerms.mpy(); + double const projT = uy*softTerms.mpx() - ux*softTerms.mpy(); + xAOD::MissingET proj(projL,projT,softTerms.sumet(), + softTerms.name(),softTerms.source()); + + return proj; + } + + xAOD::MissingET METSystematicsTool::calcPtHard(xAOD::MissingETContainer const * const cont) const + { + ATH_MSG_VERBOSE(__PRETTY_FUNCTION__ ); + + //get truth container + xAOD::MissingETContainer const * truthCont = nullptr; + if(evtStore()->retrieve(truthCont, m_truthCont).isFailure()){ + ATH_MSG_ERROR( m_truthCont << " container empty or does not exist, calcPtHard returning zero."); + return xAOD::MissingET(0.,0.,0.); + } + + xAOD::MissingET const * truthptr = (*truthCont)[MissingETBase::Source::truthNonInt()]; + if(truthptr == nullptr){ + ATH_MSG_ERROR( m_truthObj << " is not in " << m_truthCont << ". calcPtHard returing zero." ); + return xAOD::MissingET(0.,0.,0.); + } + + xAOD::MissingET ptHard = (*truthptr); + + //loop over all non soft or total terms + for(xAOD::MissingETContainer::const_iterator iMET=cont->begin();iMET!=cont->end(); ++iMET ) { + if( !( MissingETBase::Source::isSoftTerm( (*iMET)->source() ) || + MissingETBase::Source::isTotalTerm( (*iMET)->source() ) ) ) { + ptHard -= **iMET; + } + } + + ptHard *= 1./1000. ; + + ATH_MSG_DEBUG("PtHard value is " << std::to_string(ptHard.met()) ); + + return ptHard; + } + + StatusCode METSystematicsTool::extractHistoPath(std::string & histfile , + std::string & systpath , + std::string & histpath, + std::string & suffix , + SystType const & type + ) + { + ATH_MSG_VERBOSE (__PRETTY_FUNCTION__); + TEnv reader; + + std::string configpath = ""; + std::string configfile = ""; + + switch(type){ + case SOFTCALO : configfile = m_configSoftCaloFile; + configpath = PathResolver::find_file(m_configPrefix+m_configSoftCaloFile, "CALIBPATH", PathResolver::RecursiveSearch) ; + break; + case SOFTTRK : configfile = m_configSoftTrkFile; + configpath = PathResolver::find_file(m_configPrefix+m_configSoftTrkFile, "CALIBPATH", PathResolver::RecursiveSearch) ; + break; + case JETTRK : configfile = m_configJetTrkFile; + configpath = PathResolver::find_file(m_configPrefix+m_configJetTrkFile, "CALIBPATH", PathResolver::RecursiveSearch) ; + break; + default : configpath = ""; + } + + ATH_MSG_INFO( "Searching for configFile: " << configfile); + ATH_MSG_INFO( "PWD: " << gSystem->Getenv("PWD") ) ; + ATH_MSG_INFO( "CALIBPATH: " << gSystem->Getenv("CALIBPATH") ); + + if(configpath.empty() || configfile.empty() ){ + ATH_MSG_ERROR( "Path Resolver couldn't find config file"); + return StatusCode::FAILURE; + } + + // boost::filesystem::path configpathBf(configpath); + // configdir = configpathBf.parent_path().string(); + + if( reader.ReadFile(configpath.c_str(),EEnvLevel(0)) < 0) { + ATH_MSG_ERROR( "Couldn't read the track config file!" ); + return StatusCode::FAILURE; + } + + ATH_MSG_INFO( "Configuring from file " << configpath ); + + + switch(type){ + case SOFTCALO : + histfile = reader.GetValue( "Conf.InputFile" , ""); + systpath = reader.GetValue( "GlobalSyst.sourcedir" , "" ); + break; + case SOFTTRK : + histfile = reader.GetValue( "Conf.InputFile" , ""); + systpath = reader.GetValue( "PtHardSyst.sourcedir" , "" ); + suffix = reader.GetValue( "PtHardSyst.suffix" , "" ); + break; + case JETTRK : + histfile = reader.GetValue( "JetTrkSyst.InputFile" , ""); + systpath = "/"; + break; + default : break; + } + + ATH_MSG_INFO( "Will read histograms from " << histfile ); + + histpath = PathResolver::find_file(histfile, "CALIBPATH", PathResolver::RecursiveSearch) ; + ATH_MSG_INFO("Extracted histogram path " << histpath); + + if(histfile.empty() || systpath.empty() || histpath.empty() ){ + ATH_MSG_ERROR("Failed to correctly set histfile path, or path to histograms inside of the histfile" ); + return StatusCode::FAILURE; + } + return StatusCode::SUCCESS; + } + + + //stolen from JetUncertainties + xAOD::EventInfo const * METSystematicsTool::getDefaultEventInfo() const + { + ATH_MSG_VERBOSE (__PRETTY_FUNCTION__ ); + + const xAOD::EventInfo* eInfoConst = nullptr; + if (evtStore()->retrieve(eInfoConst,m_eventInfo).isFailure()){ + ATH_MSG_ERROR("Failed to retrieve default EventInfo object"); + return nullptr; + } + + return eInfoConst; + } + + int METSystematicsTool::getNPV() const{ + ATH_MSG_VERBOSE (__PRETTY_FUNCTION__ ); + const xAOD::VertexContainer* vertices = nullptr; + + if (evtStore()->retrieve(vertices,m_vertexCont).isFailure()){ + ATH_MSG_ERROR("Failed to retrieve default NPV value from PrimaryVertices"); + return 0; + } + + int NPV = 0; + xAOD::VertexContainer::const_iterator itr; + for (itr = vertices->begin(); itr != vertices->end(); ++itr) + if ( (*itr)->vertexType() != 0) NPV++; + + return NPV; + } + + void METSystematicsTool::setRandomSeed(int seed){ + ATH_MSG_VERBOSE(__PRETTY_FUNCTION__); + m_rand.SetSeed(seed); + } +} + +// LocalWords: SOFTTRK diff --git a/Reconstruction/MET/METUtilities/cmt/Makefile.RootCore b/Reconstruction/MET/METUtilities/cmt/Makefile.RootCore index 3481b8accfe0d2f03c12187bce0d7cea0c53ef45..c35a9ef030d0d9930a92c0aad9d88cc3a5ee3bba 100644 --- a/Reconstruction/MET/METUtilities/cmt/Makefile.RootCore +++ b/Reconstruction/MET/METUtilities/cmt/Makefile.RootCore @@ -12,14 +12,15 @@ PACKAGE_OBJFLAGS = PACKAGE_LDFLAGS = PACKAGE_BINFLAGS = PACKAGE_LIBFLAGS = -PACKAGE_DEP = AsgTools xAODBase xAODTracking xAODMissingET xAODEgamma xAODTruth xAODJet xAODMuon xAODTau xAODPFlow METInterface +PACKAGE_DEP = AsgTools xAODBase xAODTracking xAODMissingET xAODEgamma xAODTruth xAODJet xAODMuon xAODTau xAODPFlow METInterface xAODEventInfo PATInterfaces PathResolver InDetTrackSelectionTool PACKAGE_TRYDEP = PACKAGE_CLEAN = PACKAGE_NOGRID = -PACKAGE_PEDANTIC = 0 +PACKAGE_PEDANTIC = 1 PACKAGE_NOOPT = 0 PACKAGE_NOCC = 0 +# Temporarily disable reflex dicts to work with IInDetTrackSelectionTool "enum class" PACKAGE_REFLEX = 0 include $(ROOTCOREDIR)/Makefile-common diff --git a/Reconstruction/MET/METUtilities/cmt/requirements b/Reconstruction/MET/METUtilities/cmt/requirements index 502d9f2e6decf3ed7755073f89b6f2e6d49ff785..c9c7abfc3d57ba65e6a606475db6ae4d2ef0c03a 100644 --- a/Reconstruction/MET/METUtilities/cmt/requirements +++ b/Reconstruction/MET/METUtilities/cmt/requirements @@ -10,18 +10,25 @@ use GaudiInterface GaudiInterface-* External ## framework dependencies use AsgTools AsgTools-* Control/AthToolSupport +use AtlasROOT AtlasROOT-* External +use PATInterfaces PATInterfaces-* PhysicsAnalysis/AnalysisCommon ## put here your package dependencies... use METInterface METInterface-* Reconstruction/MET use xAODJet xAODJet-* Event/xAOD - use xAODTracking xAODTracking-* Event/xAOD +use xAODMissingET xAODMissingET-* Event/xAOD + +use InDetTrackSelectionTool InDetTrackSelectionTool-* InnerDetector/InDetRecTools private -use xAODMissingET xAODMissingET-* Event/xAOD use xAODEgamma xAODEgamma-* Event/xAOD use xAODMuon xAODMuon-* Event/xAOD use xAODTau xAODTau-* Event/xAOD +use xAODCore xAODCore-* Event/xAOD +use xAODBase xAODBase-* Event/xAOD + +use PathResolver PathResolver-* Tools # needed for track momentum errors use EventPrimitives EventPrimitives-* Event @@ -34,11 +41,23 @@ use AthenaBaseComps AthenaBaseComps-* Control public branches src src/components doc python share METUtilities Root +#library METUtilitiesLib ../Root/METHelper.cxx +#apply_pattern named_installed_library library=METUtilitiesLib + +#todo probably remove when we get files in PathResolver space +apply_pattern declare_calib files="../data/*.config ../data/*.root" + private -library METUtilities *.cxx ../Root/*.cxx -s=components *.cxx +library METUtilities ../src/*.cxx ../Root/*.cxx -s=components *.cxx apply_pattern component_library apply_pattern declare_joboptions files="*.py" apply_pattern declare_python_modules files="*.py" -end_private +#end_private + +#macro_append METUtilities_dependencies " METUtilitiesLib" + +#private +use AtlasReflex AtlasReflex-* External -no-auto-imports +apply_pattern lcgdict dict=METUtilities selectionfile=selection.xml headerfiles="../METUtilities/METUtilitiesDict.h" diff --git a/Reconstruction/MET/METUtilities/data/data12_8TeV/final2/METRefFinal.config b/Reconstruction/MET/METUtilities/data/data12_8TeV/final2/METRefFinal.config new file mode 100644 index 0000000000000000000000000000000000000000..7f4ff4e89518f9aa243796434c16a5639f8bd4e7 --- /dev/null +++ b/Reconstruction/MET/METUtilities/data/data12_8TeV/final2/METRefFinal.config @@ -0,0 +1,9 @@ +############################################################################## +# Configuration file for MET_RefFinal +############################################################################## +Conf.DoJetBins: false +Conf.InputFile: METUtilities/data12_8TeV/final2/SoftTermsSyst_METRefFinal_STVF.root +Conf.DoGlobalSyst: true + +GlobalSyst.sourcedir: /syst2012/RefFinal/globalsyst +PtHardSyst.sourcedir: /syst2012/RefFinal/pthardsyst diff --git a/Reconstruction/MET/METUtilities/data/data12_8TeV/final2/METRefFinalSTVF.config b/Reconstruction/MET/METUtilities/data/data12_8TeV/final2/METRefFinalSTVF.config new file mode 100644 index 0000000000000000000000000000000000000000..b37465984713ced39cd58eb0a3af2be351d6e273 --- /dev/null +++ b/Reconstruction/MET/METUtilities/data/data12_8TeV/final2/METRefFinalSTVF.config @@ -0,0 +1,9 @@ +############################################################################## +# Configuration file for MET_RefFinal_STVF +############################################################################## +Conf.DoJetBins: false +Conf.InputFile: METUtilities/data12_8TeV/final2/SoftTermsSyst_METRefFinal_STVF.root +Conf.DoGlobalSyst: true + +GlobalSyst.sourcedir: /syst2012/STVF/globalsyst +PtHardSyst.sourcedir: /syst2012/STVF/pthardsyst diff --git a/Reconstruction/MET/METUtilities/data/data12_8TeV/final2/METTrack_2012.config b/Reconstruction/MET/METUtilities/data/data12_8TeV/final2/METTrack_2012.config new file mode 100644 index 0000000000000000000000000000000000000000..deb81d9fe040a107d328928616aa407e64972486 --- /dev/null +++ b/Reconstruction/MET/METUtilities/data/data12_8TeV/final2/METTrack_2012.config @@ -0,0 +1,11 @@ +############################################################################## +# Configuration file for MET_Track +############################################################################## +Conf.InputFile: SoftTermsSyst_METTrack.root +Conf.DoGlobalSyst: false + +PtHardSyst.sourcedir: / +PtHardSyst.suffix: 2012 + +Conf.DoJetTrkSyst: true +JetTrkSyst.InputFile: JetTrackSyst.root diff --git a/Reconstruction/MET/METUtilities/data/data12_8TeV/final2/TrackSoftTerms.config b/Reconstruction/MET/METUtilities/data/data12_8TeV/final2/TrackSoftTerms.config new file mode 100644 index 0000000000000000000000000000000000000000..dfcaad759f130155b737113b85d85c5c51d62e2c --- /dev/null +++ b/Reconstruction/MET/METUtilities/data/data12_8TeV/final2/TrackSoftTerms.config @@ -0,0 +1,10 @@ +############################################################################## +# Configuration file for MET_Track +############################################################################## +Conf.InputFile: METUtilities/data12_8TeV/final2/SoftTermsSyst_TrackSoftTerms.root +Conf.DoGlobalSyst: false + +PtHardSyst.sourcedir: / +PtHardSyst.suffix: TST + +Conf.DoJetTrkSyst: false diff --git a/Reconstruction/MET/METUtilities/python/METMakerConfig.py b/Reconstruction/MET/METUtilities/python/METMakerConfig.py index 1dd63dd6f142e9f2ea1f3da2aff3f70663ca09c1..400cf4384bfea6a248facee2f62d401aae7051c9 100644 --- a/Reconstruction/MET/METUtilities/python/METMakerConfig.py +++ b/Reconstruction/MET/METUtilities/python/METMakerConfig.py @@ -9,9 +9,10 @@ def getMETMakerAlg(suffix): from AthenaCommon.AlgSequence import AlgSequence topSequence = AlgSequence() + doPFlow = 'PFlow' in suffix metMaker = CfgMgr.met__METMaker('METMaker_'+suffix, - JetPtCut=0e3, - DoJetJvfCut=False); + DoPFlow=doPFlow, + JetPtCut=0e3); ToolSvc += metMaker makerAlg = CfgMgr.met__METMakerAlg('METMakerAlg_'+suffix, diff --git a/Reconstruction/MET/METUtilities/src/METMakerAlg.cxx b/Reconstruction/MET/METUtilities/src/METMakerAlg.cxx index 71e2b0bd559afa6e3179d0ed30cae89e35626d24..ceef7555afc38a0f5316135c4e067209228fd479 100644 --- a/Reconstruction/MET/METUtilities/src/METMakerAlg.cxx +++ b/Reconstruction/MET/METUtilities/src/METMakerAlg.cxx @@ -152,8 +152,6 @@ namespace met { // Electrons if(!m_eleColl.empty()) { - MissingET* metEle = new MissingET(0.,0.,0.,"RefEle",MissingETBase::Source::electron()); - newMet->push_back(metEle); uniques.clear(); ConstDataVector<ElectronContainer> metElectrons(SG::VIEW_ELEMENTS); for(const auto& el : *elCont) { @@ -161,15 +159,17 @@ namespace met { metElectrons.push_back(el); } } - ATH_CHECK( m_metmaker->rebuildMET(metEle, metElectrons.asDataVector(), metMap, uniques) ); + if( m_metmaker->rebuildMET("RefEle", xAOD::Type::Electron, newMet, + metElectrons.asDataVector(), + metMap, uniques).isFailure() ) { + ATH_MSG_WARNING("Failed to build electron term."); + } ATH_MSG_DEBUG("Selected " << metElectrons.size() << " MET electrons. " << uniques.size() << " are non-overlapping."); } // Photons if(!m_gammaColl.empty()) { - MissingET* metGamma = new MissingET(0.,0.,0.,"RefGamma",MissingETBase::Source::photon()); - newMet->push_back(metGamma); uniques.clear(); ConstDataVector<PhotonContainer> metPhotons(SG::VIEW_ELEMENTS); for(const auto& ph : *phCont) { @@ -177,15 +177,17 @@ namespace met { metPhotons.push_back(ph); } } - ATH_CHECK( m_metmaker->rebuildMET(metGamma, metPhotons.asDataVector(), metMap, uniques) ); + if( m_metmaker->rebuildMET("RefGamma", xAOD::Type::Photon, newMet, + metPhotons.asDataVector(), + metMap, uniques).isFailure() ) { + ATH_MSG_WARNING("Failed to build photon term."); + } ATH_MSG_DEBUG("Selected " << metPhotons.size() << " MET photons. " << uniques.size() << " are non-overlapping."); } // Taus if(!m_tauColl.empty()) { - MissingET* metTau = new MissingET(0.,0.,0.,"RefTau",MissingETBase::Source::tau()); - newMet->push_back(metTau); uniques.clear(); ConstDataVector<TauJetContainer> metTaus(SG::VIEW_ELEMENTS); for(const auto& tau : *tauCont) { @@ -193,15 +195,17 @@ namespace met { metTaus.push_back(tau); } } - ATH_CHECK( m_metmaker->rebuildMET(metTau, metTaus.asDataVector(), metMap, uniques) ); + if( m_metmaker->rebuildMET("RefTau", xAOD::Type::Tau, newMet, + metTaus.asDataVector(), + metMap, uniques).isFailure() ){ + ATH_MSG_WARNING("Failed to build tau term."); + } ATH_MSG_DEBUG("Selected " << metTaus.size() << " MET taus. " << uniques.size() << " are non-overlapping."); } // Muons if(!m_muonColl.empty()) { - MissingET* metMuon = new MissingET(0.,0.,0.,"Muons",MissingETBase::Source::muon()); - newMet->push_back(metMuon); uniques.clear(); ConstDataVector<MuonContainer> metMuons(SG::VIEW_ELEMENTS); for(const auto& mu : *muonCont) { @@ -209,27 +213,28 @@ namespace met { metMuons.push_back(mu); } } - ATH_CHECK( m_metmaker->rebuildMET(metMuon, metMuons.asDataVector(), metMap, uniques) ); + if( m_metmaker->rebuildMET("Muons", xAOD::Type::Muon, newMet, + metMuons.asDataVector(), + metMap, uniques).isFailure() ) { + ATH_MSG_WARNING("Failed to build muon term."); + } ATH_MSG_DEBUG("Selected " << metMuons.size() << " MET muons. " << uniques.size() << " are non-overlapping."); } - MissingET* metJet = new MissingET(0.,0.,0.,"RefJet",MissingETBase::Source::jet()); - MissingET* metSoftCl = new MissingET(0.,0.,0.,"SoftClus", - MissingETBase::Source::softEvent()|MissingETBase::Source::cluster()); - MissingET* metSoftTrk = new MissingET(0.,0.,0.,"PVSoftTrk", - MissingETBase::Source::softEvent()|MissingETBase::Source::idTrack()); - newMet->push_back(metJet); - newMet->push_back(metSoftTrk); - newMet->push_back(metSoftCl); - ATH_CHECK( m_metmaker->rebuildJetMET(metJet, jetCont, metMap, uniques, - metSoftCl, (*coreMet)["SoftClusCore"], - metSoftTrk, (*coreMet)["PVSoftTrkCore"]) ); + if( m_metmaker->rebuildJetMET("RefJet", "SoftClus", "PVSoftTrk", newMet, + jetCont, coreMet, metMap, false, uniques ).isFailure() ) { + ATH_MSG_WARNING("Failed to build jet and soft terms."); + } ATH_MSG_DEBUG("Of " << jetCont->size() << " jets, " - << uniques.size() << " are non-overlapping."); + << uniques.size() << " are non-overlapping."); - ATH_CHECK( m_metmaker->buildMETSum("FinalTrk", newMet, MissingETBase::Source::Track) ); - ATH_CHECK( m_metmaker->buildMETSum("FinalClus", newMet, MissingETBase::Source::LCTopo) ); + if( m_metmaker->buildMETSum("FinalTrk", newMet, MissingETBase::Source::Track).isFailure() ){ + ATH_MSG_WARNING("Building MET FinalTrk sum failed."); + } + if( m_metmaker->buildMETSum("FinalClus", newMet, MissingETBase::Source::LCTopo).isFailure() ) { + ATH_MSG_WARNING("Building MET FinalClus sum failed."); + } return StatusCode::SUCCESS; } diff --git a/Reconstruction/MET/METUtilities/src/components/METUtilities_entries.cxx b/Reconstruction/MET/METUtilities/src/components/METUtilities_entries.cxx index 55cd120b05fb40522d99a53e254a40312639e952..f71ec04976c5e1f60351413352dff6ebce69bfb5 100644 --- a/Reconstruction/MET/METUtilities/src/components/METUtilities_entries.cxx +++ b/Reconstruction/MET/METUtilities/src/components/METUtilities_entries.cxx @@ -3,6 +3,7 @@ // Top level tool #include "METUtilities/METMaker.h" #include "METUtilities/METRebuilder.h" +#include "METUtilities/METSystematicsTool.h" // Algs #include "../METUtilAlg.h" #include "../METMakerAlg.h" @@ -11,6 +12,7 @@ using namespace met; DECLARE_TOOL_FACTORY(METMaker) DECLARE_TOOL_FACTORY(METRebuilder) +DECLARE_TOOL_FACTORY(METSystematicsTool) // DECLARE_ALGORITHM_FACTORY(METUtilAlg) DECLARE_ALGORITHM_FACTORY(METMakerAlg) @@ -18,6 +20,7 @@ DECLARE_ALGORITHM_FACTORY(METMakerAlg) DECLARE_FACTORY_ENTRIES(METReconstruction) { DECLARE_TOOL(METMaker) DECLARE_TOOL(METRebuilder) + DECLARE_TOOL(METSystematicsTool) // DECLARE_ALGORITHM(METUtilAlg) DECLARE_ALGORITHM(METMakerAlg) diff --git a/Reconstruction/MET/METUtilities/util/ut_metSystematicTool.cxx b/Reconstruction/MET/METUtilities/util/ut_metSystematicTool.cxx new file mode 100755 index 0000000000000000000000000000000000000000..347bf73a9bac418d3dc6d33dde47b8ba686a6028 --- /dev/null +++ b/Reconstruction/MET/METUtilities/util/ut_metSystematicTool.cxx @@ -0,0 +1,563 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifdef ROOTCORE +#include "xAODRootAccess/Init.h" +#include "xAODRootAccess/TEvent.h" +#include "xAODRootAccess/TStore.h" +#endif + +// FrameWork includes +#include "AsgTools/ToolHandle.h" +#include "AsgTools/AsgTool.h" +#include "xAODBase/IParticleContainer.h" +#include "xAODMissingET/MissingETAuxContainer.h" + +#include "assert.h" +#include "TFile.h" + +#define private public //ugly hacking way to do this... +#define protected public +#include "METUtilities/METSystematicsTool.h" + + +static std::string const toolname = "tool"; + +bool testDynamicCast(){std::cout << __PRETTY_FUNCTION__ << std::endl; + + xAOD::MissingET * testMET = new xAOD::MissingET("myMET"); + + const xAOD::IParticleContainer* cont = dynamic_cast<const xAOD::IParticleContainer*>(testMET->container()); + assert(cont == nullptr ); + const xAOD::MissingETContainer* mcont = dynamic_cast<const xAOD::MissingETContainer*>(cont); + assert(mcont == nullptr); + const xAOD::MissingETContainer* cont2 = dynamic_cast<const xAOD::MissingETContainer*>(testMET->container()); + assert(cont2 == nullptr); + + delete testMET; + + return true; +} + +bool testAddObjToCont(){std::cout << __PRETTY_FUNCTION__ << std::endl; + + xAOD::MissingETContainer * testMETcont = new xAOD::MissingETContainer(); + assert(testMETcont != nullptr);//passes + + xAOD::MissingET * myObj = new xAOD::MissingET("myMETobj"); + testMETcont->push_back(myObj); + assert(testMETcont->at(0) == myObj);//passes + assert(testMETcont->size() == 1); //one element passes + + + const xAOD::MissingETContainer* cont = dynamic_cast<const xAOD::MissingETContainer*>(myObj->container()); + // dynamic_cast<const xAOD::IParticleContainer*>(myObj->container()); + assert(cont != nullptr ); + std::cout << "checked that the cont isn't empty" << std::endl; + assert(cont == testMETcont); + std::cout << "checked the cast worked, and we have the same final address as original" << std::endl; + assert(*cont == *testMETcont); + std::cout << "checked the cast worked, and we have the same final object as original" << std::endl; + + delete myObj; + + return true; +} + + +bool testDefaultHistosFilled(){std::cout << __PRETTY_FUNCTION__ << std::endl; + met::METSystematicsTool tool(toolname); + tool.msg().setLevel(MSG::VERBOSE); + assert(tool.initialize()); + + assert(tool.shiftpara_pthard_njet_mu!=nullptr); + assert(tool.resopara_pthard_njet_mu !=nullptr); + assert(tool.resoperp_pthard_njet_mu !=nullptr); + assert(tool.jet_systRpt_pt_eta ==nullptr); + + return true; +} + +bool testJetTrkHistosFilled(){std::cout << __PRETTY_FUNCTION__ << std::endl; + met::METSystematicsTool tool(toolname); + tool.msg().setLevel(MSG::VERBOSE); + tool.setProperty("ConfigJetTrkFile" ,"METTrack_2012.config"); + assert(tool.initialize()); + + assert(tool.shiftpara_pthard_njet_mu!=nullptr); + assert(tool.resopara_pthard_njet_mu !=nullptr); + assert(tool.resoperp_pthard_njet_mu !=nullptr); + assert(tool.jet_systRpt_pt_eta !=nullptr); + + return true; +} + + +bool testNoCollection(){std::cout << __PRETTY_FUNCTION__ << std::endl;//test source checking + met::METSystematicsTool tool(toolname); + tool.msg().setLevel(MSG::VERBOSE); + assert(tool.initialize()); + xAOD::MissingET softMET("softMet", MissingETBase::Source::softEvent()); + xAOD::MissingET softMETcopy(softMET); + + assert( ! tool.applyCorrection(softMET)) ; //this should fail because the softMET isn't owned by a collection, so we can't applyCorrection on it + // assert(scaledMET.mpx() != 0); + // assert(scaledMET.mpy() != 0); + // assert(scaledMET.sumet() != 0); + assert(softMET == softMETcopy); + + + return true; +} + +bool testNonSoftTermFailure(){std::cout << __PRETTY_FUNCTION__ << std::endl;//test source checking + + met::METSystematicsTool tool(toolname); + tool.msg().setLevel(MSG::VERBOSE); + assert(tool.initialize()); + + xAOD::MissingET jetMET("jetMet", MissingETBase::Source::jet()); + + assert( !tool.applyCorrection(jetMET));//we fail because our object isn't owned by a collection + //assert(scaledMET.mpx() == 0); + // assert(scaledMET.mpy() == 0); + // assert(scaledMET.sumet() == 0); + //assert(jetMET.source() == MissingETBase::Source::unknown()); + + + return true; +} + +bool testSoftTermSuccess(){std::cout << __PRETTY_FUNCTION__ << std::endl; + met::METSystematicsTool tool(toolname); + tool.msg().setLevel(MSG::VERBOSE); + assert(tool.initialize()); + + std::cout << "creating test containers" << std::endl; + xAOD::MissingETContainer * testMETcont = new xAOD::MissingETContainer(); + xAOD::MissingETAuxContainer* testMETcontAux = new xAOD::MissingETAuxContainer(); + testMETcont->setStore( testMETcontAux ); + assert(testMETcont != nullptr);//passes + + std::cout << "creating test object" << std::endl; + xAOD::MissingET * myObj = new xAOD::MissingET("myMETobj", MissingETBase::Source::softEvent()); + testMETcont->push_back(myObj); + assert(testMETcont->at(0) == myObj);//passes + assert(testMETcont->size() == 1); //one element passes + + std::cout << "applying correction" << std::endl; + assert( tool.applyCorrection(*myObj));//this time, we should pass because our object is owned + + delete testMETcont; + + return true; +} + +bool testRetrieveMETTerms(){std::cout << __PRETTY_FUNCTION__ << std::endl; + met::METSystematicsTool tool(toolname); + tool.msg().setLevel(MSG::VERBOSE); + assert(tool.initialize()); + + xAOD::MissingETContainer const * cont = nullptr; + + if( ! tool.evtStore()->retrieve(cont,"MET_RefFinal") ) return true; + + const xAOD::MissingET* refgamma = (*cont)["RefGamma"]; + + assert(refgamma != nullptr);//we can get met terms out + + return true; +} + + +bool testAddAffectingSystematic(){std::cout << __PRETTY_FUNCTION__ << std::endl; + met::METSystematicsTool tool(toolname); + tool.msg().setLevel(MSG::VERBOSE); + assert(tool.initialize()); + + CP::SystematicSet recSys = tool.recommendedSystematics(); + + std::cout << "sys set" << recSys.name() << std::endl; + + CP::SystematicVariation test("testSystematic");//recommendedSystematic + assert( !tool.isAffectedBySystematic(test)); + + std::cout << "test set " << test.name() << std::endl; + assert( tool.addAffectingSystematic(test , true /*recommended */) ); + + std::cout << "sys set 2 " << tool.recommendedSystematics().name() << std::endl; + assert( tool.isAffectedBySystematic(test)); + + + return true; +} + +bool testProjectST(){std::cout << __PRETTY_FUNCTION__ << std::endl; + met::METSystematicsTool tool(toolname); + tool.msg().setLevel(MSG::VERBOSE); + assert(tool.initialize()); + + xAOD::MissingET yaxisSoftTerm (0., 15., 50.); + xAOD::MissingET xaxisHardTerm (10., 0., 50.); + + xAOD::MissingET proj = tool.projectST(yaxisSoftTerm, xaxisHardTerm); + + std::cout << "projection "<< proj.mpx() << ", " << proj.mpy() << "," << proj.sumet() << std::endl; + assert(proj.mpx() == yaxisSoftTerm.mpx() && + proj.mpy() == (-1)*yaxisSoftTerm.mpy() && //in this case, our "projection" will do a reflection over the x (ptHard) axis + proj.sumet()== yaxisSoftTerm.sumet() + ); + + xAOD::MissingET projectBack = tool.projectST(proj, xaxisHardTerm); + + assert(projectBack == yaxisSoftTerm); + // assert(projectBack.mpx() == yaxisSoftTerm.mpx() && + // projectBack.mpy() == yaxisSoftTerm.mpy() && //reflect back + // projectBack.sumet()== yaxisSoftTerm.sumet() + // ); + + return true; + +} + +bool testProjectST2(){std::cout << __PRETTY_FUNCTION__ << std::endl; + met::METSystematicsTool tool(toolname); + tool.msg().setLevel(MSG::VERBOSE); + assert(tool.initialize()); + + xAOD::MissingET yaxisSoftTerm (0., 15., 50.); + xAOD::MissingET xaxisHardTerm (10., 0., 50.); + + xAOD::MissingET proj = tool.softTrkSyst_scale(yaxisSoftTerm, xaxisHardTerm, 0. ); + std::cout << "projection: "<< proj.mpx() << ", " << proj.mpy() << "," << proj.sumet() << std::endl; + assert(proj == yaxisSoftTerm); + // assert(proj.mpx() == yaxisSoftTerm.mpx() && + // proj.mpy() == yaxisSoftTerm.mpy() && + // proj.sumet()== yaxisSoftTerm.sumet() + // ); + + double shift = 1.; + proj = tool.softTrkSyst_scale(yaxisSoftTerm, xaxisHardTerm, shift); // + std::cout << "projection: "<< proj.mpx() << ", " << proj.mpy() << "," << proj.sumet() << std::endl; + assert(proj.mpx() == (yaxisSoftTerm.mpx()+shift) && + proj.mpy() == yaxisSoftTerm.mpy() && + proj.sumet()== yaxisSoftTerm.sumet() + ); + + return true; +} + +// bool testPtHardHistosFilled0(){std::cout << __PRETTY_FUNCTION__ << std::endl; +// met::METSystematicsTool tool(toolname); +// tool.msg().setLevel(MSG::VERBOSE); +// assert(tool.initialize()); + +// xAOD::MissingET yaxisSoftTerm (15., 0., 50.); +// xAOD::MissingET xaxisHardTerm (10., 0., 50.); + +// xAOD::MissingET proj = tool.pthardsyst_reso(yaxisSoftTerm, xaxisHardTerm, 0., 0. , 0.);//test with no shift, should end the same +// std::cout << "projection "<< proj.mpx() << ", " << proj.mpy() << "," << proj.sumet() << std::endl; +// assert(proj.mpx() == yaxisSoftTerm.mpx() && +// proj.mpy() == yaxisSoftTerm.mpy() && +// proj.sumet()== yaxisSoftTerm.sumet() +// ); + +// proj = tool.pthardsyst_reso(yaxisSoftTerm, xaxisHardTerm, 1., 0. , 0.); // +// std::cout << "projection "<< proj.mpx() << ", " << proj.mpy() << "," << proj.sumet() << std::endl; +// assert(proj.mpx() == yaxisSoftTerm.mpx() && +// proj.mpy() == yaxisSoftTerm.mpy() && +// proj.sumet()== yaxisSoftTerm.sumet() +// ); + +// return true; +// } + +bool testCorrectedCopy(){std::cout << __PRETTY_FUNCTION__ << std::endl; + met::METSystematicsTool tool(toolname); + tool.msg().setLevel(MSG::VERBOSE); + assert(tool.initialize()); + + std::cout << "creating test containers" << std::endl; + xAOD::MissingETContainer * testMETcont = new xAOD::MissingETContainer(); + xAOD::MissingETAuxContainer* testMETcontAux = new xAOD::MissingETAuxContainer(); + testMETcont->setStore( testMETcontAux ); + assert(testMETcont != nullptr);//passes + + std::cout << "creating test object" << std::endl; + xAOD::MissingET * myObj = new xAOD::MissingET("myMETobj", MissingETBase::Source::softEvent()); + testMETcont->push_back(myObj); + assert(testMETcont->at(0) == myObj); + std::cout << "checked object at array position 0, total size is " << testMETcont->size() << std::endl; + assert(testMETcont->size() == 1); //one element passes + + xAOD::MissingET * correctedCopy = nullptr; + assert( tool.correctedCopy(*myObj, correctedCopy));//this time, we should pass because our object is owned + + assert(correctedCopy != nullptr); + + delete testMETcont; + + return true; +} + +bool testRetrieveTruth(){std::cout << __PRETTY_FUNCTION__ << std::endl; + xAOD::MissingETContainer * testMETcont = new xAOD::MissingETContainer(); + assert(testMETcont != nullptr); + + std::cout << "getting missingET test" << std::endl; + + xAOD::MissingET * met = (*testMETcont)["NonInt"]; + assert(met == nullptr ); + + delete testMETcont; + + return true; +} + +bool testSoftTrkPtHardScale(){std::cout << __PRETTY_FUNCTION__ << std::endl; + met::METSystematicsTool tool(toolname); + tool.msg().setLevel(MSG::VERBOSE); + assert(tool.initialize()); + + CP::SystematicVariation softTrk_ScaleUp("MET_SoftTrk_ScaleUp"); + // assert( tool.addAffectingSystematic(softTrk_scaleUp , true /*recommended */) ); + + CP::SystematicSet applySys(softTrk_ScaleUp.name()); + // assert( tool.applySystematicVariation(softTrk_scaleUp)); + assert( tool.applySystematicVariation(applySys ) ); + + std::cout << "applied systematics are : " << tool.appliedSystematicsString() << std::endl; + + std::cout << "creating test containers" << std::endl; + xAOD::MissingETContainer * testMETcont = new xAOD::MissingETContainer(); + xAOD::MissingETAuxContainer* testMETcontAux = new xAOD::MissingETAuxContainer(); + testMETcont->setStore( testMETcontAux ); + assert(testMETcont != nullptr);//passes + + std::cout << "creating test object" << std::endl; + xAOD::MissingET * myObj = new xAOD::MissingET(10., 15., 20., "myMETobj", MissingETBase::Source::softEvent()); + testMETcont->push_back(myObj); + + xAOD::MissingET * correctedCopy = nullptr; + assert( tool.correctedCopy(*myObj, correctedCopy)); + + assert(correctedCopy != nullptr); + assert(tool.applyCorrection(*myObj)); + std::cout << "correctedCopy :" << correctedCopy->mpx() << ", " << correctedCopy->mpy() << ", " << correctedCopy->sumet() << std::endl; + std::cout << "myObj :" << myObj->mpx() << ", " << myObj->mpy() << ", " << myObj->sumet() << std::endl; + assert( (int)correctedCopy->mpx() == (int)myObj->mpx() && //todo relook at this, maybe not a great way of doing this + (int)correctedCopy->mpy() == (int)myObj->mpy() && + (int)correctedCopy->sumet() == (int)myObj->sumet()); + + delete testMETcont; + + return true; +} + +bool testSoftTrkPtHardResoPara(){std::cout << __PRETTY_FUNCTION__ << std::endl; + met::METSystematicsTool tool(toolname); + tool.msg().setLevel(MSG::VERBOSE); + assert(tool.initialize()); + + CP::SystematicVariation softTrk_ResoPara("MET_SoftTrk_ResoPara"); + // assert( tool.addAffectingSystematic(softTrk_scaleUp , true /*recommended */) ); + + CP::SystematicSet applySys(softTrk_ResoPara.name()); + // assert( tool.applySystematicVariation(softTrk_scaleUp)); + assert( tool.applySystematicVariation(applySys ) ); + + std::cout << "applied systematics are : " << tool.appliedSystematicsString() << std::endl; + + std::cout << "creating test containers" << std::endl; + xAOD::MissingETContainer * testMETcont = new xAOD::MissingETContainer(); + xAOD::MissingETAuxContainer* testMETcontAux = new xAOD::MissingETAuxContainer(); + testMETcont->setStore( testMETcontAux ); + assert(testMETcont != nullptr);//passes + + std::cout << "creating test object" << std::endl; + xAOD::MissingET * myObj = new xAOD::MissingET(10., 15., 20., "myMETobj", MissingETBase::Source::softEvent()); + testMETcont->push_back(myObj); + xAOD::MissingET * correctedCopy = nullptr; + tool.setRandomSeed(199); + assert( tool.correctedCopy(*myObj, correctedCopy)); + assert(correctedCopy != nullptr); + + tool.setRandomSeed(199);//reset the random seed to redo the same calculation + assert(tool.applyCorrection(*myObj)); + + std::cout << "correctedCopy :" << correctedCopy->mpx() << ", " << correctedCopy->mpy() << ", " << correctedCopy->sumet() << std::endl; + std::cout << "myObj :" << myObj->mpx() << ", " << myObj->mpy() << ", " << myObj->sumet() << std::endl; + assert( (int)correctedCopy->mpx() == (int)myObj->mpx() && //todo relook at this, maybe not a great way of doing this + (int)correctedCopy->mpy() == (int)myObj->mpy() && + (int)correctedCopy->sumet() == (int)myObj->sumet()); + + delete testMETcont; + + return true; +} + + +bool testSoftTrkPtHardResoPerp(){std::cout << __PRETTY_FUNCTION__ << std::endl; + met::METSystematicsTool tool(toolname); + tool.msg().setLevel(MSG::VERBOSE); + assert(tool.initialize()); + + CP::SystematicVariation softTrk_ResoPerp("MET_SoftTrk_ResoPerp"); + // assert( tool.addAffectingSystematic(softTrk_scaleUp , true /*recommended */) ); + + CP::SystematicSet applySys(softTrk_ResoPerp.name()); + // assert( tool.applySystematicVariation(softTrk_scaleUp)); + assert( tool.applySystematicVariation(applySys ) ); + + std::cout << "applied systematics are : " << tool.appliedSystematicsString() << std::endl; + + std::cout << "creating test containers" << std::endl; + xAOD::MissingETContainer * testMETcont = new xAOD::MissingETContainer(); + xAOD::MissingETAuxContainer* testMETcontAux = new xAOD::MissingETAuxContainer(); + testMETcont->setStore( testMETcontAux ); + assert(testMETcont != nullptr);//passes + + std::cout << "creating test object" << std::endl; + xAOD::MissingET * myObj = new xAOD::MissingET(10., 15., 20., "myMETobj", MissingETBase::Source::softEvent()); + testMETcont->push_back(myObj); + xAOD::MissingET * correctedCopy = nullptr; + std::cout << "setting random seed" << std::endl; + tool.setRandomSeed(199); + assert( tool.correctedCopy(*myObj, correctedCopy)); + assert(correctedCopy != nullptr); + + tool.setRandomSeed(199);//reset the random seed to redo the same calculation + assert(tool.applyCorrection(*myObj)); + + std::cout << "correctedCopy :" << correctedCopy->mpx() << ", " << correctedCopy->mpy() << ", " << correctedCopy->sumet() << std::endl; + std::cout << "myObj :" << myObj->mpx() << ", " << myObj->mpy() << ", " << myObj->sumet() << std::endl; + assert( (int)correctedCopy->mpx() == (int)myObj->mpx() && //todo relook at this, maybe not a great way of doing this + (int)correctedCopy->mpy() == (int)myObj->mpy() && + (int)correctedCopy->sumet() == (int)myObj->sumet()); + + delete testMETcont; + + return true; +} + + + +bool testSoftCaloPtScale(){std::cout << __PRETTY_FUNCTION__ << std::endl; + met::METSystematicsTool tool(toolname); + tool.msg().setLevel(MSG::VERBOSE); + assert(tool.initialize()); + + CP::SystematicVariation softCalo_scaleUp("MET_SoftCalo_ScaleUp"); + // assert( tool.addAffectingSystematic(softCalo_scaleUp , true /*recommended */) ); + + CP::SystematicSet recSys(softCalo_scaleUp.name()); + // assert( tool.applySystematicVariation(softCalo_scaleUp)); + assert( tool.applySystematicVariation(recSys) ); + + std::cout << "applied systematics are : " << tool.appliedSystematicsString() << std::endl; + + std::cout << "creating test containers" << std::endl; + xAOD::MissingETContainer * testMETcont = new xAOD::MissingETContainer(); + xAOD::MissingETAuxContainer* testMETcontAux = new xAOD::MissingETAuxContainer(); + testMETcont->setStore( testMETcontAux ); + assert(testMETcont != nullptr);//passes + + std::cout << "creating test object" << std::endl; + xAOD::MissingET * myObj = new xAOD::MissingET(10., 15., 20., "myMETobj", MissingETBase::Source::softEvent()); +// xAOD::MissingET * myObj = new xAOD::MissingET("myMETobj", MissingETBase::Source::softEvent()); + testMETcont->push_back(myObj); + + xAOD::MissingET * correctedCopy = nullptr; + assert( tool.correctedCopy(*myObj, correctedCopy)); + + assert(correctedCopy != nullptr); + assert(tool.applyCorrection(*myObj)); + std::cout << "correctedCopy :" << correctedCopy->mpx() << ", " << correctedCopy->mpy() << ", " << correctedCopy->sumet() << std::endl; + std::cout << "myObj :" << myObj->mpx() << ", " << myObj->mpy() << ", " << myObj->sumet() << std::endl; + assert( (int)correctedCopy->mpx() == (int)myObj->mpx() && //todo relook at this, maybe not a great way of doing this + (int)correctedCopy->mpy() == (int)myObj->mpy() && + (int)correctedCopy->sumet() == (int)myObj->sumet()); + + delete testMETcont; + + return true; +} + +bool testJetTrackMET(){std::cout << __PRETTY_FUNCTION__ << std::endl; + met::METSystematicsTool tool(toolname); + tool.msg().setLevel(MSG::VERBOSE); + tool.setProperty("ConfigJetTrkFile" ,"METTrack_2012.config"); + + assert(tool.initialize()); + + + + return true; +} + + +bool testNoConfigFiles(){std::cout << __PRETTY_FUNCTION__ << std::endl; + met::METSystematicsTool tool(toolname); + tool.msg().setLevel(MSG::VERBOSE); + tool.setProperty("ConfigSoftTrkFile",""); + tool.setProperty("ConfigJetTrkFile",""); + tool.setProperty("ConfigSoftCaloFile",""); + + assert(tool.initialize()); + + return true; +} + +int main(){std::cout << __PRETTY_FUNCTION__ << std::endl; + // Initialise the application: + xAOD::Init() ; + + // CP::CorrectionCode::enableFailure(); + // StatusCode::enableFailure(); + // xAOD::TReturnCode::enableFailure(); + TString const fileName = "AOD.pool.root"; + // Info( APP_NAME, "Opening file: %s", fileName.Data() ); + std::auto_ptr< TFile > ifile( TFile::Open( fileName, "READ" ) ); + assert( ifile.get() ); + + // Create a TEvent object: + // xAOD::TEvent event( xAOD::TEvent::kClassAccess ); + xAOD::TEvent * event; + event = new xAOD::TEvent( xAOD::TEvent::kClassAccess ); + assert( event->readFrom( ifile.get() ) ); + // event = wk()->xaodEvent(); + + // Info( APP_NAME, "Number of events in the file: %i", + // static_cast< int >( event.getEntries() ) ); + + // Create a transient object store. Needed for the tools. + xAOD::TStore store; + + assert(testAddObjToCont()); + assert(testRetrieveTruth()); + assert(testDynamicCast()); + assert(testDefaultHistosFilled()); + assert(testJetTrkHistosFilled()); + // assert(testNoCollection()); + assert(testNonSoftTermFailure()); + assert(testSoftTermSuccess());//todo readd this one? + assert(testRetrieveMETTerms()); + assert(testAddAffectingSystematic()); + assert(testProjectST()); + assert(testProjectST2()); + assert(testCorrectedCopy()); + + assert(testSoftTrkPtHardScale()); + assert(testSoftTrkPtHardResoPerp()); + assert(testSoftTrkPtHardResoPara()); + assert(testSoftCaloPtScale()); + + assert(testJetTrackMET()); + assert(testNoConfigFiles()); + + delete event; + + return 0; + + }