From 0c13f82f8ff17bc66fa5837cdf181c98a9e23a90 Mon Sep 17 00:00:00 2001
From: Pierre-Antoine Delsart <delsart@in2p3.fr>
Date: Thu, 11 Sep 2014 12:17:33 +0200
Subject: [PATCH] 'fix checkreq' (JetMomentTools-00-02-27)
---
.../JetMomentTools/JetBadChanCorrTool.h | 185 +++++++
.../JetMomentTools/JetCaloEnergies.h | 27 +
.../JetMomentTools/JetCaloQualityTool.h | 66 +++
.../JetMomentTools/JetIsolationTool.h | 147 ++++++
.../JetMomentTools/JetLArHVTool.h | 58 +++
.../JetMuonSegmentMomentsTool.h | 38 ++
.../JetMomentTools/JetPtAssociationTool.h | 66 +++
.../JetMomentTools/JetTrackMomentsTool.h | 74 +++
.../JetMomentTools/JetVertexFractionTool.h | 70 +++
.../JetMomentTools/JetWidthTool.h | 34 ++
.../Root/JetMuonSegmentMomentsTool.cxx | 31 ++
.../Root/JetPtAssociationTool.cxx | 123 +++++
.../Root/JetTrackMomentsTool.cxx | 161 ++++++
.../Root/JetVertexFractionTool.cxx | 115 +++++
.../Jet/JetMomentTools/Root/JetWidthTool.cxx | 53 ++
.../Jet/JetMomentTools/cmt/Makefile.RootCore | 24 +
.../Jet/JetMomentTools/cmt/requirements | 31 ++
.../JetMomentTools/python/GhostAssociation.py | 148 ++++++
.../python/JetCopyMomentsAlg.py | 72 +++
.../python/JetMomentsConfigHelpers.py | 363 ++++++++++++++
.../python/SetupJetMomentTools.py | 302 ++++++++++++
.../JetMomentTools/src/JetBadChanCorrTool.cxx | 461 ++++++++++++++++++
.../JetMomentTools/src/JetCaloEnergies.cxx | 57 +++
.../JetMomentTools/src/JetCaloQualityTool.cxx | 177 +++++++
.../JetMomentTools/src/JetIsolationTool.cxx | 456 +++++++++++++++++
.../Jet/JetMomentTools/src/JetLArHVTool.cxx | 105 ++++
.../src/components/JetMomentTools_entries.cxx | 38 ++
.../src/components/JetMomentTools_load.cxx | 5 +
28 files changed, 3487 insertions(+)
create mode 100644 Reconstruction/Jet/JetMomentTools/JetMomentTools/JetBadChanCorrTool.h
create mode 100644 Reconstruction/Jet/JetMomentTools/JetMomentTools/JetCaloEnergies.h
create mode 100644 Reconstruction/Jet/JetMomentTools/JetMomentTools/JetCaloQualityTool.h
create mode 100644 Reconstruction/Jet/JetMomentTools/JetMomentTools/JetIsolationTool.h
create mode 100644 Reconstruction/Jet/JetMomentTools/JetMomentTools/JetLArHVTool.h
create mode 100644 Reconstruction/Jet/JetMomentTools/JetMomentTools/JetMuonSegmentMomentsTool.h
create mode 100644 Reconstruction/Jet/JetMomentTools/JetMomentTools/JetPtAssociationTool.h
create mode 100644 Reconstruction/Jet/JetMomentTools/JetMomentTools/JetTrackMomentsTool.h
create mode 100644 Reconstruction/Jet/JetMomentTools/JetMomentTools/JetVertexFractionTool.h
create mode 100644 Reconstruction/Jet/JetMomentTools/JetMomentTools/JetWidthTool.h
create mode 100644 Reconstruction/Jet/JetMomentTools/Root/JetMuonSegmentMomentsTool.cxx
create mode 100644 Reconstruction/Jet/JetMomentTools/Root/JetPtAssociationTool.cxx
create mode 100644 Reconstruction/Jet/JetMomentTools/Root/JetTrackMomentsTool.cxx
create mode 100644 Reconstruction/Jet/JetMomentTools/Root/JetVertexFractionTool.cxx
create mode 100644 Reconstruction/Jet/JetMomentTools/Root/JetWidthTool.cxx
create mode 100644 Reconstruction/Jet/JetMomentTools/cmt/Makefile.RootCore
create mode 100755 Reconstruction/Jet/JetMomentTools/cmt/requirements
create mode 100644 Reconstruction/Jet/JetMomentTools/python/GhostAssociation.py
create mode 100644 Reconstruction/Jet/JetMomentTools/python/JetCopyMomentsAlg.py
create mode 100644 Reconstruction/Jet/JetMomentTools/python/JetMomentsConfigHelpers.py
create mode 100644 Reconstruction/Jet/JetMomentTools/python/SetupJetMomentTools.py
create mode 100644 Reconstruction/Jet/JetMomentTools/src/JetBadChanCorrTool.cxx
create mode 100644 Reconstruction/Jet/JetMomentTools/src/JetCaloEnergies.cxx
create mode 100644 Reconstruction/Jet/JetMomentTools/src/JetCaloQualityTool.cxx
create mode 100644 Reconstruction/Jet/JetMomentTools/src/JetIsolationTool.cxx
create mode 100644 Reconstruction/Jet/JetMomentTools/src/JetLArHVTool.cxx
create mode 100644 Reconstruction/Jet/JetMomentTools/src/components/JetMomentTools_entries.cxx
create mode 100644 Reconstruction/Jet/JetMomentTools/src/components/JetMomentTools_load.cxx
diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetBadChanCorrTool.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetBadChanCorrTool.h
new file mode 100644
index 00000000000..e0abbe5288b
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetBadChanCorrTool.h
@@ -0,0 +1,185 @@
+// this file is -*- C++ -*-
+
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef JETMOMENTTOOLS_JETBADCHANCORRTOOL_H
+#define JETMOMENTTOOLS_JETBADCHANCORRTOOL_H
+
+/** @class JetBadChanCorrTol.
+ Calculates correction for dead cells
+
+ It produces and stores three jet moment : BchCorrCell, BchCorrDotx, BchCorrJet
+ after computing the energy fraction of dead cells.
+
+ IMPORTANT : this run2 implementation is less complete than before.
+ -> it can not make use of the older IHadronicCalibrationTool (wasn't used by default anyway)
+ and therefore doesn't take the cryostat correction into account.
+
+ Also it now requires the existence of a MissingCaloCellsMap in SG.
+
+ Lots of commented blocks are still around until we're sure they're not needed anymore.
+
+ @author Yousuke Kataoka
+ @date (first implementation) December , 2009
+ @author P.A. Delsart
+ @date (faster implementation) March , 2011
+ @date (Run 2 re-implemantation) March , 2013
+
+*/
+
+
+#include "JetRec/JetModifierBase.h"
+
+#include "JetRecTools/MissingCellListTool.h"
+
+#include "CaloIdentifier/CaloCell_ID.h"
+
+#include "TH1D.h"
+
+
+class Identifier;
+class CaloDetDescrManager;
+
+class ITileBadChanTool;
+class ILArBadChanTool;
+class ITHistSvc;
+
+class JetBadChanCorrTool: public JetModifierBase
+{
+ ASG_TOOL_CLASS0(JetBadChanCorrTool);
+public:
+ JetBadChanCorrTool(const std::string& name);
+
+ virtual ~JetBadChanCorrTool();
+
+ virtual StatusCode initialize();
+
+ virtual int modifyJet( xAOD::Jet& pJet) const ;
+
+ virtual StatusCode setupEvent();
+
+
+protected:
+ int correctionFromCellsInJet( xAOD::Jet* j, const jet::CaloCellFastMap * badCellMap) const;
+ int correctionFromCellsInCone( xAOD::Jet* j, const jet::CaloCellFastMap * badCellMap) const;
+
+ int correctionFromClustersBadCells( xAOD::Jet* j) const;
+
+ private:
+
+ ServiceHandle<ITHistSvc> m_thistSvc;
+
+ // const CaloDetDescrManager* m_caloDDM;
+ // const CaloCell_ID* m_calo_id;
+
+ // limit to calculate moments
+ int m_nBadCellLimit;
+
+
+
+
+
+ // for jet-level correction
+ std::string m_streamName;
+ std::string m_profileName;
+ std::string m_profileTag;
+
+ bool m_useCone;
+ // double m_coneDr;
+
+ // for cell-level correction
+ bool m_useCalibScale;
+ //ToolHandle<IHadronicCalibrationTool> m_calibTool;
+
+ //ToolHandle<IMissingCellListTool> m_missingCellToolHandle;
+ //MissingCellListTool* m_missingCellTool;
+ std::string m_missingCellMapName;
+
+ bool m_forceMissingCellCheck;
+
+ bool m_useClusters;
+
+ // jet profiles
+ class ProfileData {
+ public:
+ ProfileData(TH1D* th, int sample,
+ double ptMin=0, double ptMax=9999,
+ double etaMin=0, double etaMax=5.0,
+ double phiMin=-M_PI, double phiMax=M_PI):
+ th(th), sample(sample),
+ ptMin(ptMin), ptMax(ptMax),
+ etaMin(etaMin), etaMax(etaMax),
+ phiMin(phiMin), phiMax(phiMax) {}
+
+ virtual ~ProfileData(){}
+
+ bool match(double pt, int sample, double eta, double phi) const {
+ return ( pt>=ptMin && pt<ptMax
+ && sample==this->sample
+ && fabs(eta)>=etaMin && fabs(eta)<etaMax
+ && phi>=phiMin && phi<phiMax);
+ }
+
+ double frac(double dr) const {
+ int idr = th->FindBin(dr);
+ return th->GetBinContent(idr);
+ }
+ private:
+ TH1D* th;
+ int sample;
+ double ptMin;
+ double ptMax;
+ double etaMin;
+ double etaMax;
+ double phiMin;
+ double phiMax;
+ };
+ std::vector<ProfileData> m_profileDatas[CaloCell_ID::Unknown];//24
+
+ double getProfile(double pt, double dr, int sample, double eta, double phi) const;
+
+};
+#endif
+// DoxygenDocumentation
+/*! @class JetBadChanCorrTool
+ @brief JetBadChanCorrTool
+
+The energy fraction of dead cells are estimated by jet-level estimator using jet profile.
+Some cells are already corrected in cell-level using neighboring cells, so this fraction is calculated.
+The estimation by jet-level for those cells corrected in cell-level is also stored to see the difference or combined the corrections.
+
+BCH_CORR_JET= jet level estimation of the energy fraction of dead cells
+BCH_CORR_CELL= cell level correction, which is already corrected
+BCH_CORR_JET_FORCELL= jet level estimation for the corrected cells in cell level
+
+The relation and usage of the moment is
+jet = jet already corrected in cell level
+jet x ( 1 - BCH_CORR_CELL ) = jet without cell level correction
+jet x ( 1 - BCH_CORR_CELL + BCH_CORR_JET ) = jet corrected by jet level estimation
+jet x ( 1 + BCH_CORR_JET - BCH_CORR_JET_FORCELL ) = jet with cell level correction and corrected in jet level for not covered cells
+
+
+<b>Properties:</b>
+<table style="text-align: left; width: 80%;" border="1">
+<tr><td> <b>Name</b> </td> <td><b>Type</b> </td> <td><b>Default</b> </td> <td><b>Description</b></td></tr>
+
+<tr><td>ProcessNJetMax </td> <td> </td> <td>999999 </td> <td></td></tr>
+<tr><td>NBadCellLimit </td> <td> </td> <td>10000 </td> <td>limit to calculate the moments, if exceed, giving up the calculation ie) when Calorimeter off</td></tr>
+<tr><td>AttachCorrCell </td> <td> </td> <td>True </td> <td>switch to attach BCH_CORR_CELL</td></tr>
+<tr><td>AttachCorrDOTX </td> <td> </td> <td>True </td> <td>switch to attach BCH_CORR_DOTX</td></tr>
+<tr><td>AttachCorrJet </td> <td> </td> <td>True </td> <td>switch to attach BCH_CORR_JET</td></tr>
+<tr><td>AttachCorrJetForCell </td> <td> </td> <td>True </td> <td>switch to attach BCH_CORR_JET_FORCELL</td></tr>
+
+<tr><td>StreamName </td> <td> </td> <td>"/JetBadChanCorrTool/" </td> <td>name of stream name</td></tr>
+<tr><td>ProfileName </td> <td> </td> <td>"JetBadChanCorrTool.root" </td> <td>name of profile data</td></tr>
+<tr><td>ProfileTag </td> <td> </td> <td>"" </td> <td>tag to select profile,"" for all in file</td></tr>
+<tr><td>UseCone </td> <td> </td> <td>True </td> <td>boundary of correction, cone or only associated cells</td></tr>
+<tr><td>ConeDr </td> <td> </td> <td>0.4 </td> <td>boundary of correction in case UseCone=True</td></tr>
+<tr><td>UseCalibScale </td> <td> </td> <td>False </td> <td>scale for calculation of the contribution of bad cells</td></tr>
+<tr><td>CellCalibrator </td> <td> </td> <td> </td> <td>calibration tool for cell-level correction</td></tr>
+
+</table>
+*/
+
diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetCaloEnergies.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetCaloEnergies.h
new file mode 100644
index 00000000000..c7ee80fc727
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetCaloEnergies.h
@@ -0,0 +1,27 @@
+// this file is -*- C++ -*-
+
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef JETMOMENTTOOLS_JETCALOENERGIES_H
+#define JETMOMENTTOOLS_JETCALOENERGIES_H
+
+#include "JetRec/JetModifierBase.h"
+
+class JetCaloEnergies : public JetModifierBase {
+ ASG_TOOL_CLASS0(JetCaloEnergies);
+public:
+
+ JetCaloEnergies(const std::string & t);
+
+
+ virtual int modifyJet(xAOD::Jet& ) const ;
+
+
+};
+
+
+#undef ASG_DERIVED_TOOL_CLASS
+#endif
+
diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetCaloQualityTool.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetCaloQualityTool.h
new file mode 100644
index 00000000000..cfeab4afea8
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetCaloQualityTool.h
@@ -0,0 +1,66 @@
+// this file is -*- C++ -*-
+
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+
+/**
+ @class JetCaloQualityTool
+ Calculates calorimeter based variables for jet quality
+ @author Nikola Makovec
+ @author P-A Delsart
+ @date (first implementation) October , 2009
+ @date (run 2 implementation) February , 2014
+*/
+
+#ifndef JETREC_JETCALOQUALITYTOOL_H
+#define JETREC_JETCALOQUALITYTOOL_H
+
+
+#include "JetRec/JetModifierBase.h"
+
+#include "JetUtils/JetCaloCalculations.h"
+
+#include <vector>
+#include <string>
+
+
+class JetCaloQualityTool: public JetModifierBase {
+ ASG_TOOL_CLASS0(JetCaloQualityTool);
+
+public:
+ JetCaloQualityTool(const std::string & name);
+
+ virtual int modifyJet(xAOD::Jet& ) const ;
+
+ virtual StatusCode initialize();
+
+ private:
+
+ int m_LArQualityCut;
+ int m_TileQualityCut;
+
+ bool m_doFracSamplingMax;
+
+ bool m_doN90;
+ bool m_doLArQuality;
+ bool m_doTileQuality;
+ bool m_doTiming;
+ bool m_doHECQuality;
+ bool m_doNegativeE;
+ bool m_doAverageLArQF;
+ bool m_timingOnlyPosE;
+ bool m_computeFromCluster;
+ bool m_doJetCentroid;
+
+
+ std::vector <double> m_timingCellTimeCuts;
+ std::vector <double> m_timingClusterTimeCuts;
+
+ mutable jet::JetCaloCalculations m_jetCalculations;
+
+
+};
+#endif
+
diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetIsolationTool.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetIsolationTool.h
new file mode 100644
index 00000000000..f1bbd89509f
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetIsolationTool.h
@@ -0,0 +1,147 @@
+///////////////////////// -*- C++ -*- /////////////////////////////
+
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// JetIsolationTool.h
+// Header file for class JetIsolationTool
+// Author: P-A. Delsart
+///////////////////////////////////////////////////////////////////
+#ifndef JETMOMENTTOOLS_JETISOLATIONTOOL_H
+#define JETMOMENTTOOLS_JETISOLATIONTOOL_H 1
+
+#include "JetRec/JetModifierBase.h"
+#include "JetUtils/TiledEtaPhiMap.h"
+#include "AsgTools/CLASS_DEF.h"
+
+///////////////////////////////////////////////////////////////////
+/// \class JetIsolationTool
+///
+/// This JetModifier calculates isolation variables for jets.
+///
+/// The isolation variables are calculated from a list of input particles
+/// which are close to the jet but not part of its constituents.
+/// There are multiple options to find this list and then to calculate the
+/// variables.
+/// They are passed to the tools by giving a list of string identifiers through
+/// the IsolationCalculations property as in
+///
+/// IsolationCalculations = ["IsoKR:11:Perp","IsoKR:11:Par", "IsoFixedCone:10:SumPt",]
+///
+/// where the strings have the form "ISOCRITERIA:NN:VAR"
+/// - ISOCRITERIA : encode the method used to find isolation particles (IsoKR, IsoDelta, IsoFixedCone, IsoFixedArea, Iso6To8)
+/// - NN : encode the value of the main parameter of ISOCRITERIA : the actual parameter used is NN/10. (same convetion as jet radius in JetContainer names)
+/// - VAR : encode the variable to be calculated (Per, Par, SumPt, P)
+///
+/// The input particles container from which constituents come from must be specified
+/// in the ConstituentContainer property. This must correspond to an IParticleContainer in the event store.
+/// (longer term : make it possible to pass PseudoJetContainer)
+///
+///
+/// There are several technical difficulties
+/// - the multiplicity of possible variables
+/// - the dependency on the calibration state of the constituents
+/// - the calculation time (looking up the full input container for each jet constituents can be *very* slow)
+/// These are solved at the cost of increased multiplicity (internal helper classes, usage of a special fast lookup map)
+///
+/// WARNING : currently works well only for LCTopoJets, TrackJets, TruthJets AND small R jets (R ~<0.8)
+///
+///////////////////////////////////////////////////////////////////////
+
+namespace jet {
+ namespace JetIsolation {
+ class IsolationCalculator;
+ }
+
+ /// \class ParticlePosition
+ /// Object describing the position of a particle in (eta,phi) and usable within a
+ /// TiledEtaPhiMap (see below).
+ struct ParticlePosition {
+ ParticlePosition(const xAOD::IParticle* p, int sigstate=-1): m_eta(p->eta()), m_phi(p->phi()) ,m_part(p), m_sigState(sigstate){};
+ ParticlePosition(double x=0, double y=0):m_eta(x),m_phi(y),m_part(NULL), m_sigState(-1){}
+ double x() const {return m_eta;}
+ double y() const {return m_phi;}
+ void setX(double x){m_eta=x;}
+ void setY(double x){m_phi=x;}
+
+ const xAOD::IParticle* particle()const {return m_part;}
+ void setParticle(const xAOD::IParticle* p){m_part=p;}
+
+ struct DR2 {
+ double operator()(const ParticlePosition &p1,const ParticlePosition &p2) const {
+ return JetTiledMap::utils::DR2(p1.x(),p1.y(), p2.x(), p2.y() );
+ }
+ };
+
+
+ protected:
+ double m_eta,m_phi;
+ const xAOD::IParticle* m_part ;
+ int m_sigState; // if any
+ };
+
+ /// This map is a geometric helper. It greatly speeds up the retrieval
+ /// of a set of point around a given position in (eta,phi).
+ /// It is saved in the event store, so it's not rebuild for every JetContainer
+ /// (could be factorized in a separate AlgTool if needed)
+ typedef JetTiledMap::TiledEtaPhiMap<ParticlePosition> ParticleFastMap;
+
+}
+
+
+class JetIsolationTool
+ : public JetModifierBase
+{
+ ASG_TOOL_CLASS0(JetIsolationTool);
+
+ ///////////////////////////////////////////////////////////////////
+ // Public methods:
+ ///////////////////////////////////////////////////////////////////
+ public:
+ friend class jet::JetIsolation::IsolationCalculator;
+
+ // Copy constructor:
+
+ /// Constructor with parameters:
+ JetIsolationTool(const std::string &myname);
+
+ /// Destructor:
+ virtual ~JetIsolationTool();
+
+ // Athena algtool's Hooks
+ virtual StatusCode initialize();
+ virtual StatusCode finalize();
+
+ virtual int modify(xAOD::JetContainer& jets) const;
+ virtual int modifyJet(xAOD::Jet& ) const{return 0;};
+
+
+
+ ///////////////////////////////////////////////////////////////////
+
+
+protected:
+
+
+ const jet::ParticleFastMap * retrieveInputMap() const ;
+
+ private:
+
+
+ std::vector<std::string> m_isolationCodes;
+ float m_deltaRmax;
+ std::string m_inputPseudoJets;
+
+ std::vector<jet::JetIsolation::IsolationCalculator*> m_isoCalculators;
+
+
+ bool m_rejectNegE;
+
+
+};
+
+
+CLASS_DEF(jet::ParticleFastMap , 106240765 , 1)
+
+#endif //> !JETMOMENTTOOLS_JETISOLATIONTOOL_H
diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetLArHVTool.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetLArHVTool.h
new file mode 100644
index 00000000000..a9abdf5dd41
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetLArHVTool.h
@@ -0,0 +1,58 @@
+// this file is -*- C++ -*-
+
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+
+/**
+ @class JetLArHVTool
+ Calculates calorimeter based variables for jet quality
+ @author Nikola Makovec
+ @author P-A Delsart
+ @date (first implementation) October , 2009
+ @date (run 2 implementation) February , 2014
+*/
+
+#ifndef JETMOMENTTOOLS_JETLARHVTOOL_H
+#define JETMOMENTTOOLS_JETLARHVTOOL_H
+
+#include "JetRec/JetModifierBase.h"
+#include "AsgTools/ToolHandle.h"
+
+
+#include "LArElecCalib/ILArHVCorrTool.h"
+
+class JetLArHVTool: public JetModifierBase {
+ ASG_TOOL_CLASS0(JetLArHVTool);
+
+public:
+ JetLArHVTool(const std::string & name);
+
+ virtual int modifyJet(xAOD::Jet& ) const ;
+
+ virtual StatusCode initialize();
+
+ private:
+
+ bool m_useCells;
+ std::string m_channelKey;
+ std::string m_keyHVScaleCorr;
+
+#ifdef ASGTOOL_ATHENA
+ ToolHandle<ILArHVCorrTool> m_hvCorrTool;
+#endif
+
+};
+
+#endif
+// DoxygenDocumentation
+/*! @class JetLArHVMoment
+ @brief JetLArHVMoment
+
+
+Calculates the proportion of cells in a bad LAr HV area, and gives an information about the energy affected.
+(the second moment should change later).
+
+</table>
+*/
diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetMuonSegmentMomentsTool.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetMuonSegmentMomentsTool.h
new file mode 100644
index 00000000000..b8de3c3ff6a
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetMuonSegmentMomentsTool.h
@@ -0,0 +1,38 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// JetMuonSegmentMomentsTool.h
+
+#ifndef JETMOMENTTOOLS_JETMUONSEGMENTMOMENTSTOOL_H
+#define JETMOMENTTOOLS_JETMUONSEGMENTMOMENTSTOOL_H
+
+/// Steven Schramm \n
+/// February 2014
+///
+/// Tool to calculate general muon segment based jet moments
+/// Primarily intended for punch-through studies and corrections
+///
+/// Moments to calculate:
+/// None yet...
+
+#include "JetRec/JetModifierBase.h"
+
+class JetMuonSegmentMomentsTool : public JetModifierBase
+{
+ ASG_TOOL_CLASS(JetMuonSegmentMomentsTool,IJetModifier);
+
+ public:
+ // Constructor from tool name
+ JetMuonSegmentMomentsTool(const std::string& name);
+
+ // Inherited methods to modify a jet
+ virtual int modifyJet(xAOD::Jet& jet) const;
+
+ private:
+ // Configurable parameters
+ std::string m_assocMuonSegName;
+};
+
+#endif
+
diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetPtAssociationTool.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetPtAssociationTool.h
new file mode 100644
index 00000000000..c339a3d15dd
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetPtAssociationTool.h
@@ -0,0 +1,66 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// JetPtAssociationTool.h
+
+#ifndef JetMomentTools_JetPtAssociationTool_H
+#define JetMomentTools_JetPtAssociationTool_H
+
+/// David Adams \n
+/// April 2014
+///
+/// Tool to associate jets using the fraction of AP (associated particle) pT
+/// shared with the constituents of jets in a matching container.
+/// Properties:
+/// AssociationName - Name of the association
+/// InputContainer - Name of the matching jet container
+
+#include "JetRec/JetModifierBase.h"
+#include "xAODJet/JetContainer.h"
+
+class JetPtAssociationTool : public JetModifierBase {
+ ASG_TOOL_CLASS(JetPtAssociationTool, IJetModifier);
+
+public:
+
+ typedef std::vector<const xAOD::IParticle*> APVector;
+ typedef std::vector<unsigned> IndexVector;
+ typedef std::vector<float> FloatVector;
+
+public:
+
+ /// Constructor from tool name.
+ JetPtAssociationTool(std::string myname);
+
+ /// Initialization.
+ StatusCode initialize();
+
+ /// Inherited method to modify a jet.
+ /// Extract and record the element link and pT fraction
+ /// for the jet with the highest such fraction.
+ int modifyJet(xAOD::Jet& jet) const;
+
+ /// Return the matched pT sum for each jet in a collection.
+ /// The matching is done with match().
+ /// aps - input AP vector
+ /// jet - Jet to be matched
+ /// ptsums - output vector with matched pT fraction for each jet
+ int ptfrac(const APVector& aps, const xAOD::JetContainer& jets, FloatVector& ptsums) const;
+
+ /// Return the vector of AP's that appear as
+ /// constituents of the given jet.
+ /// aps - input AP vector
+ /// jet - Jet to be matched
+ /// apvs - output vector of matched AP indices
+ int match(const APVector& aps, const xAOD::Jet& jet, APVector& apvs) const;
+
+private: // data
+
+ /// Properties.
+ std::string m_aname;
+ std::string m_conname;
+
+};
+
+#endif
diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetTrackMomentsTool.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetTrackMomentsTool.h
new file mode 100644
index 00000000000..7212e47f0cc
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetTrackMomentsTool.h
@@ -0,0 +1,74 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// JetTrackMomentsTool.h
+
+#ifndef JETMOMENTTOOLS_JETTRACKMOMENTSTOOL_H
+#define JETMOMENTTOOLS_JETTRACKMOMENTSTOOL_H
+
+/// Steven Schramm \n
+/// February 2014
+///
+/// Tool to calculate general track-based jet moments
+///
+/// Moments to calculate:
+/// NumTrkPtXXX
+/// SumPtTrkPtXXX
+/// TrackWidthPtXXX
+/// where the track pT threhsold is XXX MeV.
+
+#include "JetRec/JetModifierBase.h"
+#include "JetEDM/TrackVertexAssociation.h"
+
+#include "xAODTracking/TrackParticle.h"
+#include "xAODTracking/Vertex.h"
+#include "xAODTracking/VertexContainer.h"
+
+#include <vector>
+#include <string>
+
+
+class JetTrackMomentsTool : public JetModifierBase {
+ ASG_TOOL_CLASS(JetTrackMomentsTool,IJetModifier);
+
+public:
+
+ // Constructor from tool name
+ JetTrackMomentsTool(const std::string& name);
+
+ // Inherited methods to modify a jet
+ // Calls getTrackMoments and puts the results in the jet
+ virtual int modifyJet(xAOD::Jet& jet) const;
+
+private:
+
+ // Configurable parameters
+ std::string m_vertexContainer;
+ std::string m_assocTracksName;
+ std::string m_tva;
+ std::vector<float> m_minTrackPt;
+
+ // Private struct to make it unambiguous what each value is (rather than a vector)
+ // Doubles for calculation for now - will be written as float in the aux store
+ struct TrackMomentStruct { int numTrk; double sumPtTrk; double trackWidth; };
+
+ // Local method to calculate NumTrk, SumPtTrk, and TrackWidth for all vertices
+ const std::vector<TrackMomentStruct>
+ getTrackMoments(const xAOD::Jet& jet, const xAOD::VertexContainer* vertices,
+ const float minTrackPt, const std::vector<const xAOD::TrackParticle*>& tracks,
+ const jet::TrackVertexAssociation* tva) const;
+
+ // Local method to calculate NumTrk, SumPtTrk, and TrackWidth for one vertex
+ TrackMomentStruct
+ getTrackMoments(const xAOD::Jet&, const xAOD::Vertex* vertex, const float minTrackPt,
+ const std::vector<const xAOD::TrackParticle*>& tracks,
+ const jet::TrackVertexAssociation* tva) const;
+
+ // Parse the float to get a moment base name
+ const std::string getMomentBaseName(const float minTrackPt) const;
+
+};
+
+#endif
+
diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetVertexFractionTool.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetVertexFractionTool.h
new file mode 100644
index 00000000000..9ef955bdf85
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetVertexFractionTool.h
@@ -0,0 +1,70 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// JetVertexFractionTool.h
+
+#ifndef JETMOMENTTOOLS_JETVERTEXFRACTIONTOOL_H
+#define JETMOMENTTOOLS_JETVERTEXFRACTIONTOOL_H
+
+/// Steven Schramm \n
+/// February 2014
+///
+/// Tool to calculate the jet vertex fraction (JVF)
+/// JVF is a vector<float> with respect to all vertices
+///
+/// Calculation requires three main types of information
+/// 1. Vertex container for the event (from evtStore), which the JVF vector is aligned with
+/// 2. Tracks associated to the input jet (in the jet aux store)
+/// 3. Track vertex association object (from evtStore)
+///
+/// Using this information, the procedure can be broken into three main steps:
+/// 1. Retrieve necessary information
+/// 2. Helper method to create the full vector given the information
+/// 3. Helper method to create the JVF value for one vertex given the information
+
+#include "xAODTracking/Vertex.h"
+#include "xAODTracking/VertexContainer.h"
+#include "xAODTracking/TrackParticle.h"
+#include "xAODTracking/TrackParticleContainer.h"
+
+#include "JetRec/JetModifierBase.h"
+#include "JetEDM/TrackVertexAssociation.h"
+
+#include <vector>
+#include <string>
+
+class JetVertexFractionTool : public JetModifierBase
+{
+ ASG_TOOL_CLASS(JetVertexFractionTool,IJetModifier);
+
+ public:
+ // Constructor from tool name
+ JetVertexFractionTool(const std::string& name);
+
+ // Inherited methods to modify a jet
+ // Calls getJetVertexFraction and puts the result in the jet
+ virtual int modifyJet(xAOD::Jet& jet) const;
+
+ // Local method to calculate and return the JVF vector
+ const std::vector<float> getJetVertexFraction(const xAOD::VertexContainer*, const std::vector<const xAOD::TrackParticle*>&, const jet::TrackVertexAssociation*) const;
+
+ // Local method to calculate the JVF for a given vertex
+ float getJetVertexFraction(const xAOD::Vertex*, const std::vector<const xAOD::TrackParticle*>&, const jet::TrackVertexAssociation*) const;
+
+ // Local method to determine the highest JVF vertex and get an ElementLink to it
+ ElementLink<xAOD::VertexContainer> getMaxJetVertexFraction(const xAOD::VertexContainer*, const std::vector<float>&) const;
+
+ private:
+ // Configurable parameters
+ std::string m_verticesName;
+ std::string m_assocTracksName;
+ std::string m_tvaName;
+
+ std::vector<float> getEmptyJetVertexFraction(const xAOD::VertexContainer*) const;
+
+};
+
+
+#endif
+
diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetWidthTool.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetWidthTool.h
new file mode 100644
index 00000000000..9dc6f4ac9c1
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetWidthTool.h
@@ -0,0 +1,34 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// JetWidthTool.h
+
+#ifndef JetMomentTools_JetWidthTool_H
+#define JetMomentTools_JetWidthTool_H
+
+/// Steven Schramm and David Adams \n
+/// February 2014
+///
+/// Tool to calculate the jet width.
+
+#include "JetRec/JetModifierBase.h"
+
+class JetWidthTool : public JetModifierBase {
+ ASG_TOOL_CLASS(JetWidthTool, IJetModifier);
+
+public:
+
+ // Constructor from tool name.
+ JetWidthTool(std::string myname);
+
+ // Inherited method to modify a jet.
+ // Calls width and puts the result on the jet.
+ virtual int modifyJet(xAOD::Jet& jet) const;
+
+ // Local method to calculate and return the width.
+ double width(const xAOD::Jet& jet) const;
+
+};
+
+#endif
diff --git a/Reconstruction/Jet/JetMomentTools/Root/JetMuonSegmentMomentsTool.cxx b/Reconstruction/Jet/JetMomentTools/Root/JetMuonSegmentMomentsTool.cxx
new file mode 100644
index 00000000000..b781cb906e9
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/Root/JetMuonSegmentMomentsTool.cxx
@@ -0,0 +1,31 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+
+#include "JetMomentTools/JetMuonSegmentMomentsTool.h"
+
+JetMuonSegmentMomentsTool::JetMuonSegmentMomentsTool(const std::string& name)
+ : JetModifierBase(name)
+ , m_assocMuonSegName("GhostMuonSegment")
+{
+ declareProperty("AssociatedMuonSegments", m_assocMuonSegName);
+}
+
+int JetMuonSegmentMomentsTool::modifyJet(xAOD::Jet& jet) const {
+ // Get muon segments associated to the jet
+ // Note that there are often no segments
+ std::vector<const xAOD::IParticle*> segments;
+ bool havesegs = jet.getAssociatedObjects(m_assocMuonSegName, segments);
+
+ if ( ! havesegs ) {
+ ATH_MSG_ERROR("Muon segments not found.");
+ return 1;
+ } else {
+ ATH_MSG_DEBUG("Jet has " << segments.size() << " muon segments.");
+ }
+
+ return 0;
+}
+
+
diff --git a/Reconstruction/Jet/JetMomentTools/Root/JetPtAssociationTool.cxx b/Reconstruction/Jet/JetMomentTools/Root/JetPtAssociationTool.cxx
new file mode 100644
index 00000000000..cea7a27d6f1
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/Root/JetPtAssociationTool.cxx
@@ -0,0 +1,123 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// JetPtAssociationTool.cxx
+
+#include "JetMomentTools/JetPtAssociationTool.h"
+
+using std::string;
+using xAOD::IParticle;
+using xAOD::Jet;
+using xAOD::JetContainer;
+using xAOD::JetConstituent;
+using xAOD::JetConstituentVector;
+
+//**********************************************************************
+
+JetPtAssociationTool::JetPtAssociationTool(std::string myname)
+: JetModifierBase(myname) {
+ declareProperty("AssociationName", m_aname);
+ declareProperty("InputContainer", m_conname);
+}
+
+//**********************************************************************
+
+StatusCode JetPtAssociationTool::initialize() {
+ return StatusCode::SUCCESS;
+}
+
+//**********************************************************************
+
+int JetPtAssociationTool::modifyJet(xAOD::Jet& jet) const {
+ ATH_MSG_DEBUG("Processing jet " << jet.index());
+ APVector apins;
+ // Retrieve the association vector.
+ if ( ! jet.getAssociatedObjects(m_aname, apins) ) {
+ ATH_MSG_WARNING("Jet does not have association vector " << m_aname);
+ return 1;
+ }
+ // Retrieve the container of jets to be matched.
+ const JetContainer* pjets = evtStore()->retrieve<const JetContainer>(m_conname);
+ if ( pjets == 0 ) {
+ ATH_MSG_WARNING("Matching jet container not found: " << m_conname);
+ return 2;
+ }
+ // Match associated particle to jets.
+ FloatVector ptfs;
+ if ( ptfrac(apins, *pjets, ptfs) ) {
+ ATH_MSG_WARNING("Match of association vector to jets failed");
+ return 3;
+ }
+ // Find the best match.
+ float frac = 0.0;
+ unsigned int ijet_matched = ptfs.size();
+ for ( unsigned int ijet=0; ijet<ptfs.size(); ++ijet ) {
+ ATH_MSG_VERBOSE(" Pt fraction[" << ijet << "]: " << ptfs[ijet]);
+ if ( ptfs[ijet] > frac ) {
+ frac = ptfs[ijet];
+ ijet_matched = ijet;
+ }
+ }
+ ElementLink<JetContainer> el;
+ if ( ijet_matched < ptfs.size() ) {
+ el = ElementLink<JetContainer>(*pjets, ijet_matched);
+ }
+ string sfrac = m_aname + "AssociationFraction";
+ string slink = m_aname + "AssociationLink";
+ ATH_MSG_DEBUG(" Setting " << sfrac << " = " << frac);
+ jet.setAttribute(sfrac, frac);
+ ATH_MSG_DEBUG(" Setting " << slink << " = "
+ << el.dataID() << "[" << el.index() << "]");
+ jet.setAttribute(slink, el);
+ return 0;
+}
+
+//**********************************************************************
+
+int JetPtAssociationTool::
+ptfrac(const APVector& apins, const xAOD::JetContainer& jets, FloatVector& ptfs) const {
+ ptfs.clear();
+ ptfs.resize(jets.size(), 0.0);
+ double ptden = 0.0;
+ for ( const IParticle* ppar : apins ) {
+ if ( ppar == 0 ) return -1;
+ ptden += ppar->pt();
+ }
+ if ( ptden < 0.0 ) return -2;
+ if ( ptden < 1.e-20 ) return 0;
+ ATH_MSG_DEBUG("Match jet count: " << jets.size());
+ for ( unsigned int ijet=0; ijet<jets.size(); ++ijet ) {
+ const Jet* pjet = jets[ijet];
+ if ( pjet == 0 ) return -3;
+ APVector apouts;
+ match(apins, *pjet, apouts);
+ double ptsum = 0.0;
+ for ( const IParticle* ppar : apouts ) {
+ if ( ppar == 0 ) return -4;
+ ptsum += ppar->pt();
+ }
+ ptfs[ijet] = ptsum/ptden;
+ }
+ return 0;
+}
+
+//**********************************************************************
+
+int JetPtAssociationTool::
+match(const APVector& aps, const xAOD::Jet& jet, APVector& apvs) const {
+ apvs.clear();
+ const JetConstituentVector cons = jet.getConstituents();
+ for ( const JetConstituent* pcon : cons ) {
+ const IParticle* pparcon = pcon->rawConstituent();
+ for ( const IParticle* pparap : aps ) {
+ if ( pparap == pparcon ) {
+ apvs.push_back(pparcon);
+ break;
+ }
+ }
+ }
+ return 0;
+}
+
+//**********************************************************************
diff --git a/Reconstruction/Jet/JetMomentTools/Root/JetTrackMomentsTool.cxx b/Reconstruction/Jet/JetMomentTools/Root/JetTrackMomentsTool.cxx
new file mode 100644
index 00000000000..a5d7f7840c2
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/Root/JetTrackMomentsTool.cxx
@@ -0,0 +1,161 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+
+#include "JetMomentTools/JetTrackMomentsTool.h"
+#include <sstream>
+#include "JetUtils/JetDistances.h"
+
+JetTrackMomentsTool::JetTrackMomentsTool(const std::string& name)
+ : JetModifierBase(name)
+ , m_vertexContainer("")
+ , m_assocTracksName("")
+ , m_tva("")
+ , m_minTrackPt()
+{
+ declareProperty("VertexContainer",m_vertexContainer);
+ declareProperty("AssociatedTracks",m_assocTracksName);
+ declareProperty("TrackVertexAssociation",m_tva);
+ declareProperty("TrackMinPtCuts",m_minTrackPt);
+}
+
+int JetTrackMomentsTool::modifyJet(xAOD::Jet& jet) const {
+
+ // Get input vertex collection
+ const xAOD::VertexContainer* vertexContainer = nullptr;
+ if ( evtStore()->retrieve(vertexContainer,m_vertexContainer).isFailure()
+ || vertexContainer == nullptr ) {
+ ATH_MSG_ERROR("Could not retrieve the VertexContainer from evtStore: "
+ << m_vertexContainer);
+ return 1;
+ }
+
+ // Get the track-vertex association
+ const jet::TrackVertexAssociation* tva = nullptr;
+ if ( evtStore()->retrieve(tva,m_tva).isFailure() || tva==nullptr ) {
+ ATH_MSG_ERROR("Could not retrieve the TrackVertexAssociation from evtStore: "
+ << m_tva);
+ return 2;
+ }
+
+#if 0
+ // Get the tracks associated to the jet
+ // Note that there may be no tracks
+ const std::vector<const xAOD::IParticle*> tracksAsParticles;
+ bool havetracks = jet.getAssociatedObjects(m_assocTracksName, tracksAsParticles);
+
+ if ( ! havetracks ) {
+ ATH_MSG_WARNING("Associted tracks not found");
+ }
+
+ // Do the dynamic_cast once for the particles instead of repeatedly
+ std::vector<const xAOD::TrackParticle*> tracks;
+ tracks.resize(tracksAsParticles.size());
+ for (size_t iTrack = 0; iTrack < tracksAsParticles.size(); ++iTrack)
+ {
+ const xAOD::TrackParticle* track = dynamic_cast<const xAOD::TrackParticle*>(tracksAsParticles.at(iTrack));
+ if (!track)
+ {
+ ATH_MSG_ERROR("Error converting index " << iTrack << " from IParticle to TrackParticle");
+ return 3;
+ }
+ tracks[iTrack] = track;
+ }
+#endif
+
+ // Retrieve the associated tracks.
+ std::vector<const xAOD::TrackParticle*> tracks;
+ bool havetracks = jet.getAssociatedObjects(m_assocTracksName, tracks);
+ if ( ! havetracks ) ATH_MSG_WARNING("Associated tracks not found");
+ ATH_MSG_DEBUG("Successfully retrieved track particles");
+
+ // For each track cut, get the associated moments
+ for (size_t iCut = 0; iCut < m_minTrackPt.size(); ++iCut) {
+ // Get info
+ const float minPt = m_minTrackPt.at(iCut);
+ const std::vector<TrackMomentStruct> moments =
+ getTrackMoments(jet,vertexContainer,minPt,tracks,tva);
+ const std::string baseName = getMomentBaseName(minPt);
+ // Collate info
+ std::vector<int> numTrkVec; numTrkVec.resize(moments.size());
+ std::vector<float> sumPtTrkVec; sumPtTrkVec.resize(moments.size());
+ std::vector<float> trackWidthVec; trackWidthVec.resize(moments.size());
+ for ( size_t iVertex = 0; iVertex < moments.size(); ++iVertex ) {
+ numTrkVec[iVertex] = moments.at(iVertex).numTrk;
+ sumPtTrkVec[iVertex] = moments.at(iVertex).sumPtTrk;
+ trackWidthVec[iVertex] = moments.at(iVertex).trackWidth;
+ }
+ // Set info
+ jet.setAttribute("NumTrk"+baseName,numTrkVec);
+ jet.setAttribute("SumPtTrk"+baseName,sumPtTrkVec);
+ jet.setAttribute("TrackWidth"+baseName,trackWidthVec);
+ }
+
+ return 0;
+}
+
+
+const std::vector<JetTrackMomentsTool::TrackMomentStruct> JetTrackMomentsTool::getTrackMoments(const xAOD::Jet& jet, const xAOD::VertexContainer* vertices, const float minTrackPt, const std::vector<const xAOD::TrackParticle*>& tracks, const jet::TrackVertexAssociation* tva) const
+{
+ std::vector<TrackMomentStruct> moments;
+ moments.resize(vertices->size());
+
+ for (size_t iVertex = 0; iVertex < vertices->size(); ++iVertex)
+ moments[iVertex] = getTrackMoments(jet,vertices->at(iVertex),minTrackPt,tracks,tva);
+
+ return moments;
+}
+
+JetTrackMomentsTool::TrackMomentStruct JetTrackMomentsTool::getTrackMoments(const xAOD::Jet& jet, const xAOD::Vertex* vertex, const float minTrackPt, const std::vector<const xAOD::TrackParticle*>& tracks, const jet::TrackVertexAssociation* tva) const
+{
+ // Prepare the moments
+ TrackMomentStruct moments;
+ moments.numTrk = 0;
+ moments.sumPtTrk = 0;
+ moments.trackWidth = 0;
+
+ // Prepare const vars for jet eta/phi
+ const float jetEta = jet.eta();
+ const float jetPhi = jet.phi();
+
+ // Loop over the tracks
+ for (size_t iTrack = 0; iTrack < tracks.size(); ++iTrack)
+ {
+ const xAOD::TrackParticle* track = tracks.at(iTrack);
+ const float trackPt = track->pt();
+
+ // Skip the track if it doesn't pass the cut
+ if (trackPt < minTrackPt)
+ continue;
+
+ // Skip the track if it's not from the vertex in question
+ if( vertex == nullptr || vertex != tva->associatedVertex(track) ) continue ;
+
+ // Calculate necessary info for the moments
+ const double deltaR = jet::JetDistances::deltaR(jetEta, jetPhi, track->eta(), track->phi() );
+
+ // Adjust values as necessary for this track
+ moments.numTrk += 1;
+ moments.sumPtTrk += trackPt;
+ moments.trackWidth += deltaR * trackPt;
+ }
+
+ // Finish processing the moments
+ moments.trackWidth = moments.sumPtTrk > 0 ? moments.trackWidth / moments.sumPtTrk : -1;
+
+ return moments;
+}
+
+const std::string JetTrackMomentsTool::getMomentBaseName(const float minTrackPt) const
+{
+ int value = round(minTrackPt);
+ if (fabs(value - minTrackPt) > 0.1)
+ ATH_MSG_WARNING("Cut float and int disagree: " << minTrackPt << " float vs " << value << " int");
+
+ std::ostringstream sout;
+ sout << "Pt" << value;
+ return sout.str();
+}
+
+
diff --git a/Reconstruction/Jet/JetMomentTools/Root/JetVertexFractionTool.cxx b/Reconstruction/Jet/JetMomentTools/Root/JetVertexFractionTool.cxx
new file mode 100644
index 00000000000..72948bd303f
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/Root/JetVertexFractionTool.cxx
@@ -0,0 +1,115 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+
+#include "JetMomentTools/JetVertexFractionTool.h"
+
+JetVertexFractionTool::JetVertexFractionTool(const std::string& name)
+: JetModifierBase(name)
+, m_verticesName("")
+, m_assocTracksName("")
+, m_tvaName("")
+{
+ declareProperty("VertexContainer",m_verticesName);
+ declareProperty("AssociatedTracks",m_assocTracksName);
+ declareProperty("TrackVertexAssociation",m_tvaName);
+}
+
+int JetVertexFractionTool::modifyJet(xAOD::Jet& jet) const {
+
+ // Get the vertices container
+ const xAOD::VertexContainer* vertices = NULL;
+ if ( evtStore()->retrieve(vertices,m_verticesName).isFailure() ) {
+ ATH_MSG_ERROR("Could not retrieve the VertexContainer from evtStore: " << m_verticesName);
+ return 1;
+ }
+ ATH_MSG_DEBUG("Successfully retrieved VertexContainer from evtStore: " << m_verticesName);
+
+ // Get the tracks associated to the jet
+ // Note that there may be no tracks - this is both normal and an error case
+ // In this case, just fill a vector with zero and don't set the highest vtx moment
+ std::vector<const xAOD::TrackParticle*> tracks;
+ if ( ! jet.getAssociatedObjects(m_assocTracksName, tracks) ) {
+ ATH_MSG_WARNING("Associated tracks not found.");
+ }
+
+ // Get the TVA object
+ const jet::TrackVertexAssociation* tva = NULL;
+ if (evtStore()->retrieve(tva,m_tvaName).isFailure()) {
+ ATH_MSG_ERROR("Could not retrieve the TrackVertexAssociation from evtStore: " << m_tvaName);
+ return 3;
+ }
+ ATH_MSG_DEBUG("Successfully retrieved TrackVertexAssociation from evtStore: " << m_tvaName);
+
+ // Get and set the JVF vector
+ const std::vector<float> jvf = getJetVertexFraction(vertices,tracks,tva);
+ jet.setAttribute("JVF",jvf);
+
+ // Get and set the highest JVF vertex
+ jet.setAttribute("HighestJVFVtx",getMaxJetVertexFraction(vertices,jvf));
+
+ // Done
+ return 0;
+}
+
+
+
+std::vector<float> JetVertexFractionTool::getEmptyJetVertexFraction(const xAOD::VertexContainer* vertices) const
+{
+ std::vector<float> jvf;
+ jvf.resize(vertices->size());
+ for (size_t iVertex = 0; iVertex < vertices->size(); ++iVertex)
+ jvf.at(iVertex) = 0;
+
+ return jvf;
+}
+
+const std::vector<float> JetVertexFractionTool::getJetVertexFraction(const xAOD::VertexContainer* vertices, const std::vector<const xAOD::TrackParticle*>& tracks, const jet::TrackVertexAssociation* tva) const
+{
+ std::vector<float> jvf = getEmptyJetVertexFraction(vertices);
+
+ for (size_t iVertex = 0; iVertex < vertices->size(); ++iVertex)
+ jvf.at(iVertex) = getJetVertexFraction(vertices->at(iVertex),tracks,tva);
+
+ return jvf;
+}
+
+
+float JetVertexFractionTool::getJetVertexFraction(const xAOD::Vertex* vertex, const std::vector<const xAOD::TrackParticle*>& tracks, const jet::TrackVertexAssociation* tva) const
+{
+ float sumTrackPV = 0;
+ float sumTrackAll = 0;
+
+ for (size_t iTrack = 0; iTrack < tracks.size(); ++iTrack) {
+ const xAOD::TrackParticle* track = tracks.at(iTrack);
+ const xAOD::Vertex* ptvtx = tva->associatedVertex(track);
+ if ( ptvtx != nullptr ) {
+ if ( ptvtx->index() == vertex->index() ) sumTrackPV += track->pt();
+ }
+ sumTrackAll += track->pt();
+ }
+
+ return sumTrackAll>0 ? sumTrackPV/sumTrackAll : -1;
+}
+
+
+ElementLink<xAOD::VertexContainer> JetVertexFractionTool::getMaxJetVertexFraction(const xAOD::VertexContainer* vertices, const std::vector<float>& jvf) const
+{
+ size_t maxIndex = 0;
+ float maxVal = -100;
+
+ for (size_t iVertex = 0; iVertex < jvf.size(); ++iVertex)
+ if (jvf.at(iVertex) > maxVal)
+ {
+ maxIndex = iVertex;
+ maxVal = jvf.at(iVertex);
+ }
+
+ ElementLink<xAOD::VertexContainer> link = ElementLink<xAOD::VertexContainer>(*vertices,vertices->at(maxIndex)->index());
+
+ return link;
+}
+
+
+
diff --git a/Reconstruction/Jet/JetMomentTools/Root/JetWidthTool.cxx b/Reconstruction/Jet/JetMomentTools/Root/JetWidthTool.cxx
new file mode 100644
index 00000000000..a006afce73e
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/Root/JetWidthTool.cxx
@@ -0,0 +1,53 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// JetWidthTool.cxx
+
+#include "JetMomentTools/JetWidthTool.h"
+
+#include "xAODJet/JetConstituentVector.h"
+#include "JetUtils/JetDistances.h"
+
+//**********************************************************************
+
+JetWidthTool::JetWidthTool(std::string myname)
+: JetModifierBase(myname) { }
+
+//**********************************************************************
+
+int JetWidthTool::modifyJet(xAOD::Jet& jet) const {
+ jet.setAttribute("Width", width(jet));
+ return 0;
+}
+
+//**********************************************************************
+
+double JetWidthTool::width(const xAOD::Jet& jet) const {
+
+ // Get the constituents of the jet
+ const xAOD::JetConstituentVector constituents = jet.getConstituents();
+ xAOD::JetConstituentVector::iterator iter = constituents.begin();
+ xAOD::JetConstituentVector::iterator iEnd = constituents.end();
+
+ // TODO: Switch to using helper function once JetUtils has been updated
+ // Set the width
+ // jetWidth = JetKinematics::ptWeightedWidth(iter,iEnd,&jet);
+
+ // Calculate the pt weighted width
+ const double jetEta = jet.eta();
+ const double jetPhi = jet.phi();
+ double weightedWidth = 0;
+ double ptSum = 0;
+
+ for ( ; iter != iEnd; ++iter) {
+ const double dR = jet::JetDistances::deltaR(jetEta, jetPhi, iter->eta(), iter->phi() );
+ const double pt = iter->pt();
+ weightedWidth += dR * pt;
+ ptSum += pt;
+ }
+
+ return ptSum > 0 ? weightedWidth/ptSum : -1;
+}
+
+//**********************************************************************
diff --git a/Reconstruction/Jet/JetMomentTools/cmt/Makefile.RootCore b/Reconstruction/Jet/JetMomentTools/cmt/Makefile.RootCore
new file mode 100644
index 00000000000..08a084251e3
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/cmt/Makefile.RootCore
@@ -0,0 +1,24 @@
+# this makefile also gets parsed by shell scripts
+# therefore it does not support full make syntax and features
+# edit with care
+
+# for full documentation check:
+# https://twiki.cern.ch/twiki/bin/viewauth/Atlas/RootCore#Package_Makefile
+
+PACKAGE = JetMomentTools
+PACKAGE_PRELOAD =
+PACKAGE_CXXFLAGS =
+PACKAGE_OBJFLAGS =
+PACKAGE_LDFLAGS =
+PACKAGE_BINFLAGS =
+PACKAGE_LIBFLAGS =
+PACKAGE_DEP = xAODJet JetEDM JetInterface JetRec
+PACKAGE_TRYDEP =
+PACKAGE_CLEAN =
+PACKAGE_NOGRID =
+PACKAGE_PEDANTIC = 0
+PACKAGE_NOOPT = 0
+PACKAGE_NOCC = 0
+PACKAGE_REFLEX = 0
+
+include $(ROOTCOREDIR)/Makefile-common
diff --git a/Reconstruction/Jet/JetMomentTools/cmt/requirements b/Reconstruction/Jet/JetMomentTools/cmt/requirements
new file mode 100755
index 00000000000..031621bd8ed
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/cmt/requirements
@@ -0,0 +1,31 @@
+package JetMomentTools
+
+author P-A Delsart <delsart at in2p3 fr>
+
+private
+use xAODCaloEvent xAODCaloEvent-* Event/xAOD
+use CaloDetDescr CaloDetDescr-* Calorimeter
+use CaloEvent CaloEvent-* Calorimeter
+use PathResolver PathResolver-* Tools
+use AthenaKernel AthenaKernel-* Control
+
+public
+use AtlasPolicy AtlasPolicy-*
+use AtlasROOT AtlasROOT-* External
+use GaudiInterface GaudiInterface-* External
+use CaloIdentifier CaloIdentifier-* Calorimeter
+use xAODJet xAODJet-* Event/xAOD
+use JetEDM JetEDM-* Reconstruction/Jet
+use JetRec JetRec-* Reconstruction/Jet
+use JetUtils JetUtils-* Reconstruction/Jet
+use JetRecTools JetRecTools-* Reconstruction/Jet
+use xAODTracking xAODTracking-* Event/xAOD
+use AsgTools AsgTools-* Control/AthToolSupport
+use LArElecCalib LArElecCalib-* LArCalorimeter
+
+library JetMomentTools *.cxx ../Root/*.cxx -s=components *.cxx
+apply_pattern component_library
+
+apply_pattern declare_runtime files="*.root *.txt"
+
+#macro_prepend asNeeded_linkopt "-lJetEDM "
diff --git a/Reconstruction/Jet/JetMomentTools/python/GhostAssociation.py b/Reconstruction/Jet/JetMomentTools/python/GhostAssociation.py
new file mode 100644
index 00000000000..645c1e94ee9
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/python/GhostAssociation.py
@@ -0,0 +1,148 @@
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+from AthenaCommon.AlgSequence import AlgSequence
+def interpretJetName(jetcollName, finder = None,input=None, mainParam=None):
+ # first step : guess the finder, input , mainParam, if needed
+ if finder is None:
+ for a in [ 'AntiKt','CamKt','Kt', 'Cone','SISCone','CMSCone']:
+ if jetcollName.startswith(a):
+ finder = a
+ break
+ if finder is None:
+ print "Error could not guess jet finder type in ",jetcollName
+ return
+
+ if input is None:
+ for i in ['LCTopo','Tower','Topo', "Truth"]:
+ if i in jetcollName:
+ input = i
+ if i== "Tower":
+ if 'Ghost' in jetcollName:
+ input = 'Tower'
+ else:
+ input = "TopoTower"
+ break
+ if input is None:
+ print "Error could not guess input type in ",jetcollName
+ return
+
+ if mainParam is None:
+ # get the 2 chars following finder :
+ mp = jetcollName[len(finder):len(finder)+2]
+ mp = mp[0] if not mp[1] in '0123456789' else mp
+ try :
+ mainParam = float(mp)/10.
+ except ValueError :
+ print "Error could not guess main parameter in ",jetcollName
+ return
+
+ return finder, input, mainParam
+
+
+
+
+def setupGhostAssociator(assocName, assocType, assocKeys="", seq=AlgSequence()):
+ from JetMomentTools.JetMomentToolsConf import GhostAssociator_JetINav4MomAssociation_ , GhostTrackAssociator , GhostAssociator_TruthPtAssociation_
+
+ # make sure the associated container key is given as a list
+ assocKeys = [assocKeys] if not isinstance(assocKeys, list) else assocKeys
+
+ # get the correct associator
+ if assocType == "TrackAssociation":
+ assocTool = GhostTrackAssociator(assocName, AssociationName = assocName, AssociatedContainers = assocKeys )
+ assocTool.UseTrackSelector = True
+ from JetRec.TrackSelectionForJets import getDefaultJetVtxTrackHelper
+ assocTool.TrackSelector = getDefaultJetVtxTrackHelper()
+ elif assocType == "TruthAssociation":
+ from JetSimTools.JetTruthParticleFilterGetter import JetTruthParticleFilterGetter
+ if assocKeys == [""]:
+ assocKeys = [JetTruthParticleFilterGetter(seq=seq)._outputName]
+ assocTool = GhostAssociator_JetINav4MomAssociation_(assocName, AssociationName = assocName, AssociatedContainers = assocKeys, AssociatedPtMoment="pt_truth" )
+ elif assocType == "TruthPtAssociation":
+ from JetSimTools.JetTruthParticleFilterGetter import JetTruthParticleFilterGetter
+ if assocKeys == [""]:
+ assocKeys = [JetTruthParticleFilterGetter(seq=seq)._outputName]
+ assocTool = GhostAssociator_TruthPtAssociation_(assocName, AssociationName = assocName, AssociatedContainers = assocKeys )
+ else:
+ #print "Can't associate to ",assocType, " will use a generic association JetINav4MomAssociation"
+ assocTool = GhostAssociator_JetINav4MomAssociation_(assocName, AssociationName = assocName, AssociatedContainers = assocKeys )
+
+ return assocTool
+
+
+def defaultGhostAssociator( assocName, seq=AlgSequence() ):
+ """Shorcuts for creating usual ghost associators """
+ if assocName == "TrackAssoc":
+ from JetRec.TrackSelectionForJets import trackParticleContainer
+ return setupGhostAssociator( "TrackAssoc", "TrackAssociation", [trackParticleContainer] , seq=seq)
+ elif assocName == "TruthPt":
+ return setupGhostAssociator( "TruthPt", "TruthPtAssociation", seq=seq)
+ elif assocName == "TruthAssoc":
+ return setupGhostAssociator( "TruthAssoc", "TruthAssociation", seq=seq)
+ elif assocName == "TruthBHadron":
+ return setupGhostAssociator( "TruthBHadron", "", ["TruthBHadrons"], seq=seq)
+
+
+
+def setupGhostAssociationTool(jetcollName, ghostsAssoc = [], finder = None,input=None, mainParam=None, seq = AlgSequence()):
+ from JetMomentTools.JetMomentToolsConf import GhostAssociationTool
+ from JetRec.JetGetters import getStandardInputTools, getStandardFinderTools
+
+ # first step : guess the finder, input , mainParam, if needed
+ finder, input, mainParam = interpretJetName( jetcollName, finder, input, mainParam)
+ replaceNegEnergy = 'Ghost' in input
+
+ #print 'Will create an association by re-running with this jet alg : ', finder, mainParam, input
+
+ # Create the tool
+ t= GhostAssociationTool("JetGhostAsso")
+ ghostAssocTools = []
+ for g in ghostsAssoc:
+ if isinstance(g,str):
+ ghostAssocTools.append( defaultGhostAssociator(g,seq=seq) )
+ else: ## Assume g is a GhostAssociator
+ ghostAssocTools.append( g)
+
+ t.GhostAssociators = ghostAssocTools
+ t.JetFindingSequence = getStandardInputTools( input, replaceNegEnergy =replaceNegEnergy )+getStandardFinderTools( finder , mainParam)
+ t.JetFindingSequence[0].NoCleanup = True # Necessary !
+
+ return t
+
+
+def addGhostAssociation(jetcollName, ghostAssoc = [], finder = None,input=None, mainParam=None, seq=AlgSequence()):
+ from JetRec.JetMomentGetter import make_JetMomentGetter
+
+ t = setupGhostAssociationTool( jetcollName, ghostAssoc, finder ,input, mainParam, seq = seq)
+ g=make_JetMomentGetter(jetcollName,[t],seq=seq)
+ return g.getJetMomentCalc()
+
+
+def setupGhostAssociationTool2(jetcollName, assocName, assocType, assocKeys="", finder = None,input=None, mainParam=None):
+ from JetMomentTools.JetMomentToolsConf import GhostAssociationTool
+ from JetRec.JetGetters import getStandardInputTools, getStandardFinderTools
+
+ # first step : guess the finder, input , mainParam, if needed
+ finder, input, mainParam = interpretJetName( jetcollName, finder, input, mainParam)
+
+ #print 'Will create an association by re-running with this jet alg : ', finder, mainParam, input
+
+ # Create the tool
+ t= GhostAssociationTool("JetAsso"+assocName)
+ t.GhostAssociators = [ setupGhostAssociator(assocName , assocType, assocKeys) ]
+ # t.GhostAssociators = [ associatorClass(assocName, AssociationName = assocName, AssociatedContainers = assocKeys ) ]
+ t.JetFindingSequence = getStandardInputTools( input )+getStandardFinderTools( finder , mainParam)
+ t.JetFindingSequence[0].NoCleanup = True # Necessary !
+
+ return t
+
+
+
+
+
+def addGhostAssociation2(jetcollName, assocName, assocType, assocKeys="", finder = None,input=None, mainParam=None):
+ from JetRec.JetMomentGetter import make_JetMomentGetter
+
+ t = setupGhostAssociationTool( jetcollName, assocName, assocType, assocKeys, finder ,input, mainParam)
+ g=make_JetMomentGetter(jetcollName,[t])
+ return g.getJetMomentCalc()
diff --git a/Reconstruction/Jet/JetMomentTools/python/JetCopyMomentsAlg.py b/Reconstruction/Jet/JetMomentTools/python/JetCopyMomentsAlg.py
new file mode 100644
index 00000000000..34de1861786
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/python/JetCopyMomentsAlg.py
@@ -0,0 +1,72 @@
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+# @file: JetMomentTools/python/JetCopyMomentsAlg.py.py
+# @purpose: <Copy some moments from a collection to another >
+# @author: P-A Delsart <delsart in2p3 fr>
+
+__doc__ = 'some documentation here'
+__version__ = '$Revision: 1.5 $'
+__author__ = 'Sebastien Binet <binet@cern.ch>'
+
+import AthenaCommon.SystemOfUnits as Units
+import AthenaPython.PyAthena as PyAthena
+from AthenaPython.PyAthena import StatusCode
+import ROOT
+
+class JetCopyMomentsAlg (PyAthena.Alg):
+ 'put some documentation here'
+ def __init__(self, name='JetCopyMomentsAlg', **kw):
+ ## init base class
+ kw['name'] = name
+ super(JetCopyMomentsAlg, self).__init__(**kw)
+
+ ## properties and data members
+ self.Moments = kw.get('Moments', []) # default value
+ self.JetCollectionSrc = kw.get('JetCollectionSrc', "") # default value
+ self.JetCollectionDest = kw.get('JetCollectionDest', "") # default value
+ return
+
+ def initialize(self):
+ self.msg.info('==> initialize...')
+ self.sg = PyAthena.py_svc ('StoreGateSvc')
+ return StatusCode.Success
+
+ def execute(self):
+ jcoll_src = list(self.sg.retrieve('JetCollection', self.JetCollectionSrc))
+ jcoll_dest = list(self.sg.retrieve('JetCollection', self.JetCollectionDest))
+
+ # sort on constit scale to help with ordering.
+ jcoll_dest.sort(key= lambda x: x.pt(2) )
+ jcoll_src.sort(key= lambda x: x.pt(2) )
+ getMoment = ROOT.Jet.getMoment
+ moments = self.Moments
+
+ dr = ROOT.JetDistances.deltaR
+ needUnorderedMatch = False
+ for i, (jS, jD) in enumerate(zip(jcoll_src, jcoll_dest)):
+ if dr(jS,jD)>0.05:
+ needUnorderedMatch = True
+ break
+ for m in moments:
+ jD.setMoment(m,jS.getMoment(m))
+
+ if needUnorderedMatch :
+ # this may happen because of final Pt cut and a calibration changing order w.r.t lcscale.
+ # fall back to a double loop with deltaR check on remaining jets
+ jcoll_src = jcoll_src[i:]
+ jcoll_dest = jcoll_dest[i:]
+ for jS in jcoll_src:
+ for jD in jcoll_dest:
+ if dr(jS,jD)<0.05:
+ for m in moments:
+ jD.setMoment(m,jS.getMoment(m))
+
+
+
+ return StatusCode.Success
+
+ def finalize(self):
+ self.msg.info('==> finalize...')
+ return StatusCode.Success
+
+ # class JetCopyMomentsAlg
diff --git a/Reconstruction/Jet/JetMomentTools/python/JetMomentsConfigHelpers.py b/Reconstruction/Jet/JetMomentTools/python/JetMomentsConfigHelpers.py
new file mode 100644
index 00000000000..616f98b1f95
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/python/JetMomentsConfigHelpers.py
@@ -0,0 +1,363 @@
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+"""
+A set of helper function to configure calculations of additionnal jet moments.
+
+"""
+from JetRec.JetAlgConfiguration import checkAndUpdateOptions, standardJetOption
+from RecExConfig.RecFlags import rec
+from AthenaCommon.SystemOfUnits import GeV
+from AthenaCommon.AlgSequence import AlgSequence
+from AthenaCommon.Logging import logging
+from AthenaCommon.JobProperties import jobproperties
+
+from JetRec.JetRecFlags import jetFlags
+from JetRec.JetFinderConfig import setupAreaCalculation
+
+from JetMomentTools.SetupJetMomentTools import *
+
+
+## A dedicated logger
+_momlog = logging.getLogger( 'JetMomentConfig' )
+
+
+# We try to schedule cpu intensive moment calculation only once. This requires to keep track of
+# moments already scheduled for a given collection. Use this dictionnary.
+_specialMomentAlgDict = {}
+
+# Some moments request can translate in modification of option to already scheduled moment tool.
+# so we keep only one main algorithm (JetAlgorithm or JetMomentCalculator) per jet collection, in order to be able to
+# retrieve & modify the relevant moment tool.
+# However this can cause dependency issue :
+# for example, truth ghostAssoc requires JetTruthParticleFilter to run before. If a 1st call doesn't require truth ghost assoc, then later a 2nd call does, then JetTruthParticleFilter will be scheduled *after* the moment calculator alg...
+# workaround : move by jand JetTruthParticleFilter at a proper place in the topSequence or the PreD3PDSequence.
+def moveAlgInSequence(alg, seq, newpos):
+ l = seq.getChildren()
+ delattr(seq, alg.getName())
+ seq.insert(newpos, alg)
+
+
+def checkRequiredAlgPosition(seq, momAlg, requiredAlg):
+ topSeq = AlgSequence()
+
+ if requiredAlg in topSeq:
+ truthInd = topSeq.getChildren().index(requiredAlg)
+ if momAlg in topSeq:
+ index = topSeq.getChildren().index(momAlg)
+ else:
+ # then it is in seq, and seq must be in topSeq
+ index = topSeq.getChildren().index(seq)
+ # make sure requiredAlg is before index
+ if truthInd > index:
+ moveAlgInSequence(requiredAlg, topSeq, index)
+ else:
+ # will move requiredAlg in topSeq before seq and momAlg
+ # where is momAlg ?
+ seqInd = topSeq.getChildren().index(seq)
+ if momAlg in topSeq:
+ index = topSeq.getChildren().index(momAlg)
+ else:
+ index = seqInd
+ delattr(seq, requiredAlg.getName())
+ topSeq.insert(index, requiredAlg)
+
+
+
+def checkTruthParticleFilterPosition(seq, momAlg):
+ from JetSimTools.JetTruthParticleFilterGetter import JetTruthParticleFilterGetter
+ truthAlg = JetTruthParticleFilterGetter().alg
+ checkRequiredAlgPosition(seq, momAlg, truthAlg)
+
+
+
+def retrieveFinderTool(toolSeq):
+ from JetRec.JetRecConf import JetConeFinderTool, JetFastJetFinderTool
+ for t in toolSeq:
+ if isinstance(t, JetFastJetFinderTool) or isinstance(t, JetConeFinderTool):
+ return t
+
+
+import inspect
+def whosdaddy():
+ l = [ ll[3] for ll in inspect.stack()[2:] ]
+ return l
+
+
+
+
+
+# Helper functions to schedule jet moments.
+# They are typically wrappers around functions in SetupJetMomentTools which are as generic as possible
+# The wrappers below also perform tasks related to the particular jet collection or configuration sequence.
+def scheduleMFTool(toolType,jetcollname, jf, alg, seq ,*l):
+ from JetMomentTools.SetupJetMomentTools import getMatchingFractionTool
+ tool, requiredAlg = getMatchingFractionTool(jetcollname, toolType, seq)
+ checkRequiredAlgPosition(seq, alg, requiredAlg)
+ return tool
+
+from JetMomentTools.JetMomentToolsConf import JetWidthTool, NumTowerInJetTool, JetEinSamplingTool
+
+
+# The dictionnary below lists the known standard moments and how to compute them.
+# keys in the dict are shorcuts, they might represent several moments
+# values are function used to actually schedule the moment calculations
+defaultMoments = {
+ 'ghostAssoc' : lambda *l : None,
+ 'area' : lambda jcoll,jetfinder,*l : setupAreaCalculation(jetfinder, "ActiveArea", addArea4Vec=True, addSecondArea=jetFlags.jetPerformanceJob() ) ,
+ 'quality' : lambda jetcollname,*l : getCaloQualityTool(doCellBasedVars = (jetcollname=="AntiKt4TopoEMJets" or "Tower" in jetcollname ), computeFromCluster = "Tower" not in jetcollname ),
+ 'ktDeltaR' : lambda *l : getsubJetsMomentTool(),
+ 'isolation' : lambda *l : getJetIsolationTool("JetIsolation"),
+ 'jvf' : lambda jetcollname,*l : getJetVertexAssociationTool(toolName=jetcollname+'JVAtool') ,
+ 'trackMoments' : lambda *l : getJetTracksMomentTool(),
+ 'truthMF' : lambda *l : scheduleMFTool("Truth",*l) ,
+ 'trackMF' : lambda *l : scheduleMFTool("Track",*l) ,
+ 'larHV' : lambda *l : getJetLArHVMomentTool(),
+ 'width' : lambda *l : JetWidthTool(),
+ 'eInSampling' : lambda *l : JetEinSamplingTool(),
+ 'numTowerJet' : lambda *l : NumTowerInJetTool(),
+ 'badChanCorr' : lambda jc,jf,alg,seq,jetAlgConfigDict : getJetBadChanCorrTool(**jetAlgConfigDict),
+ 'origin' : lambda jetcollname,*l : getORIGINMoments(jetcollname),
+ }
+
+
+## The list of standard moments keys used in reconstruction.
+standardMoments = standardJetOption("jetMoments")
+
+class TestWithMsg(object):
+ """Simply a boolean and a str message explaining why the boolean is False
+ Also supports & and | operation (logical and or) with concatenation of message"""
+ jetcollname=""
+ def __init__(self, test, failmsg):
+ self.test, self.failmsg= test, failmsg
+ def check(self, mName):
+ if self.test:
+ return True
+ _momlog.info( "No "+mName+" moments for "+self.jetcollname+". Reason : "+self.failmsg)
+ return False
+ def __and__(self, o ):
+ test = self.test and o.test
+ failmsg = [ t.failmsg for t in (self,o) if not t.test ]
+ return TestWithMsg( test, ' AND '.join(failmsg) )
+ def __or__(self, o ):
+ test = self.test or o.test
+ failmsg = '' if test else self.failmsg+' AND '+o.failmsg
+ return TestWithMsg( test, failmsg )
+
+
+def checkMomentAvailability(jetcollname, momentList, jetAlgConfigDict={}):
+ """ Input : a jet collection name (jetcollname), list of moment keys (momentList), a jet alg configuration dict (as in JetRec.JetAlgConfiguration)
+ returns a filtered list of moment keys where moments incompatible with the jet collection or the job conditions are removed.
+ """
+ # we'll need to know what the input file is
+ from JetRec.JetGetters import guessInputFileType
+ inputFileType = guessInputFileType()
+
+
+
+ from AthenaCommon.DetFlags import DetFlags
+ from AthenaCommon.AppMgr import ToolSvc
+ from AthenaCommon.GlobalFlags import globalflags
+ from JetRec.TrackSelectionForJets import tracksAvailableForJets
+
+ TestWithMsg.jetcollname = jetcollname
+ # short cut flags
+ # use TestWithMsg rather than simply bool only to be able to display a reason for the rejection of the moment.
+ accessCells = TestWithMsg( (inputFileType=="RDO") or (inputFileType=="ESD") , "No CaloCell available")
+ caloJet = TestWithMsg( ("Topo" in jetcollname) or ("Tower" in jetcollname), "Not calo jets")
+ topoJet = TestWithMsg( "Topo" in jetcollname, "Not TopoCluster jets")
+ lctopoJet = TestWithMsg( "LCTopo" in jetcollname, "Not LC TopoCluster jets")
+ notTruthJet = TestWithMsg( 'Truth' not in jetcollname, "Is Truth jets")
+ tracksAvailable = TestWithMsg( tracksAvailableForJets(), "Tracks selection for jets impossible")
+
+ trueTest = TestWithMsg(True,"")
+
+ # build the conditions at which moments can be built
+ compatibilityDict = {
+ 'area' : TestWithMsg('Cone' not in jetcollname, "no area for cone jets"),
+ 'ghostAssoc' : notTruthJet,
+ 'width' : trueTest,
+ 'numTowerJet' : topoJet & accessCells,
+ 'quality' : caloJet,
+ 'badChanCorr' : caloJet & accessCells,
+ 'constitCalib': lctopoJet & TestWithMsg( not jetAlgConfigDict.get('calibName','').startswith("LC:CONST") , "Calibration differ from LC:CONST") ,
+ 'jvf' : caloJet & tracksAvailable,
+ 'larHV' : caloJet & TestWithMsg( (inputFileType=='RDO') and DetFlags.dcs.LAr_on() and globalflags.DataSource()=='data', "No proper LAr info available"),
+ 'eInSampling' : (caloJet & accessCells) | (topoJet),
+ 'truthMF' : TestWithMsg( rec.doTruth() , "Truth not available") & notTruthJet,
+ 'trackMF' : notTruthJet & tracksAvailable,
+ 'trackMoments': notTruthJet & tracksAvailable,
+
+ }
+
+ # filter moments :
+ availableMoments = set([ m for m in momentList if compatibilityDict.get(m, trueTest).check(m) ])
+ return availableMoments
+
+
+
+
+def findMainAlgForJetCollection(jetcollname, seq = AlgSequence()):
+ """Given a JetCollection name (jetcollname) in input, returns an algorithm computing moments for this collection.
+ The functions explore all know possibilities for finding such algs.
+
+ The return value is (alg, set) where
+ - alg is a JetAlgorithm or a JetMomentCalculator. alg can be None if nothing found or unusable in the case a JetGetter was called with the disable option.
+ - set : is the set of known moments associated to the alg
+ """
+ # try to find a moment calculator for this collection :
+ alg , existing_moments = _specialMomentAlgDict.get( jetcollname, (None,set()) )
+ unknownAlg = alg is None
+
+ if alg is None:
+ # try to find a getter :
+ from JetRec.JetGetters import retrieveConfiguredJetAlg
+ alg = retrieveConfiguredJetAlg(jetcollname, retrieveGetter=True)
+ if alg is not None:
+ # check if usable :
+ if not alg.usable():
+ return ('unusable', set())
+
+ # get the real JetAlgorithm
+ alg = alg.jetAlgorithmHandle()
+ # in this case, no moment have been associated yet
+ existing_moments = set()
+ if alg is None:
+ # check if a collection exists in input file
+ from RecExConfig.ObjKeyStore import objKeyStore
+ if 'JetCollection#'+jetcollname in objKeyStore['inputFile'].list("JetCollection"):
+ # then create a new moment alg for this collection
+ from JetRec.JetMomentGetter import getJetMomentAlgorithm
+ alg = getJetMomentAlgorithm(jetcollname, seq = seq)
+ else:
+ # likely to fail...
+ # try to invoke a standard JetAlgorithm
+ from JetRec.JetRec_defaults import standardJetAlgorithm
+ try :
+ alg = standardJetAlgorithm(jetcollname, seq = seq)
+ except:
+ alg = None #
+ if unknownAlg and alg is not None:
+ # we created a new alg, remember it
+ _specialMomentAlgDict[jetcollname] = (alg, set() )
+ return alg, existing_moments
+
+from JetRec.JetRecConf import JetAlgorithm
+from JetMomentTools.JetMomentToolsConf import JetMomentCalculator
+def addStandardMoments(jetcollname, moments=standardMoments, seq = AlgSequence(), jetAlgConfigDict={}):
+ """Add a list of jet moment calculation for a given jet collection
+ - jetcollname : string, the name of the jetcollection
+ - moments : list of string or function. string are expected to be some of the defaultMoment shortcuts. functions must return a configured moment tool
+ - seq : the AlgSequence to which the calculation is attached
+ (- jetAlgConfigDict : dictionnary, used to pass JetAlgorithm configuration options when this addStandardMoments is called while configuring the JetAlgorithm.)
+ Returns nothing.
+ Effect : retrieve (or invoke if needed) the moment calculation algorithm dedicated to jetcollname (this can be a JetAlgorithm or a JetMomentCalculator),
+ then add to this algorithm the moments requested in defaultMoment (only if they are not already scheduled).
+ """
+
+ if moments == []:
+ # nothing to do. Return now to avoid any inconvenience
+ return
+
+ # retrieve the main algorithm for this jet collection
+ alg, existing_moments = findMainAlgForJetCollection(jetcollname, seq)
+
+
+ # check if the alg is usable ---------
+ if alg == 'unusable':
+ _momlog.info("No usable algorithm for "+jetcollname+". No jet moments added")
+ return
+ if alg is None :
+ # didn't find any alg. We can't add moments
+ _momlog.error("Can't create momens for "+jetcollname+" : no such collection in input or no algorithms scheduled for it")
+ _momlog.error(" ---> possible solution, invoke make_StandardJetGetter() before to schedule a JetAlgorithm for "+jetcollname)
+ return
+
+ isJetMomentCalc = isinstance(alg, JetMomentCalculator)
+
+ # filter moment list according to jet collection and other conditions
+ missing_moments = checkMomentAvailability( jetcollname, set(moments) - existing_moments , jetAlgConfigDict )
+
+ # enforce ghostAssoc, thus jet refinding, if area requested (could separate the 2 options with more config) when using a JetMomentCalculator
+ if 'area' in missing_moments and ('ghostAssoc' not in existing_moments) and isJetMomentCalc:
+ missing_moments.add('ghostAssoc')
+
+
+ jetfinder = None
+ # ghost association is very particular since it needs jet finding.
+ if 'ghostAssoc' in missing_moments:
+ gAssoc = []
+ from JetRec.TrackSelectionForJets import tracksAvailableForJets
+ if tracksAvailableForJets() :
+ gAssoc += ["TrackAssoc"]
+ if rec.doTruth():
+ gAssoc+=[ "TruthAssoc"]
+ if isJetMomentCalc:
+ from JetMomentTools.GhostAssociation import setupGhostAssociationTool
+ gtool = setupGhostAssociationTool(jetcollname, gAssoc, seq=seq)
+ tmpL = alg.CalculatorTools
+ alg.CalculatorTools = [gtool] + tmpL
+ else: # it's a JetAlgorithm
+ tmpL = list(alg.AlgTools)
+ from JetRec.JetGetters import adaptToolListToGhostAssociation
+ alg.AlgTools = adaptToolListToGhostAssociation(tmpL, gAssoc, seq)
+ if rec.doTruth():
+ checkTruthParticleFilterPosition( seq, alg)
+
+ jetfinder = alg.CalculatorTools['JetGhostAsso'].JetFindingSequence["JetFastJetFinder"] if isJetMomentCalc \
+ else alg.AlgTools[0].JetFindingSequence["JetFastJetFinder"]
+
+ # we are done with ghostAssoc
+ missing_moments.remove('ghostAssoc')
+ existing_moments.add('ghostAssoc')
+ else: # still need to retrieve the jetfinder
+ if 'ghostAssoc' in existing_moments:
+ jetfinder = alg.CalculatorTools['JetGhostAsso'].JetFindingSequence["JetFastJetFinder"] if isJetMomentCalc \
+ else alg.AlgTools[0].JetFindingSequence["JetFastJetFinder"]
+ else:
+ jetfinder = retrieveFinderTool( alg.CalculatorTools if isJetMomentCalc else alg.AlgTools)
+
+ momentTools = []
+ # Loop over requested moment types and set-up the related tools
+ for momentType in missing_moments:
+ if callable(momentType):
+ # then assume a function returning a tool has been supplied
+ func = momentType
+ else:
+ func = defaultMoments[momentType]
+ # call the setup function, add the returned tool to the sequence
+ tool = func( jetcollname, jetfinder, alg, seq, jetAlgConfigDict)
+ if tool is not None: # some function just adapts existing tools
+ momentTools.append(tool)
+
+ # add the moment tools list to the relevant alg :
+ if isJetMomentCalc:
+ alg.CalculatorTools += momentTools
+ else:
+ alg.AlgTools += momentTools
+
+ _specialMomentAlgDict[jetcollname] = (alg, missing_moments.union(existing_moments) )
+ return alg
+
+
+
+
+def recommendedAreaAndJVFMoments(jetcollname, oldMomentNames=False,seq=AlgSequence()):
+ """Recompute the area, JVF and track related jet moment on jetcollname.
+
+ jetcollname : string a JetCollection StoreGate key
+ oldMomentNames : bool, if True the moment naming will be the same as in 17.2.5 and previous releases.
+ """
+
+ alg = addStandardMoments(jetcollname,moments = ['ghostAssoc','area','jvf', 'trackMoments'], seq=seq)
+
+ return alg
+
+def recommendedAreaAndJVF_StandardJetGetter( finder, mainParam, input, oldMomentNames=False, **options):
+ """Run a jet algorithm, adding the recomended area, JVF and track related jet moments calculation.
+
+ all argument similar as in make_StandardJetGetter except :
+
+ oldMomentNames : bool, if True the moment naming will be the same as in 17.2.5 and previous releases.
+ """
+ return make_StandardJetAlg(finder, mainParam, input, jetMoments = ['ghostAssoc','area','jvf', 'trackMoments'], **options)
+ #alg = specialJet( finder, mainParam, input, oldMomentNames = oldMomentNames, moments=['ghostAssoc','area','jvf', 'trackMoments'], **options)
diff --git a/Reconstruction/Jet/JetMomentTools/python/SetupJetMomentTools.py b/Reconstruction/Jet/JetMomentTools/python/SetupJetMomentTools.py
new file mode 100644
index 00000000000..08115dbac9d
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/python/SetupJetMomentTools.py
@@ -0,0 +1,302 @@
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+from JetRec.JetAlgConfiguration import checkAndUpdateOptions
+from AthenaCommon.SystemOfUnits import *
+from JetRec.JetRecFlags import jetFlags
+from AthenaCommon.AlgSequence import AlgSequence
+
+class memoize(object):
+ """memoize is intended to be used as a decorator
+ It will force the decorated function to be run only once, cache it's return value
+ and return this value on every subsequent call.
+ """
+ def __init__(self, func):
+ # func is the function being decorated
+ self.func = func
+ self.cache = None
+
+ def __call__(self):
+ # this gives the memoize object the ability to be called
+ # and thus this object behaves like a function
+
+ if self.cache is None:
+ # call the function only once.
+ self.cache = self.func()
+ return self.cache
+
+
+
+def getJetVertexAssociationTool( toolName="JetVertexAssociation", *largs, **args):
+
+ #options=checkAndUpdateOptions(finder=finder, mainParam=mainParam, input=input,)
+
+ from JetMomentTools.JetMomentToolsConf import JetVertexAssociationTool, JetOriginCorrectionTool
+
+
+ jvaTool = JetVertexAssociationTool(toolName)
+
+ jvaTool.AssociateToHighestJVF = True
+
+ from JetRec.TrackSelectionForJets import getDefaultJetVtxTrackHelper
+ jvaTool.JetVtxTrackHelper = getDefaultJetVtxTrackHelper()
+
+ # origin correction tool
+ origCorr = JetOriginCorrectionTool("JetOriginCorrection")
+ origCorr.CellCalibrator = None
+ origCorr.UseGCW = False
+ origCorr.UseClusters = True
+
+ # Official Choice of Jet/MET Group : Primary Vertex
+ origCorr.UseJVA = False
+ origCorr.UsePrimaryVertex = True
+ origCorr.UseBeamSpot = True
+ jvaTool.OriginCorrectionTool = origCorr
+
+ jvaTool.StoreTrkMoments = False
+
+ jvaTool.sumPtTrkMomentName = "sumPtTrk_pv0_500MeV"
+ jvaTool.nTrkMomentName = "nTrk_pv0_500MeV"
+ jvaTool.TrackAssocName = "TrackAssoc"
+
+ return jvaTool
+
+def getJetTracksMomentTool(**options):
+
+ from JetMomentTools.JetMomentToolsConf import JetTracksMoments
+ from JetRec.TrackSelectionForJets import getDefaultJetVtxTrackHelper
+
+ suffix = options.get("MomentSuffix", None)
+ if suffix is None:
+ minpt = options.get("MinPt", 1*GeV)
+ if minpt>=GeV :
+ suffix = '_%1.1fGeV'%(minpt/GeV,)
+ suffix = suffix.replace('.0','') # don't write '1.0GeV', simply '1GeV'
+ else:
+ suffix = '_%dMeV'%(int(minpt),)
+ options["MomentSuffix"] = suffix
+ defaultJVtx = getDefaultJetVtxTrackHelper()
+ options.setdefault('JetVtxTrackHelper',defaultJVtx) # update the JetVtxTrackHelper option if not given
+ trkTool = JetTracksMoments(
+ **options
+ )
+ return trkTool
+
+
+def getORIGINMoments(jetcollname):
+ from JetCalibTools.MakeCalibSequences import alternateCalibSequence
+ from JetMomentTools.JetMomentToolsConf import JetMomentsFromCalib
+ from JetMomentTools.GhostAssociation import interpretJetName
+ from JetRec.JetAlgConfiguration import checkAndUpdateOptions
+ recognizedJet = interpretJetName(jetcollname)
+ if recognizedJet is not None:
+ finder, input, mainParam = recognizedJet
+ calib = 'LC:ORIGIN' if 'LC' in input else 'EM:ORIGIN'
+ d = checkAndUpdateOptions(input=input)
+ return JetMomentsFromCalib( "OriginRecalculator" , Calibrator = alternateCalibSequence(calib, d) )
+
+
+def getDefaultMissingCellTool(rmax=1.0):
+ rmax = max(rmax, 1.0)
+ if getDefaultMissingCellTool.cache is None:
+ from JetMomentTools.JetMomentToolsConf import MissingCellListTool
+ from AthenaCommon.AppMgr import ToolSvc
+
+ tool = MissingCellListTool()
+ tool.DeltaRmax = rmax
+ tool.AddBadCells = True
+ getDefaultMissingCellTool.cache = tool
+ tool.RemoveCellList = [1275120944, 1275121456, 1275137328, 1275268400, 1275268912, 1275383088, 1275383600, 1275514160, 1275514672, 1275628848, 1275661616, 1275662128, 1275792688, 1275793200, 1275940144, 1275940656, 1275956528, 1276054832, 1276055344, 1276071216, 1283493168, 1283493680, 1283640624, 1283641136, 1283820848, 1283821360, 1283935536, 1283936048, 1284017456, 1284017968, 1284164912, 1284165424, 1284312368, 1284312880, 1284427056, 1284427568, ]
+ #1275349024, 1283688480 # these tile cells were not ignored in the original implementation
+
+
+ ToolSvc += tool
+ else:
+ if rmax > getDefaultMissingCellTool.cache.DeltaRmax:
+ getDefaultMissingCellTool.cache.DeltaRmax = rmax
+
+ return getDefaultMissingCellTool.cache
+getDefaultMissingCellTool.cache = None
+
+
+def getJetBadChanCorrTool(finder, mainParam, input, **options):
+
+ options=checkAndUpdateOptions(finder=finder, mainParam=mainParam, input=input, **options)
+
+ from JetMomentTools.JetMomentToolsConf import JetBadChanCorrTool
+ from CaloClusterCorrection.StandardCellWeightCalib import getCellWeightTool
+
+ cellcalibtool = getCellWeightTool(options['finder'], options['mainParam'], input, onlyCellWeight=True)
+ coneSize = options['mainParam']
+
+ BadChanCorrT = JetBadChanCorrTool("JetBadChanCorrTool")
+ # Explicitly list flags
+ BadChanCorrT.ConeDr = coneSize
+ BadChanCorrT.CellCalibrator = cellcalibtool
+ BadChanCorrT.UseCalibScale = False
+ BadChanCorrT.AttachCorrCell = True
+ BadChanCorrT.AttachCorrDOTX = True
+ BadChanCorrT.AttachCorrJet = True
+
+ BadChanCorrT.MissingCellTool = getDefaultMissingCellTool(coneSize)
+ # Optional
+ BadChanCorrT.AttachCorrJetForCell = False
+ if options['addJetQualityMoments'] == True:
+ BadChanCorrT.AttachCorrJetForCell = True
+
+ # new setup for THistService
+ import os
+ dataPathList = os.environ[ 'DATAPATH' ].split(os.pathsep)
+ dataPathList.insert(0, os.curdir)
+ from AthenaCommon.Utils.unixtools import FindFile
+ RefFileName = FindFile( "JetBadChanCorrTool.root" ,dataPathList, os.R_OK )
+
+
+ from AthenaCommon.AppMgr import ServiceMgr
+ if not hasattr(ServiceMgr, 'THistSvc'):
+ from GaudiSvc.GaudiSvcConf import THistSvc
+ ServiceMgr += THistSvc()
+ ServiceMgr.THistSvc.Input += ["JetBadChanCorrTool DATAFILE=\'%s\' OPT=\'READ\'" % RefFileName]
+
+ return BadChanCorrT
+
+def getCaloQualityTool(computeFromCluster = False,doCellBasedVars = False,addJetQualityMoments=False):
+
+ from JetMomentTools.JetMomentToolsConf import JetCaloQualityTool
+ caloQualityT = JetCaloQualityTool("JetCaloQualityTool")
+ caloQualityT.doCellBasedVariables = doCellBasedVars
+ caloQualityT.cutOnLAr = 4000
+ caloQualityT.doLArQuality = True
+ caloQualityT.doTiming = True
+ caloQualityT.doHECQuality = True
+ caloQualityT.doNegativeE = True
+ caloQualityT.doTileQuality = False
+ caloQualityT.doAverageLArQF = True
+ #caloQualityT.timingMomentCells = [5,10] # produces jet moments: ootFracCells5, ootFracCells10 (names built on the fly from caloQualityT.timingMomentCellsName)
+ caloQualityT.timingMomentClusters = [5,10] #
+ caloQualityT.timingMomentCells = []
+ caloQualityT.doJetCentroid = True
+
+ caloQualityT.ComputeVariableFromCluster = computeFromCluster
+
+ if addJetQualityMoments :
+ caloQualityT.doConstituentBasedVariables = True
+ caloQualityT.doSamplingBasedVariables = True
+ caloQualityT.cutOnTile = 254
+ caloQualityT.doTileQuality = True
+
+
+ return caloQualityT
+
+
+def getJetIsolationTool(toolname,**options):
+ from JetMomentTools.JetMomentToolsConf import JetIsolationTool
+
+ default_options = dict(JetRParam = 0.4,
+ IsolationMoments = ["IsoKR","IsoDelta", "IsoFixedCone", "IsoFixedArea", "Iso6To8"],
+ IsolationParam = [ 2., 0.2, 0.8, 1.4, 0.8],
+ DoIsoSumPt = True,
+ )
+ default_options.update(options)
+
+
+ isot = JetIsolationTool(toolname, **default_options)
+ return isot
+
+@memoize
+def _getTruthJetIsoSelector():
+ from JetSimTools.JetSimToolsConf import JetTruthParticleSelectorTool
+ from AthenaCommon.AppMgr import ToolSvc
+
+ selector = JetTruthParticleSelectorTool("JetTruthIsoSelector")
+ selector.includeMuons = False
+ selector.useInteractingParticlesOnly = True
+
+ ToolSvc += selector
+ return selector
+
+def getTruthJetIsolationTool(toolname, doTruthPileup=False, **options):
+ """We need special options for truth jets.
+ """
+ isot = getJetIsolationTool(toolname, **options)
+
+ from RecExConfig.RecFlags import rec
+ defaultCollName = "SpclMC" if rec.readAOD() else "INav4MomTruthEvent"
+ defaultContList = [defaultCollName] if not doTruthPileup else [defaultCollName, defaultCollName+'INTIME']
+ truthCont = options.get('TruthInputContainers', defaultContList)
+
+ isot.TruthInputContainers = truthCont
+
+ # make sure we have a truth jet selector.
+ isot.TruthSelector = _getTruthJetIsoSelector()
+
+ return isot
+
+def getJetTruthLabelTool(truthCont='TruthInputContainers') :
+ from JetMomentTools.JetMomentToolsConf import JetTruthLabel
+
+ tool = JetTruthLabel("JetTruthLabel")
+ tool.TruthInputContainers = truthCont
+
+ return tool
+
+def getsubJetsMomentTool(toolname="KtDeltaR"):
+ from JetMomentTools.JetMomentToolsConf import JetSubJetsMoments
+ from JetRec.JetFastJetConfig import configFastJet
+
+ ktDrTool = JetSubJetsMoments(toolname)
+ fjtool = configFastJet("Kt",0.4, Inclusive=False, ExclusiveNjets=2, ExclusiveDcut=-1)
+ from AthenaCommon.AppMgr import ToolSvc
+ ToolSvc += fjtool
+ ktDrTool.FastJetTool = fjtool
+ ktDrTool.FromJetAssociation = "TrackAssoc"
+ ktDrTool.MinPt = 1*GeV
+ return ktDrTool
+
+
+
+def getMatchingFractionTool(jetcollname, matchType, seq=AlgSequence()):
+ """ ATTENTION this tools returns (tool, alg) where alg is a dependency for Matching.
+ It *has* to be scheduled before (caller's reponsability).
+ It return None if it can not guess what truth alg to request.
+ """
+ from JetMomentTools.GhostAssociation import interpretJetName
+ from JetMomentTools.JetMomentToolsConf import JetMatchingFractionTool
+ recognizedJet = interpretJetName(jetcollname)
+ if recognizedJet is None:
+ return None
+
+ finder, input, mainParam = recognizedJet
+ algtype = finder+str(int(mainParam*10))
+ from JetRec.JetGetters import make_StandardJetGetter
+ if matchType == "Truth":
+ from JetRec.JetRec_defaults import standardJetAlgorithm
+ alg = standardJetAlgorithm(algtype+"TruthJets")
+ if alg is None: # possible if truth collection already exist
+ # !!!! WARNNIG minPt cut is ambiguous !!!!
+ alg=make_StandardJetGetter(finder, mainParam,"Truth", minPt=4*GeV, jetMoments=[], outputCollectionName=algtype+"TruthNewJets", seq=seq).jetAlgorithmHandle()
+
+ # Then the mftool can be scheduled
+ mftool = JetMatchingFractionTool("TruthMF",MFCollection = alg.JetCollectionName, AssociatedConstitKey="TruthAssoc",
+ MomentTag="Truth",UseHighestShareFrac=True)
+ elif matchType == "Track":
+ alg=make_StandardJetGetter(finder, mainParam,"TrackZ", addJetQualityMoments=False, seq=seq).jetAlgorithmHandle()
+ # Then the mftool can be scheduled
+ mftool = JetMatchingFractionTool("TrackMF",MFCollection = alg.JetCollectionName, AssociatedConstitKey="TrackAssoc",
+ MomentTag="Track",UseHighestShareFrac=True)
+ else:
+ return None
+
+ return (mftool, alg)
+
+def getJetLArHVMomentTool():
+
+ from LArRecUtils.LArHVCorrToolDefault import LArHVCorrToolDefault
+ theLArHVCorrTool=LArHVCorrToolDefault()
+ from AthenaCommon.AppMgr import ToolSvc
+ ToolSvc+=theLArHVCorrTool
+
+ #Now declare the jet moment algorithm:
+ from JetMomentTools.JetMomentToolsConf import JetLArHVMoment
+ JetLArHV = JetLArHVMoment("JetLArHVMoment")
+ JetLArHV.HVCorrTool = theLArHVCorrTool
+ return JetLArHV
diff --git a/Reconstruction/Jet/JetMomentTools/src/JetBadChanCorrTool.cxx b/Reconstruction/Jet/JetMomentTools/src/JetBadChanCorrTool.cxx
new file mode 100644
index 00000000000..53a8335bebe
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/src/JetBadChanCorrTool.cxx
@@ -0,0 +1,461 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//#include "GaudiKernel/Service.h"
+#include "GaudiKernel/ITHistSvc.h"
+
+#include "PathResolver/PathResolver.h"
+
+#include "CaloDetDescr/CaloDetDescrManager.h"
+#include "CaloIdentifier/CaloCell_ID.h"
+#include "CaloEvent/CaloCellContainer.h"
+//#include "CaloEvent/CaloCellPrefetchIterator.h"
+//#include "CaloInterface/IHadronicCalibrationTool.h"
+
+#include "xAODCaloEvent/CaloCluster.h"
+#include "AthenaKernel/errorcheck.h"
+
+
+#include "TFile.h"
+#include "TH1.h"
+
+#include "JetMomentTools/JetBadChanCorrTool.h"
+#include "JetUtils/JetCellAccessor.h"
+#include "JetUtils/JetDistances.h"
+#include "GaudiKernel/SystemOfUnits.h"
+
+
+using Gaudi::Units::GeV;
+
+
+JetBadChanCorrTool::JetBadChanCorrTool( const std::string& name) :
+ JetModifierBase(name),
+ m_thistSvc( "THistSvc", name ),
+ // m_caloDDM(0),
+ // m_calo_id(0),
+ //m_calibTool(this),
+ // m_missingCellToolHandle("MissingCellListTool"),
+ m_missingCellMapName(),
+ m_forceMissingCellCheck(false),
+ m_useClusters(false)
+{
+
+ // limit to calculate the moments, #bad ch is typically 3000 out of 187616
+ declareProperty("NBadCellLimit", m_nBadCellLimit = 10000);
+
+
+ // for jet-level correction
+ declareProperty("StreamName", m_streamName = "/JetBadChanCorrTool/");
+ declareProperty("ProfileName", m_profileName = "JetBadChanCorrTool.root");
+ declareProperty("ProfileTag", m_profileTag = "");
+
+ declareProperty("UseCone", m_useCone = true);
+ //declareProperty("ConeDr", m_coneDr = 0.4);
+
+ // for calibration
+ declareProperty("UseCalibScale", m_useCalibScale = false);
+ // declareProperty("CellCalibrator", m_calibTool);
+
+ //declareProperty("MissingCellTool", m_missingCellToolHandle);
+ declareProperty("MissingCellMap", m_missingCellMapName= "MissingCaloCellsMap");
+ declareProperty("ForceMissingCellCheck", m_forceMissingCellCheck=false);
+ declareProperty("UseClusters",m_useClusters= false);
+}
+
+JetBadChanCorrTool::~JetBadChanCorrTool()
+{ }
+
+StatusCode JetBadChanCorrTool::initialize()
+{
+
+ if(!m_useClusters){
+ std::string fname = PathResolver::find_file(m_profileName, "DATAPATH");
+
+ if(fname==""){
+ ATH_MSG(ERROR) << "Could not get file " << m_profileName << endreq;
+ return StatusCode::FAILURE;
+ }
+ TFile tf(fname.c_str());
+ if(tf.IsOpen()==false){
+ ATH_MSG( ERROR ) << "Could not open file " << fname << endreq;
+ return StatusCode::FAILURE;
+ }
+ // already registered hists
+ std::vector<std::string> histsInSvc = m_thistSvc->getHists();
+
+ TListIter next(tf.GetListOfKeys());
+ TObject *obj;
+ while((obj = next())) {
+ // assume hname=(tag)_sample(sample)_pt(ptMin)_(ptMax)_eta(etaMin)_(etaMax)_phi(phiMin)_(phiMax)
+ std::string hname=obj->GetName();
+ if(hname.find('_')==std::string::npos) continue;
+ std::string tag = hname.substr(0,hname.find('_'));
+ std::string para = hname.substr(hname.find('_')+1);
+ std::string location = m_streamName+hname;
+
+ if(m_profileTag=="" || m_profileTag==tag){
+ // read in histo not registered yet
+ if(find(histsInSvc.begin(),histsInSvc.end(),location)==histsInSvc.end()){
+ StatusCode sc = m_thistSvc->regHist(location);
+ if(sc.isFailure()){
+ ATH_MSG( ERROR ) << "failed to read histo " << location << endreq;
+ return StatusCode::FAILURE;
+ }
+ }
+ int sample=0;
+ double ptMin=0;
+ double ptMax=9999;
+ double etaMin=0;
+ double etaMax=5.0;
+ double phiMin=-M_PI;
+ double phiMax=M_PI;
+ int ret = sscanf(para.c_str(),"sample%d_pt%lf_%lf_eta%lf_%lf_phi%lf_%lf",
+ &sample,&ptMin,&ptMax,&etaMin,&etaMax,&phiMin,&phiMax);
+
+ if(ret<1 || sample<0 || sample>=CaloCell_ID::Unknown) {
+ ATH_MSG( DEBUG ) << "Could not understand the name of hist " << obj->GetName() << endreq;
+ continue;
+ }
+
+ TH1* th=0;
+ StatusCode sc = m_thistSvc->getHist(location,th);
+ if(sc.isFailure()){
+ ATH_MSG( ERROR ) << "failed to get histo " << location << endreq;
+ return StatusCode::FAILURE;
+ }
+ m_profileDatas[sample].push_back(ProfileData((TH1D*)th,sample,ptMin,ptMax,etaMin,etaMax,phiMin,phiMax));
+ ATH_MSG( DEBUG ) << "read hist=" << th->GetName()
+ << " tag=" << tag << " sample=" << sample
+ << " ptMin=" << ptMin << " ptMax=" << ptMax
+ << " etaMin=" << etaMin << " etaMax=" << etaMax
+ << " phiMin=" << phiMin << " phiMax=" << phiMax << endreq;
+ }
+ }
+
+ // if(m_useCalibScale){
+ // CHECK( m_calibTool.retrieve() ) ;
+ // ATH_MSG( DEBUG ) << "use calibration " << m_calibTool << endreq;
+ // }
+
+ if ( m_useCone) {
+ // CHECK( detStore()->retrieve(m_caloDDM) );
+ // m_calo_id = m_caloDDM->getCaloCell_ID();
+ // ATH_MSG( DEBUG ) << "perform bad cells in cone calculations " << endreq;
+ // if(m_calo_id==0){
+ // ATH_MSG( ERROR ) << "Could not get CaloCell_ID" << endreq;
+ // return StatusCode::FAILURE;
+ // }
+ }
+
+
+ // CHECK( m_missingCellToolHandle.retrieve() );
+ //m_missingCellTool = dynamic_cast<MissingCellListTool*>( &(*m_missingCellToolHandle));
+
+ }
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode JetBadChanCorrTool::setupEvent()
+{
+
+ // CHECK(m_missingCellTool->prepareCellList());
+
+ return StatusCode::SUCCESS;
+}
+
+
+int JetBadChanCorrTool::modifyJet( xAOD::Jet& jet) const
+{
+
+ int res = 0;
+ if(m_useClusters){
+ res = correctionFromClustersBadCells( &jet);
+ }else{
+
+ const jet::CaloCellFastMap * badCellMap ;
+ StatusCode sc= evtStore()->retrieve(badCellMap, m_missingCellMapName) ;
+ if(sc.isFailure() ) {ATH_MSG_ERROR("Could not retieve bad cell map "<< m_missingCellMapName); return 1;}
+
+
+ // number of bad cells exceed the limit, set -1 and skip
+ if((int) badCellMap->cells().size() >m_nBadCellLimit){
+ jet.setAttribute<float>("BchCorrCell",-1.0);
+ if(!m_useClusters){
+ jet.setAttribute<float>("BchCorrDotx",-1.);
+ jet.setAttribute<float>("BchCorrJet",-1.);
+ jet.setAttribute<float>("BchCorrJetForCell",-1.);
+ }
+ return 0;
+ }
+
+
+ res = correctionFromCellsInJet(&jet, badCellMap);
+ // cone-based boundary
+ if(m_useCone) {
+ ATH_MSG_DEBUG( " Missing cells for cone search "<< badCellMap->size() << " jet="<< jet.index() << " R="<< jet.getSizeParameter() << " input="<< jet.getInputType() << " jet_eta="<<jet.eta() ) ;
+
+ res+=correctionFromCellsInCone(&jet, badCellMap);
+ }
+ }
+
+
+ // // for cryo.
+ // if(useCalibTool) {
+
+ // // without cell correction
+ // double emb3_raw = emb3-emb3_cell;
+ // double tile1_raw = tile1-tile1_cell;
+ // double cryoEt_raw = m_calibTool->etCryo(emb3_raw,tile1_raw);
+ // double emb3_rawotx = emb3-emb3_dotx;
+ // double tile1_rawotx = tile1-tile1_dotx;
+ // double cryoEt_rawotx = m_calibTool->etCryo(emb3_rawotx,tile1_rawotx);
+
+ // double cryoEt = m_calibTool->etCryo(emb3,tile1);
+ // corr_cell += (cryoEt-cryoEt_raw) / rawPt;
+ // corr_dotx += (cryoEt-cryoEt_rawotx) / rawPt;
+
+ // double cryoEt_jet_forcell = m_calibTool->etCryo(emb3_raw+emb3_jet_forcell,tile1_raw+tile1_jet_forcell);
+ // corr_jet_forcell += (cryoEt_jet_forcell-cryoEt_raw) / rawPt;
+
+ // double cryoEt_jet_associate = m_calibTool->etCryo(emb3_raw+emb3_jet_associate,tile1_raw+tile1_jet_associate);
+ // corr_jet_associate += (cryoEt_jet_associate-cryoEt_raw) / rawPt;
+
+ // double cryoEt_jet_cone = m_calibTool->etCryo(emb3_raw+emb3_jet_cone,tile1_raw+tile1_jet_cone);
+ // corr_jet_cone += (cryoEt_jet_cone-cryoEt_raw) / rawPt;
+
+ // total_energy += cryoEt*cosh(jet.eta());
+
+
+
+ // double calib_scale = total_energy/rawE;
+ // corr_cell /= calib_scale;
+ // corr_dotx /= calib_scale;
+ // corr_jet_associate /= calib_scale;
+ // corr_jet_cone /= calib_scale;
+ // corr_jet_forcell /= calib_scale;
+
+ // if(corr_cell <1e-7) corr_cell =0;
+ // if(corr_dotx <1e-7) corr_dotx =0;
+ // if(corr_jet_associate<1e-7) corr_jet_associate=0;
+ // if(corr_jet_cone <1e-7) corr_jet_cone =0;
+ // if(corr_jet_forcell <1e-7) corr_jet_forcell =0;
+ // }
+
+
+
+ return res;
+}
+
+
+int JetBadChanCorrTool::correctionFromCellsInJet( xAOD::Jet* jet, const jet::CaloCellFastMap * ) const {
+
+ //const jet::cellset_t & badcells = badCellMap.cells() ;
+
+ xAOD::JetFourMom_t p4 = jet->jetP4(xAOD::JetEMScaleMomentum);
+ double rawPt = p4.Pt();
+ double rawE = p4.E();
+
+ // jet moments, fraction
+ double corr_jet_associate=0;
+ //double corr_jet_cone=0;
+ double corr_jet_forcell=0;
+ double corr_cell=0;
+ double corr_dotx=0;
+
+ // jet energy
+ double total_energy=0;//for calib, with h1 weight but without jet-scale corr.
+
+ // for cryo., et at em-scale
+ // double emb3 = 0.;
+ // double tile1 = 0.;
+ // double emb3_cell = 0.;
+ // double tile1_cell = 0.;
+ // double emb3_dotx = 0.;
+ // double tile1_dotx = 0.;
+ // double emb3_jet_associate = 0.;
+ // double tile1_jet_associate = 0.;
+ // double emb3_jet_cone = 0.;
+ // double tile1_jet_cone = 0.;
+ // double emb3_jet_forcell = 0.;
+ // double tile1_jet_forcell = 0.;
+
+ // bool topoJet = dynamic_cast<const CaloCluster*>(*jet.begin()) != 0;
+
+ ///bool useCalibTool = m_useCalibScale && jet.signalState()==P4SignalState::CALIBRATED;
+
+ // const IMissingCellListTool::list_t & missingCells = m_missingCellTool->extendedList();
+
+ jet::JetCellAccessor::const_iterator cellIt = jet::JetCellAccessor::begin(jet);
+ jet::JetCellAccessor::const_iterator cellItE = jet::JetCellAccessor::end(jet);
+
+ for( ;cellIt!=cellItE; cellIt++) {
+
+ const CaloDetDescrElement * dde = (*cellIt)->caloDDE();
+ const CaloCell *cell = *cellIt;
+ double cellWeight = cellIt.weight();
+
+ double cell_energy = cell->e() * cellWeight;
+ // if(useCalibTool) {
+ // cell_energy *= m_calibTool->wtCell(cell);
+ // }
+
+ total_energy += cell_energy;
+
+ CaloCell_ID::CaloSample sampling = dde->getSampling();
+
+ // // for cryo.
+ // if (sampling == CaloCell_ID::EMB3) emb3 += cell->et();
+ // if (sampling == CaloCell_ID::TileBar0) tile1 += cell->et();
+
+ bool considerBad = cell->badcell();
+ //if(considerBad) badCellList.push_back( cell->ID() );
+ // if(m_forceMissingCellCheck) considerBad = considerBad || ( missingCells.find( cell->ID() ) == missingCells.end() );
+ // for bad channels in jet
+ if(considerBad){
+
+ // determine if this comes from a deadOTX
+ bool isOTX = false;
+ if(!(dde->is_tile()))
+ if(cell->provenance() & 0x0200)
+ isOTX = true;
+
+ double dr = jet::JetDistances::deltaR(jet->eta(),jet->phi(),dde->eta(), dde->phi());
+ double frac_cell = getProfile(rawPt, dr, sampling, dde->eta(), dde->phi());
+
+ double frac = frac_cell * cellWeight;
+ // if(useCalibTool) {
+ // frac *= m_calibTool->wtCell(cell);
+ // }
+
+ // for cryo.
+ // if (sampling == CaloCell_ID::EMB3) emb3_jet_associate += frac_cell * rawPt;
+ // if (sampling == CaloCell_ID::TileBar0) tile1_jet_associate += frac_cell * rawPt;
+
+ corr_jet_associate += frac;
+
+ // corrected in cell-level
+ if(cell_energy!=0){
+ // for cryo.
+ if(isOTX)
+ {
+ // if (sampling == CaloCell_ID::EMB3) emb3_dotx += cell->et();
+ // if (sampling == CaloCell_ID::TileBar0) tile1_dotx += cell->et();
+ corr_dotx += cell_energy ;
+ }
+ else
+ {
+ // if (sampling == CaloCell_ID::EMB3) emb3_cell += cell->et();
+ // if (sampling == CaloCell_ID::TileBar0) tile1_cell += cell->et();
+ corr_cell += cell_energy ;
+ }
+
+ // if (sampling == CaloCell_ID::EMB3) emb3_jet_forcell += frac_cell * rawPt;
+ // if (sampling == CaloCell_ID::TileBar0) tile1_jet_forcell += frac_cell * rawPt;
+
+ corr_jet_forcell += frac;
+
+ }
+ }
+ } // loop over cells
+
+ corr_cell /= rawE;
+ corr_dotx /= rawE;
+ ATH_MSG( DEBUG ) << "pt=" << jet->pt()/GeV << " eta=" << jet->eta() << " phi=" << jet->phi()
+ << " BCH_CORR_CELL=" << corr_cell
+ << " BCH_CORR_DOTX=" << corr_dotx
+ << " BCH_CORR_JET=" << corr_jet_associate
+ << " BCH_CORR_JET_FORCELL=" << corr_jet_forcell << endreq;
+
+ jet->setAttribute<float>("BchCorrCell",corr_cell);
+ jet->setAttribute<float>("BchCorrDotx",corr_dotx);
+ if(!m_useCone) jet->setAttribute<float>("BchCorrJet", corr_jet_associate);
+ jet->setAttribute<float>("BchCorrJetForCell",corr_jet_forcell);
+
+ return 0;
+}
+
+int JetBadChanCorrTool::correctionFromCellsInCone(xAOD::Jet* jet, const jet::CaloCellFastMap * badCellMap) const {
+
+ double corr_jet_cone=0;
+
+ double rawPt = jet->jetP4(xAOD::JetEMScaleMomentum).Pt();
+
+ // badCellList.insert( badCellList.end(), m_missingCellTool->begin() , m_missingCellTool->end() );
+
+
+ double jeteta = jet->eta();
+ double jetphi = jet->phi();
+ std::vector<jet::CellPosition> closeCells = badCellMap->cellsInDeltaR(jeteta,jetphi, jet->getSizeParameter() );
+ std::vector<jet::CellPosition>::iterator itr = closeCells.begin();
+ std::vector<jet::CellPosition>::iterator itrE = closeCells.end();
+
+ for(; itr!=itrE; itr++){
+
+ double cell_eta = itr->x();
+ double cell_phi = itr->phi();
+
+ // const CaloDetDescrElement * dde = m_caloDDM->get_element(itr->m_id);
+ // CaloCell_ID::CaloSample sampling = dde->getSampling();
+ CaloCell_ID::CaloSample sampling = itr->sampling();
+
+ double dr = jet::JetDistances::deltaR(jeteta,jetphi , cell_eta, cell_phi);
+ double frac_cell = getProfile(rawPt, dr, sampling, cell_eta, cell_phi);
+
+ //double frac = frac_cell;
+
+ // if(useCalibTool) {
+
+ // CaloCell temp(dde,frac_cell*rawE,0,0,CaloGain::UNKNOWNGAIN);//to pick up weight
+ // frac *= m_calibTool->wtCell(&temp);
+ // }
+
+ // if (sampling == CaloCell_ID::EMB3) emb3_jet_cone += frac_cell * rawPt;
+ // if (sampling == CaloCell_ID::TileBar0) tile1_jet_cone += frac_cell * rawPt;
+
+ corr_jet_cone += frac_cell;
+ }
+
+ jet->setAttribute<float>("BchCorrJet", corr_jet_cone );
+
+ return 0;
+}
+
+int JetBadChanCorrTool::correctionFromClustersBadCells( xAOD::Jet* jet) const {
+
+ double corrCell=0;
+
+ size_t nconstit=jet->numConstituents();
+ for(size_t i=0; i<nconstit; i++) {
+ // use the raw constituents since we're interested only in one moment :
+ const xAOD::IParticle* constit = jet->rawConstituent(i);
+ if( constit->type() != xAOD::Type::CaloCluster ) continue;
+ double badE;
+ bool v = static_cast<const xAOD::CaloCluster*>(constit)->retrieveMoment( xAOD::CaloCluster::ENG_BAD_CELLS, badE);
+ if(v) corrCell += badE;
+ }
+
+ double rawE = jet->jetP4(xAOD::JetEMScaleMomentum).E();
+ if(rawE==0)jet->setAttribute<float>("BchCorrCell", 0. );
+ else jet->setAttribute<float>("BchCorrCell", corrCell / rawE );
+ return 0;
+}
+
+
+
+
+
+double JetBadChanCorrTool::getProfile(double pt, double dr, int sample, double eta, double phi) const {
+ std::vector<ProfileData>::const_iterator itr =m_profileDatas[sample].begin();
+ std::vector<ProfileData>::const_iterator itrE =m_profileDatas[sample].end();
+ for(; itr!=itrE; itr++){
+ if((*itr).match(pt/GeV,sample,eta,phi)){
+ return (*itr).frac(dr);
+ }
+ }
+ return 0;
+}
+
+
+
diff --git a/Reconstruction/Jet/JetMomentTools/src/JetCaloEnergies.cxx b/Reconstruction/Jet/JetMomentTools/src/JetCaloEnergies.cxx
new file mode 100644
index 00000000000..06b160c4747
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/src/JetCaloEnergies.cxx
@@ -0,0 +1,57 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "JetMomentTools/JetCaloEnergies.h"
+
+#include "xAODJet/JetAccessorMap.h"
+
+#include "xAODCaloEvent/CaloCluster.h"
+#include "JetUtils/JetCaloQualityUtils.h"
+
+#include <vector>
+namespace {
+
+}
+
+JetCaloEnergies::JetCaloEnergies(const std::string & t) : JetModifierBase(t) {
+
+}
+
+
+
+int JetCaloEnergies::modifyJet( xAOD::Jet & jet)const {
+
+ static xAOD::JetAttributeAccessor::AccessorWrapper< std::vector<float> > & ePerSamplingAcc = *xAOD::JetAttributeAccessor::accessor< std::vector<float> >(xAOD::JetAttribute::EnergyPerSampling);
+
+
+ size_t numConstit = jet.numConstituents();
+ std::vector<float> & ePerSampling = ePerSamplingAcc( jet );
+ ePerSampling.resize(CaloSampling::Unknown, 0.);
+
+ for( float &e : ePerSampling) e=0.; // re-initialize
+
+ if( numConstit == 0 ) return 1;
+
+ // should find a more robust solution than using 1st constit type.
+ if( jet.rawConstituent( 0 )->type() == xAOD::Type::CaloCluster ){
+ // loop over raw constituents
+ for(size_t i=0;i<numConstit;i++){
+ for( size_t s= CaloSampling::PreSamplerB; s< CaloSampling::Unknown; s++){
+ float samplE = dynamic_cast<const xAOD::CaloCluster*>(jet.rawConstituent( i ))->eSample( (xAOD::CaloCluster::CaloSample) s);
+ //std::cout<< " ___ sampling "<< s << " adding "<< samplE << " to "<< ePerSampling[s] << std::endl;
+ ePerSampling[s] += samplE;
+ }
+ }
+
+ // for(float & cs: ePerSampling){
+ // ATH_MSG_DEBUG( jet.index() << " esampling = "<< cs ); }
+ // ATH_MSG_DEBUG("_________________________");
+ static xAOD::JetAttributeAccessor::AccessorWrapper<float>& emFracAcc = *xAOD::JetAttributeAccessor::accessor< float >(xAOD::JetAttribute::EMFrac);
+ static xAOD::JetAttributeAccessor::AccessorWrapper<float>& hecFracAcc = *xAOD::JetAttributeAccessor::accessor< float >(xAOD::JetAttribute::HECFrac);
+ emFracAcc(jet) = jet::JetCaloQualityUtils::emFraction( ePerSampling);
+ hecFracAcc(jet) = jet::JetCaloQualityUtils::hecF( &jet );
+ }
+
+ return 1;
+}
diff --git a/Reconstruction/Jet/JetMomentTools/src/JetCaloQualityTool.cxx b/Reconstruction/Jet/JetMomentTools/src/JetCaloQualityTool.cxx
new file mode 100644
index 00000000000..ed7b98a679d
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/src/JetCaloQualityTool.cxx
@@ -0,0 +1,177 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+
+#include "JetUtils/JetCaloQualityUtils.h"
+
+#include "JetMomentTools/JetCaloQualityTool.h"
+#include "xAODJet/JetAccessorMap.h"
+
+
+//#include "CaloDetDescr/CaloDepthTool.h"
+//#include "CaloIdentifier/CaloCell_ID.h"
+
+//#include "JetEvent/.h"
+
+#include <iostream>
+#include <iomanip>
+using namespace std;
+using namespace jet;
+
+
+JetCaloQualityTool::JetCaloQualityTool(const std::string& name)
+ : JetModifierBase(name)
+{
+
+ declareProperty("DoFracSamplingMax",m_doFracSamplingMax=false);
+
+ declareProperty("DoN90",m_doN90=true);
+
+ declareProperty("DoLArQuality", m_doLArQuality=true);
+ declareProperty("DoTileQuality", m_doTileQuality=false);
+ declareProperty("LArQualityCut", m_LArQualityCut=4000);
+ declareProperty("TileQualityCut", m_TileQualityCut=254);
+
+ declareProperty("DoTiming", m_doTiming=true);
+ declareProperty("TimingCellTimeCuts", m_timingCellTimeCuts);
+ declareProperty("TimingClusterTimeCuts", m_timingClusterTimeCuts);
+ declareProperty("TimingOnlyPosE", m_timingOnlyPosE=false);
+
+ declareProperty("DoNegativeE", m_doNegativeE=true);
+ declareProperty("DoHECQuality", m_doHECQuality=true);
+ declareProperty("DoAverageLArQF", m_doAverageLArQF=true);
+ declareProperty("DoJetCentroid", m_doJetCentroid=true);
+
+ declareProperty("ComputeVariableFromCluster",m_computeFromCluster=true);
+}
+
+
+int JetCaloQualityTool::modifyJet( xAOD::Jet& jet ) const
+{
+ //cout<<" "<<jet.getMomentKeys().size()<<endl;
+ //for(int i=0;i<jet.getMomentKeys().size();i++)
+ // cout<<jet.getMomentKeys()[i]<<endl;
+
+
+ if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Inside process() method" << endreq;
+
+ //==========================================================================
+ //JetCaloQualityUtils
+ //==========================================================================
+ if(m_doFracSamplingMax==true)
+ {
+ int sampling=-1;
+ double fracSamplingMax=JetCaloQualityUtils::fracSamplingMax(&jet,sampling);
+ ATH_MSG_VERBOSE("Setting " << xAOD::JetAttribute::FracSamplingMax << " to " << fracSamplingMax);
+ jet.setAttribute<float>(xAOD::JetAttribute::FracSamplingMax,fracSamplingMax );
+ jet.setAttribute<int>(xAOD::JetAttribute::FracSamplingMaxIndex,sampling );
+
+ // string samplingMaxName = "SamplingMax";
+ // ATH_MSG_VERBOSE("Setting " << samplingMaxName << " to " << sampling);
+ // jet.setAttribute<float>(samplingMaxName, sampling);
+
+ }
+
+ m_jetCalculations.process(&jet);
+ std::vector<double> cellmoments = m_jetCalculations.calculations();
+ std::vector<JetCaloCalculator*>& calculators = m_jetCalculations.calculators();
+ for(size_t i=0;i < calculators.size();i++) {
+ ATH_MSG_DEBUG( calculators[i]->name() << " -> "<< cellmoments[i] );
+ ATH_MSG_VERBOSE("Setting " << calculators[i]->name() << " to " << cellmoments[i]);
+ jet.setAttribute<float>( calculators[i]->name(), cellmoments[i] );
+ }
+
+ return StatusCode::SUCCESS;
+}
+
+
+
+StatusCode JetCaloQualityTool::initialize()
+{
+ if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Inside initialize() method" << endreq;
+
+
+ if(m_doN90)
+ {
+ JetCalcnLeadingCells* c = new jet::JetCalcnLeadingCells();
+ xAOD::JetAttribute::AttributeID id = m_computeFromCluster ? xAOD::JetAttribute::N90Constituents :
+ xAOD::JetAttribute::N90Cells ;
+
+ c->setName( xAOD::JetAttributeAccessor::name( id));
+ c->setCellLevelCalculation(! m_computeFromCluster);
+ m_jetCalculations.addCalculator( c );
+ }
+
+ if(m_doLArQuality)
+ {
+ JetCalcQuality* c = new jet::JetCalcQuality();
+ c->LArQualityCut = m_LArQualityCut;
+ c->setCellLevelCalculation(! m_computeFromCluster);
+ m_jetCalculations.addCalculator( c );
+ }
+
+ if(m_doTileQuality)
+ {
+ JetCalcQuality* c = new jet::JetCalcQuality();
+ c->TileQualityCut = m_TileQualityCut;
+ c->includeTile = true; c->includeLAr = false;
+ c->setCellLevelCalculation(! m_computeFromCluster);
+ m_jetCalculations.addCalculator( c );
+ }
+
+ if(m_doTiming)
+ {
+ JetCalcTimeCells* c = new JetCalcTimeCells();
+ c->setCellLevelCalculation(! m_computeFromCluster);
+ m_jetCalculations.addCalculator( c );
+ }
+
+ if(m_doHECQuality)
+ {
+ JetCalcQualityHEC* c = new JetCalcQualityHEC();
+ c->setCellLevelCalculation(! m_computeFromCluster);
+ m_jetCalculations.addCalculator( c );
+ }
+
+ if(m_doNegativeE)
+ {
+ JetCalcNegativeEnergy* c = new JetCalcNegativeEnergy();
+ c->setCellLevelCalculation(! m_computeFromCluster);
+ m_jetCalculations.addCalculator( c );
+ }
+
+ if(m_doAverageLArQF){
+ JetCalcAverageLArQualityF* c = new JetCalcAverageLArQualityF();
+ c->setCellLevelCalculation(! m_computeFromCluster);
+ m_jetCalculations.addCalculator( c );
+ }
+
+ if(m_doJetCentroid){
+ JetCalcCentroid* c = new jet::JetCalcCentroid();
+ c->setCellLevelCalculation(! m_computeFromCluster);
+ m_jetCalculations.addCalculator( c );
+ }
+
+ for( double & timeCut : m_timingClusterTimeCuts){
+
+ // build the moment name from the base-name and the value of the timing cut
+ std::stringstream s;
+ string timingMomentClustersName = "OotFracClusters";
+ s << std::setprecision(0) << std::fixed << timingMomentClustersName << timeCut;
+ std::string timingMomentClustersNameNow = s.str();
+
+ JetCalcOutOfTimeEnergyFraction* c = new jet::JetCalcOutOfTimeEnergyFraction();
+ c->setName(timingMomentClustersNameNow);
+ c->setCellLevelCalculation(false);
+ c->timecut = timeCut;
+ m_jetCalculations.addCalculator( c );
+ }
+ // cell timing to be done...
+ // for( double & timeCut : m_timingCellTimeCuts){
+
+ return StatusCode::SUCCESS;
+}
+
+
+
diff --git a/Reconstruction/Jet/JetMomentTools/src/JetIsolationTool.cxx b/Reconstruction/Jet/JetMomentTools/src/JetIsolationTool.cxx
new file mode 100644
index 00000000000..99f862eb51b
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/src/JetIsolationTool.cxx
@@ -0,0 +1,456 @@
+///////////////////////// -*- C++ -*- /////////////////////////////
+
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// JetIsolationTool.cxx
+// Implementation file for class JetIsolationTool
+///////////////////////////////////////////////////////////////////
+
+// JetMomentTools includes
+#include "JetMomentTools/JetIsolationTool.h"
+#include "xAODCaloEvent/CaloCluster.h"
+#include "JetUtils/JetDistances.h"
+#include <sstream>
+
+///////////////////////////////////////////////////////////////////
+// Define some isolation calculators.
+///////////////////////////////////////////////////////////////////
+namespace jet {
+
+ //////////////////////////////////////////////////////
+ /// \namespace JetIsolation
+ /// \brief Helper classes to calculate various isolation quantities
+ ///
+ /// There are many ways of calculating isolation variables. 3 aspect can be considered :
+ /// 1) The isolation zone (i.e which constituents not part of the jets are considered). In most of the cases this is simply a cone which radius depends on the jet.
+ /// 2) The kinematics computed from the contistuents in the isolation zone (sum pt, full P, ...)
+ /// 3) The signal sate at which constituents have to be chosen (i.e. LC scale or EM scale in the cluster case).
+ ///
+ /// To avoid to copy-paste similar code for the many possible combinations,
+ /// these 3 aspects are encapsulated in different classes :
+ /// 1) IsolationAreaBase and derived classes
+ /// 2) IsolationResult : which holds the 4-vector and has helpers to return the desired variables.
+ /// 3) KinematicGetter classes which hides the actual access to kinematics of constituents
+ ///
+ /// All these pieces are then combined in a IsolationCalculatorT<IsolationAreaBase, KinematicGetter> class (inheriting IsolationCalculator)
+ /// which is able to build the IsolationResult for the concrete IsolationAreaBase and KinematicGetter it receives as template paramters.
+ ///
+ /// Finally the JetIsolationTool instantiates as many IsolationCalculator as requested from its jobOption property.
+ ///
+
+ namespace JetIsolation {
+
+ typedef xAOD::IParticle::FourMom_t FourMom_t;
+
+ /// IsolationCalculator : base class for isolation calculations
+ class IsolationCalculator {
+ public:
+
+ /// Define the availble isolation variables
+ enum Kinematics {
+ SumPt, Par, Perp, P
+ };
+
+ /// names for isolation variables
+ static const std::string s_kname[4];
+
+ /// Holds the 4-vector of all constituents contributing to isolation.
+ /// accessors to corresponding isolation variables.
+ struct IsolationResult {
+ double isoSumPt(){ return m_isoSumPt;}
+ double isoPar(const TVector3 &unitV){ return m_isoP.Vect().Dot(unitV);}
+ double isoPerp(const TVector3 &unitV){ return m_isoP.Vect().Perp(unitV);}
+ FourMom_t & isoP(){ return m_isoP;}
+
+ FourMom_t m_isoP ;
+ double m_isoSumPt;
+
+ };
+
+
+ virtual ~IsolationCalculator() {}
+
+
+ virtual std::string baseName() {return "";}
+ virtual IsolationCalculator * clone(const xAOD::Jet* ){return NULL;}
+ virtual void copyFrom( const IsolationCalculator* o, const xAOD::Jet*){
+ m_attNames = o->m_attNames;
+ m_kinematics = o->m_kinematics;
+ }
+
+
+ /// compute the isolation momentum from jet and the constituents in nearbyConstit (these ones are supposed to be outside the jet).
+ virtual IsolationResult jetIsolation(const xAOD::Jet* , std::vector<jet::ParticlePosition> & ) const {return IsolationResult();}
+ virtual void setIsolationAttributes(xAOD::Jet* jet, std::vector<jet::ParticlePosition> & nearbyConstit) const {
+
+ IsolationResult result = jetIsolation(jet,nearbyConstit);
+ for(size_t i = 0; i< m_kinematics.size(); i++){
+ switch( m_kinematics[i]) {
+ case Perp: jet->setAttribute<float>( m_attNames[i], result.isoPerp(jet->p4().Vect()) ); break;
+ case SumPt: jet->setAttribute<float>( m_attNames[i], result.isoSumPt() ); break;
+ case Par: jet->setAttribute<float>( m_attNames[i], result.isoPar(jet->p4().Vect()) ); break;
+ case P: jet->setAttribute<float>( m_attNames[i], result.isoP().Vect().Mag() ); break;
+ }
+ }
+
+ }
+
+
+
+ bool scheduleKinematicCalculation(std::string kname){
+ for( size_t i=0; i<4 ;i++){
+ if( s_kname[i]==kname) {
+ m_kinematics.push_back( (Kinematics) i);
+ m_attNames.push_back( baseName()+kname );
+ return true;
+ }
+ }
+ return false;
+ }
+
+
+ void dump(){
+ std::cout << "Isolation calculator "<< baseName() <<std::endl;
+ for( size_t i =0 ; i<m_attNames.size(); i++){
+ std::cout << " - "<< m_attNames[i] << std::endl;
+ }
+ }
+
+ protected:
+ std::vector< Kinematics > m_kinematics; /// kinematics isolation vars to be computed
+ std::vector< std::string > m_attNames; /// kinematics isolation vars names
+ //std::string m_baseName;
+ };
+
+ const std::string IsolationCalculator::s_kname[4] = {"SumPt", "Par", "Perp", "P"};
+
+
+ struct IParticleKinematicGetter {
+ const FourMom_t& getP(jet::ParticlePosition& p) const { return p.particle()->p4();}
+ double getPt(jet::ParticlePosition& p) const { return p.particle()->pt();}
+ };
+
+ struct EMClusterKinematicGetter {
+ const FourMom_t& getP(jet::ParticlePosition& p) const {
+ const xAOD::CaloCluster * cl = static_cast<const xAOD::CaloCluster*>(p.particle() );
+ return cl->p4(xAOD::CaloCluster::UNCALIBRATED);
+ }
+ double getPt(jet::ParticlePosition& p) const {
+ const xAOD::CaloCluster * cl = static_cast<const xAOD::CaloCluster*>(p.particle() );
+ return cl->p4(xAOD::CaloCluster::UNCALIBRATED).Pt();
+ }
+ };
+
+
+ template<typename ISOCRITERIA, typename KINEMATICGETTER = IParticleKinematicGetter>
+ class IsolationCalculatorT : public IsolationCalculator {
+ public:
+
+ IsolationCalculatorT(double param=0.) : m_iso(param) { }
+
+ virtual IsolationResult jetIsolation(const xAOD::Jet* jet, std::vector<jet::ParticlePosition> &nearbyConstit) const {
+ IsolationResult result;
+ result.m_isoSumPt = 0;
+ int contributed = 0;
+ for(size_t i=0 ; i<nearbyConstit.size(); i++){
+ if( m_iso.inIsolationArea(jet, nearbyConstit[i]) ) {
+ result.m_isoP += m_kine.getP( nearbyConstit[i] ) ;
+ result.m_isoSumPt += m_kine.getPt(nearbyConstit[i]);
+ contributed++;
+ }
+ }
+ //std::cout << " "<< momentName()<< " contributed : "<< contributed<< " tot :"<< nearbyConstit.size() << " deltaRmax "<< m_deltaRmax << std::endl;
+ return result;
+ }
+
+ virtual std::string baseName() {return m_iso.name();}
+
+
+
+ virtual IsolationCalculator * clone(const xAOD::Jet* j){
+ IsolationCalculator* iso;
+ m_iso.setup(j);
+ if( j->getInputType() == xAOD::JetInput::EMTopo) {
+ auto* isoT = new IsolationCalculatorT<ISOCRITERIA, EMClusterKinematicGetter >();
+ isoT->m_iso = m_iso;
+ iso = isoT;
+ } else {
+ auto* isoT= new IsolationCalculatorT();
+ isoT->m_iso = m_iso;
+ isoT->m_kine = m_kine;
+ iso = isoT;
+ }
+ iso->copyFrom(this,j);
+ return iso;
+ }
+
+
+ ISOCRITERIA m_iso;
+ KINEMATICGETTER m_kine;
+ };
+
+
+ /// \class IsolationAreaBase
+ /// Defines a zone from which constituents will contribute to the isolation of a jet.
+ /// In most cases, the zone is simply a cone.
+ struct IsolationAreaBase {
+ IsolationAreaBase(double p, const std::string &n) : m_parameter(p), m_name(n){}
+
+ bool inIsolationArea(const xAOD::Jet* j, jet::ParticlePosition& part) const {
+
+ double dr2 = JetDistances::deltaR(j->eta(),j->phi(), part.x(), part.y());
+ return dr2< m_deltaRmax2;
+ }
+
+ std::string name(){
+ std::ostringstream oss; oss << m_name << int(10*m_parameter);
+ return oss.str();
+ }
+
+ double m_parameter;
+ std::string m_name;
+ double m_deltaRmax2;
+
+ };
+
+
+
+ // here we define a short cut to declare implementation of IsolationAreaBase classes.
+ // In most cases, only a definition of the max deltaR is enough to specify the area, hence these shortcuts.
+ // 1st parameter : class name
+ // 2nd parameter : code defining the max deltaR
+ // 3rd parameter (optional) : (re)-declaration of additional methods.
+ /// See below for example
+#define ISOAREA( calcName, deltaRcode , additionalDecl ) struct calcName : public IsolationAreaBase { \
+ calcName(double p) : IsolationAreaBase(p, #calcName){} \
+ virtual void setup(const xAOD::Jet* j){double jetRadius=j->getSizeParameter(); double param = m_parameter; m_deltaRmax2=deltaRcode ; m_deltaRmax2*=m_deltaRmax2; param=jetRadius*param;} additionalDecl }
+
+
+ ISOAREA( IsoKR , jetRadius*param, ) ;
+ ISOAREA( IsoDelta, jetRadius+param, ) ;
+ ISOAREA( IsoFixedCone, param, ) ;
+ ISOAREA( IsoFixedArea, sqrt(jetRadius*jetRadius+param/M_PI) , ) ;
+
+ // For Iso6To8 we need to redefine inIsolationArea
+ ISOAREA( Iso6To8, 0.8 , bool inIsolationArea(const xAOD::Jet* j, jet::ParticlePosition& constit)const ; ) ;
+ bool Iso6To8::inIsolationArea(const xAOD::Jet* j, jet::ParticlePosition& constit) const {
+ double dr2 = JetDistances::deltaR2(j->eta(), j->phi(), constit.x(), constit.y());
+ return ( (dr2<0.8*0.8) && (dr2>0.6*0.6) );
+ }
+
+
+ IsolationCalculator *createCalulator(std::string n, double parameter){
+
+ if( n == "IsoKR" ) return new IsolationCalculatorT<IsoKR>( parameter);
+ if( n == "IsoDelta" ) return new IsolationCalculatorT<IsoDelta>( parameter);
+ if( n == "IsoFixedCone" ) return new IsolationCalculatorT<IsoFixedCone>( parameter);
+ if( n == "IsoFixedArea" ) return new IsolationCalculatorT<IsoFixedArea>( parameter);
+ if( n == "Iso6To8" ) return new IsolationCalculatorT<Iso6To8>( parameter);
+ //if( n == "IsoCluster" ) return new IsoCluster(baseDr, parameter);
+
+ return NULL;
+ }
+
+
+ } // namespace JetIsolation
+
+/// helper
+ void colonSplit(const std::string &in, std::string &s1, std::string &s2, std::string &s3){
+ s2=""; s3="";
+ std::stringstream str(in);
+ std::getline(str, s1, ':' );
+ std::getline(str, s2, ':');
+ std::getline(str, s3);
+ }
+
+
+ bool isConstituent(const xAOD::IParticle* p,const xAOD::Jet* j){
+ // just need raw constituents.
+ for(size_t i = 0;i<j->numConstituents();i++){
+ if( p == j->rawConstituent(i) ) return true;
+ }
+ return false;
+ }
+
+} // namespace jet
+
+using namespace jet::JetIsolation;
+
+
+
+// Constructors
+////////////////
+JetIsolationTool::JetIsolationTool( const std::string& name)
+ : JetModifierBase(name),
+ m_deltaRmax(1.2),
+ m_rejectNegE(true)
+{
+
+ //
+ // Property declaration
+ //
+
+ declareProperty( "IsolationCalculations", m_isolationCodes);
+
+ declareProperty( "ConstituentContainer", m_inputPseudoJets = "", "The Constituents from which isolation is computed (will be retrieve automatically if blank)" );
+
+// declareProperty("RejectNegE", m_rejectNegE= true, "Ignore cluster with E<0");
+
+
+
+}
+
+// Destructor
+///////////////
+JetIsolationTool::~JetIsolationTool()
+{}
+
+// Athena algtool's Hooks
+////////////////////////////
+StatusCode JetIsolationTool::initialize()
+{
+ ATH_MSG_DEBUG ("Initializing " << name() << "...");
+
+ // a temporary map of isolation calculator
+ std::map<std::string, IsolationCalculator*> calcMap;
+
+ size_t nmom = m_isolationCodes.size();
+ /// decode each isolation calculations as entered by properties
+ for(size_t i=0;i <nmom; i++){
+ std::string isocriteria, param_s, kinematic;
+ jet::colonSplit( m_isolationCodes[i], isocriteria, param_s, kinematic);
+ std::string calcId = isocriteria+param_s;
+ IsolationCalculator* & isoC= calcMap[calcId];
+ if( isoC == NULL) {
+ isoC = createCalulator( isocriteria, std::stod( param_s)/10. );
+ }
+ if( isoC == NULL ) {
+ ATH_MSG_ERROR(" Unkown isolation criteria "<< isocriteria << " from "<< m_isolationCodes[i] );
+ return StatusCode::FAILURE;
+ }
+ bool ok = isoC->scheduleKinematicCalculation( kinematic);
+ if(!ok) {
+ ATH_MSG_ERROR(" Unkown isolation kinematic "<< kinematic << " from "<< m_isolationCodes[i] );
+ return StatusCode::FAILURE;
+ }
+ }
+
+ /// Fill the iso calculator vector from the map
+ for( auto & pair : calcMap ){
+ m_isoCalculators.push_back( pair.second );
+ ATH_MSG_DEBUG("Will use iso calculation : "<< pair.second->baseName() );
+ pair.second->dump();
+ }
+
+
+ return StatusCode::SUCCESS;
+}
+
+
+int JetIsolationTool::modify(xAOD::JetContainer& jets) const {
+
+ ATH_MSG_DEBUG("modify container...");
+ if(jets.empty() ) return 0;
+
+ const jet::ParticleFastMap * inputMap = retrieveInputMap();
+ if(inputMap == NULL) {
+ ATH_MSG_ERROR("Can't retrieve constituents container. No isolation calculations");
+ return 1;
+ }
+
+ ATH_MSG_DEBUG("modify : retrieved map");
+
+ // adapt the calculators to these jets (radius, input type, etc...)
+ // We use the 1st jet, assuming all others in the container are similar.
+ std::vector<IsolationCalculator*> calculators; // the adapted calculators.
+ for( auto * calc : m_isoCalculators ){
+ calculators.push_back( calc->clone( jets[0] ) );
+ }
+
+
+ /// Loop over jets in this collection.
+ for( xAOD::Jet* jet: jets ){
+
+ jet::ParticlePosition jetPos(jet);
+ std::vector<jet::ParticlePosition> all_nearbyC = inputMap->pointsInDr( jetPos , m_deltaRmax);
+
+ // restrict to nearby particles which are NOT jet constituents
+ // (for the expected num of constituents per jet (~30), doing a naive loop is
+ // faster than using set objects)
+ std::vector<jet::ParticlePosition> nearbyC;
+ nearbyC.reserve( 30 ); // educated guess.
+
+ for( size_t i=0;i<all_nearbyC.size(); i++){
+ jet::ParticlePosition &c= all_nearbyC[i];
+ if( jet::isConstituent( c.particle(), jet ) ) continue;
+ nearbyC.push_back( c ); // the constituent does not belong to jet
+ }
+
+ ATH_MSG_DEBUG( " Found nearbyC constit. All : "<< all_nearbyC.size() << " nearbyC :"<< nearbyC.size() << " in jet "<< jet->numConstituents() << " total "<< inputMap->size() );
+
+ // loop over calculators, calculate isolation given the close-by particles not part of the jet.
+ for( auto * calc : calculators ){
+ calc->setIsolationAttributes( jet, nearbyC );
+ }
+
+
+ }
+ ATH_MSG_DEBUG("done");
+ return 0;
+}
+
+
+const jet::ParticleFastMap * JetIsolationTool::retrieveInputMap() const {
+
+ ATH_MSG_DEBUG("retrieveInputMap ");
+
+ // Check if the map is in SG
+ std::string mapName = m_inputPseudoJets+"FastMap";
+ if( evtStore()->contains<jet::ParticleFastMap>(mapName) ){
+ const jet::ParticleFastMap * inputMap = 0;
+ StatusCode sc = evtStore()->retrieve( inputMap, mapName) ;
+ if(sc.isFailure()) return NULL; // should never happen
+ return inputMap;
+ }
+
+ ATH_MSG_DEBUG("Building the map ! ");
+
+ // else build the map !
+ jet::ParticleFastMap * inputMap = new jet::ParticleFastMap();
+ inputMap->init(m_deltaRmax);
+ // For now very simple and limited.
+ const xAOD::IParticleContainer* inputConstits=0;
+ StatusCode sc= evtStore()->retrieve( inputConstits, m_inputPseudoJets) ;
+ if(sc.isFailure()) {
+ ATH_MSG_ERROR("No input constituents "<< m_inputPseudoJets);
+ return NULL;
+ }
+
+ for( const xAOD::IParticle* p : *inputConstits ){
+ if( m_rejectNegE && ( p->e()<0) ) continue;
+ jet::ParticlePosition pos( p );
+ inputMap->insert( pos );
+ }
+
+ ATH_MSG_DEBUG("Built map size= "<< inputMap->size());
+
+ sc= evtStore()->record( inputMap, mapName) ;
+ if(sc.isFailure()) { ATH_MSG_WARNING( "Can't record fatmap "<< mapName); }
+ return inputMap;
+}
+
+
+
+
+
+StatusCode JetIsolationTool::finalize()
+{
+ ATH_MSG_DEBUG ("Finalizing " << name() << "...");
+ for( size_t i=0;i<m_isoCalculators.size(); i++){
+ delete m_isoCalculators[i];
+ }
+ return StatusCode::SUCCESS;
+}
+
diff --git a/Reconstruction/Jet/JetMomentTools/src/JetLArHVTool.cxx b/Reconstruction/Jet/JetMomentTools/src/JetLArHVTool.cxx
new file mode 100644
index 00000000000..2697e970b07
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/src/JetLArHVTool.cxx
@@ -0,0 +1,105 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "JetUtils/JetDistances.h"
+
+#include "JetMomentTools/JetLArHVTool.h"
+#ifdef ASGTOOL_ATHENA
+#include "JetUtils/JetCellAccessor.h"
+#endif
+
+JetLArHVTool::JetLArHVTool(const std::string& name)
+ : JetModifierBase(name),
+ m_useCells(true),
+#ifdef ASGTOOL_ATHENA
+ m_hvCorrTool("LArHVCorrTool")
+#endif
+{
+ declareProperty("UseCells", m_useCells=true);
+ declareProperty("HVCorrTool",m_hvCorrTool);
+#ifdef ASGTOOL_ATHENA
+ declareProperty("keyHVScaleCorr",m_keyHVScaleCorr="LArHVScaleCorr");
+#endif
+}
+
+
+StatusCode JetLArHVTool::initialize()
+{
+// Retrieve HVCorrTool
+ StatusCode sc = m_hvCorrTool.retrieve();
+ if (sc.isFailure()) {
+ ATH_MSG_ERROR ( "Unable to find tool for LArHVCorrTool");
+ return StatusCode::FAILURE;
+ }
+
+
+ return StatusCode::SUCCESS;
+
+}
+
+
+int JetLArHVTool::modifyJet( xAOD::Jet& jet ) const
+{
+ double energy=0;
+ double energyHVaff=0;
+ int numCells=0;
+ int numCellsHVaff=0;
+
+ double eFrac = 0;
+ double nFrac = 0;
+
+
+
+ if( m_useCells){
+#ifdef ASGTOOL_ATHENA
+ jet::JetCellAccessor::const_iterator it = jet::JetCellAccessor::begin(&jet);
+ jet::JetCellAccessor::const_iterator itE = jet::JetCellAccessor::end(&jet);
+ for( ;it!=itE; it++) {
+ //e += it->e()*it.weight();
+ const CaloCell* thisCell = *it;
+ Identifier cellID=thisCell->ID();
+ numCells++;
+ energy+=fabs(thisCell->e());
+
+ //keep only the LAr cells in the jet:
+ if(thisCell->caloDDE()->is_tile())continue;
+
+ //float hvdev = 0;
+
+ //retrieve offline correction from DB:
+ float hvcorr = m_hvCorrTool->Scale(cellID);
+
+ // //retrieve online correction from MetaData (available only in raw data)
+ // float hvonline = m_dd_HVScaleCorr->HVScaleCorr(cellID);
+
+ //Correction should be between (0 and 2)
+ if (!(hvcorr>0. && hvcorr<100.)) continue;
+
+
+ //check if the difference between correction and 1 is greater than 2 per thousand, then we are in a LAr HV affected area.
+ if(fabs(hvcorr)-1.>0.002){
+ ATH_MSG_DEBUG ( "HV affected" << hvcorr );
+ energyHVaff+=fabs(thisCell->e());
+ numCellsHVaff++;
+ }//end of non nominal HV area
+
+#endif
+ }//cells
+
+ if(energy>0) eFrac = energyHVaff/energy;
+ if(numCells>0) nFrac = double(numCellsHVaff)/double(numCells);
+
+ } // use cells
+
+
+
+ // set the attributes
+ jet.setAttribute<float>("LArBadHVEnergyFrac", eFrac);
+ jet.setAttribute<float>("LArBadHVNCellFrac", nFrac);
+
+
+ return 0;
+}
+
+
diff --git a/Reconstruction/Jet/JetMomentTools/src/components/JetMomentTools_entries.cxx b/Reconstruction/Jet/JetMomentTools/src/components/JetMomentTools_entries.cxx
new file mode 100644
index 00000000000..a5fefd200f9
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/src/components/JetMomentTools_entries.cxx
@@ -0,0 +1,38 @@
+// JetRec_entries.cxx
+#include "GaudiKernel/DeclareFactoryEntries.h"
+
+#include "JetMomentTools/JetCaloEnergies.h"
+#include "JetMomentTools/JetCaloQualityTool.h"
+#include "JetMomentTools/JetWidthTool.h"
+#include "JetMomentTools/JetBadChanCorrTool.h"
+#include "JetMomentTools/JetVertexFractionTool.h"
+#include "JetMomentTools/JetTrackMomentsTool.h"
+#include "JetMomentTools/JetMuonSegmentMomentsTool.h"
+#include "JetMomentTools/JetPtAssociationTool.h"
+#include "JetMomentTools/JetIsolationTool.h"
+#include "JetMomentTools/JetLArHVTool.h"
+
+DECLARE_TOOL_FACTORY(JetCaloEnergies)
+DECLARE_TOOL_FACTORY(JetCaloQualityTool)
+DECLARE_TOOL_FACTORY(JetWidthTool)
+DECLARE_TOOL_FACTORY(JetBadChanCorrTool)
+DECLARE_TOOL_FACTORY(JetVertexFractionTool)
+DECLARE_TOOL_FACTORY(JetTrackMomentsTool)
+DECLARE_TOOL_FACTORY(JetMuonSegmentMomentsTool)
+DECLARE_TOOL_FACTORY(JetPtAssociationTool)
+DECLARE_TOOL_FACTORY(JetIsolationTool)
+DECLARE_TOOL_FACTORY(JetLArHVTool)
+
+DECLARE_FACTORY_ENTRIES(JetRec) {
+ DECLARE_TOOL(JetCaloEnergies)
+ DECLARE_TOOL(JetCaloQualityTool)
+ DECLARE_TOOL(JetWidthTool)
+ DECLARE_TOOL(JetBadChanCorrTool)
+ DECLARE_TOOL(JetMuonSegmentMomentsTool)
+ DECLARE_TOOL(JetTrackMomentsTool)
+ DECLARE_TOOL(JetVertexFractionTool)
+ DECLARE_TOOL(JetPtAssociationTool)
+ DECLARE_TOOL(JetIsolationTool)
+ DECLARE_TOOL(JetLArHVTool)
+}
+
diff --git a/Reconstruction/Jet/JetMomentTools/src/components/JetMomentTools_load.cxx b/Reconstruction/Jet/JetMomentTools/src/components/JetMomentTools_load.cxx
new file mode 100644
index 00000000000..0e6ec8f4140
--- /dev/null
+++ b/Reconstruction/Jet/JetMomentTools/src/components/JetMomentTools_load.cxx
@@ -0,0 +1,5 @@
+// JetRec_load.cxx
+
+#include "GaudiKernel/LoadFactoryEntries.h"
+
+LOAD_FACTORY_ENTRIES(JetMomentTools)
--
GitLab