From 87854fbd46e911a4bab5f114446c8662e06f6164 Mon Sep 17 00:00:00 2001 From: Pierre-Antoine Delsart <delsart@in2p3.fr> Date: Fri, 5 Dec 2014 16:31:39 +0100 Subject: [PATCH] 'JetLArHVTool now uses cluster moments. ATLJETMET-15' (JetMomentTools-00-02-42) 2014-12-05 <delsart@localhost> * JetLArHVTool : now uses cluster moments. Moved to Root/ * ATLJETMET-15 * JetIsolationTool moved to Root * JetMomentTools-00-02-42 2014-12-04 <delsart@localhost> * JetOriginCorrectionTool also sets the "OriginVertex" associated objects. * ATLJETMET-146 * JetMomentTools-00-02-41 2014-12-02 David Adams * JetMomentTools-00-02-40 * ATLJETMET-134 * Root/JetVertexFractionTool.cxx - Change attribute HighestJVFVertex to replace JVF with the property JVFName. ... (Long ChangeLog diff - truncated) --- .../JetMomentTools/JetCaloEnergies.h | 5 + .../JetMomentTools/JetCaloQualityTool.h | 47 ++-- .../JetMomentTools/JetIsolationTool.h | 26 +- .../JetMomentTools/JetLArHVTool.h | 9 - .../JetMomentTools/JetOriginCorrectionTool.h | 48 ++++ .../JetMomentTools/JetVertexFractionTool.h | 66 +++-- .../JetMomentTools/JetVertexTaggerTool.h | 125 +++++++++ .../JetMomentTools/Root/JetCaloEnergies.cxx | 117 ++++++++ .../Root/JetCaloQualityTool.cxx | 124 +++++++++ .../{src => Root}/JetIsolationTool.cxx | 35 ++- .../Jet/JetMomentTools/Root/JetLArHVTool.cxx | 46 +++ .../Root/JetOriginCorrectionTool.cxx | 49 ++++ .../Root/JetVertexFractionTool.cxx | 137 +++++---- .../Root/JetVertexTaggerTool.cxx | 263 ++++++++++++++++++ .../Jet/JetMomentTools/cmt/Makefile.RootCore | 2 +- .../Jet/JetMomentTools/cmt/requirements | 16 +- .../python/JetMomentsConfigHelpers.py | 1 + .../share/JVTlikelihood_20140805.root | Bin 0 -> 19745 bytes .../src/JetCaloCellQualityTool.cxx | 87 ++++++ .../src/JetCaloCellQualityTool.h | 53 ++++ .../JetMomentTools/src/JetCaloEnergies.cxx | 57 ---- .../JetMomentTools/src/JetCaloQualityTool.cxx | 177 ------------ .../Jet/JetMomentTools/src/JetLArHVTool.cxx | 105 ------- .../src/components/JetMomentTools_entries.cxx | 9 + 24 files changed, 1135 insertions(+), 469 deletions(-) create mode 100644 Reconstruction/Jet/JetMomentTools/JetMomentTools/JetOriginCorrectionTool.h create mode 100644 Reconstruction/Jet/JetMomentTools/JetMomentTools/JetVertexTaggerTool.h create mode 100644 Reconstruction/Jet/JetMomentTools/Root/JetCaloEnergies.cxx create mode 100644 Reconstruction/Jet/JetMomentTools/Root/JetCaloQualityTool.cxx rename Reconstruction/Jet/JetMomentTools/{src => Root}/JetIsolationTool.cxx (93%) create mode 100644 Reconstruction/Jet/JetMomentTools/Root/JetLArHVTool.cxx create mode 100644 Reconstruction/Jet/JetMomentTools/Root/JetOriginCorrectionTool.cxx create mode 100644 Reconstruction/Jet/JetMomentTools/Root/JetVertexTaggerTool.cxx create mode 100644 Reconstruction/Jet/JetMomentTools/share/JVTlikelihood_20140805.root create mode 100644 Reconstruction/Jet/JetMomentTools/src/JetCaloCellQualityTool.cxx create mode 100644 Reconstruction/Jet/JetMomentTools/src/JetCaloCellQualityTool.h delete mode 100644 Reconstruction/Jet/JetMomentTools/src/JetCaloEnergies.cxx delete mode 100644 Reconstruction/Jet/JetMomentTools/src/JetCaloQualityTool.cxx delete mode 100644 Reconstruction/Jet/JetMomentTools/src/JetLArHVTool.cxx diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetCaloEnergies.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetCaloEnergies.h index c7ee80fc727..97e707b1270 100644 --- a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetCaloEnergies.h +++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetCaloEnergies.h @@ -8,6 +8,7 @@ #define JETMOMENTTOOLS_JETCALOENERGIES_H #include "JetRec/JetModifierBase.h" +#include <vector> class JetCaloEnergies : public JetModifierBase { ASG_TOOL_CLASS0(JetCaloEnergies); @@ -18,7 +19,11 @@ public: virtual int modifyJet(xAOD::Jet& ) const ; +protected: + void fillEperSamplingCluster(xAOD::Jet &jet, std::vector<float> & ePerSampling ) const ; + void fillEperSamplingPFO(xAOD::Jet &jet, std::vector<float> & ePerSampling ) const ; + }; diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetCaloQualityTool.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetCaloQualityTool.h index cfeab4afea8..39a155c3089 100644 --- a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetCaloQualityTool.h +++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetCaloQualityTool.h @@ -4,7 +4,6 @@ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration */ - /** @class JetCaloQualityTool Calculates calorimeter based variables for jet quality @@ -12,12 +11,20 @@ @author P-A Delsart @date (first implementation) October , 2009 @date (run 2 implementation) February , 2014 -*/ + This tool calculate calorimeter quantities from a jet and saves them as attributes. + It takes a list of string as a property 'Calculations' : the list of moments to calculate.<br> + From this list it internally prepares a list of calculators from JetUtils/JetCaloQualityUtils.h<br> + + List of known calculations (see also the implementation of initialize()) : + LArQuality, Timing, NegativeE, Centroid, N90Constituents, BchCorrCell, FracSamplingMax + + This class performs cluster-based calculation. For similar cell-based calculation, see JetCaloCellQualityUtils.h + +*/ #ifndef JETREC_JETCALOQUALITYTOOL_H #define JETREC_JETCALOQUALITYTOOL_H - #include "JetRec/JetModifierBase.h" #include "JetUtils/JetCaloCalculations.h" @@ -28,7 +35,7 @@ class JetCaloQualityTool: public JetModifierBase { ASG_TOOL_CLASS0(JetCaloQualityTool); - + public: JetCaloQualityTool(const std::string & name); @@ -36,31 +43,21 @@ public: virtual StatusCode initialize(); - private: - - int m_LArQualityCut; - int m_TileQualityCut; + protected: + /// Names of calo quantities to compute and to add as attributes + std::vector<std::string> m_calculationNames; - 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; + /// Time cuts for Out-of-time calo quantities. + std::vector <double> m_timingTimeCuts; + /// This objects holds a list of cluster-based calculators + jet::JetCaloCalculations m_jetCalculations; - std::vector <double> m_timingCellTimeCuts; - std::vector <double> m_timingClusterTimeCuts; + // internal pointer to m_jetCalculations (this pointer is also used in the cell-based derived tool) + jet::JetCaloCalculations * m_calcProcessor; - mutable jet::JetCaloCalculations m_jetCalculations; - - + + bool m_doFracSamplingMax; // internal }; #endif diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetIsolationTool.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetIsolationTool.h index f1bbd89509f..4873b96b89f 100644 --- a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetIsolationTool.h +++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetIsolationTool.h @@ -18,7 +18,7 @@ /////////////////////////////////////////////////////////////////// /// \class JetIsolationTool /// -/// This JetModifier calculates isolation variables for jets. +/// Calculate isolation variables for jets and set them as jet attributes. /// /// The isolation variables are calculated from a list of input particles /// which are close to the jet but not part of its constituents. @@ -30,10 +30,13 @@ /// 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) +/// - ISOCRITERIA : encode the method used to find isolation particles (IsoKR, IsoDelta, IsoFixedCone, IsoFixedArea, Iso6To8, signification detailed below) /// - 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) +/// - VAR : encode the variable to be calculated (Per, Par, SumPt, P, significations detailed below) /// +/// The resulting attribute name is simply the identifier string with ':' removed. +/// Example "IsoKR:11:Perp" -> "IsoKR11Perp" +/// /// 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) @@ -43,10 +46,23 @@ /// - 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) +/// These are solved at the cost of increased complexity (internal helper classes, usage of a special fast lookup map), details in the .cxx file. /// /// WARNING : currently works well only for LCTopoJets, TrackJets, TruthJets AND small R jets (R ~<0.8) /// +/// +/// Isolation Criteria ("param" below is the main parameter) : +/// - IsoKR : iso area == cone of size __jetRadius*param__ +/// - IsoDelta : iso area == cone of size __jetRadius+param__ +/// - IsoFixedCone : iso are == cone of size __param__ +/// - IsoFixedArea : iso are == cone of size __sqrt(jetRadius*jetRadius+param/M_PI)__ +/// - Iso6To8 : iso area == annulus between R=0.6 and R=0.8 +/// +/// Isolation Variables (iso4vec = sum of 4-vec in isolation area, not part of jet) +/// - Perp : iso4Vec.Perp( jetDirection ) +/// - Par : iso4Vec.Dot( jetDirection ) +/// - SumPt : sum Pt of 4-vec contibuting to iso4Vec +/// - P : iso4Vec.Vect().Mag() /////////////////////////////////////////////////////////////////////// namespace jet { @@ -133,6 +149,8 @@ protected: float m_deltaRmax; std::string m_inputPseudoJets; + /// the list of isolation calculation object (they are actually used only as template objects + /// from which the actual calculators are build & adapted to the jet object, see implementation) std::vector<jet::JetIsolation::IsolationCalculator*> m_isoCalculators; diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetLArHVTool.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetLArHVTool.h index a9abdf5dd41..2f71d9d8c55 100644 --- a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetLArHVTool.h +++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetLArHVTool.h @@ -21,8 +21,6 @@ #include "AsgTools/ToolHandle.h" -#include "LArElecCalib/ILArHVCorrTool.h" - class JetLArHVTool: public JetModifierBase { ASG_TOOL_CLASS0(JetLArHVTool); @@ -35,13 +33,6 @@ public: private: - bool m_useCells; - std::string m_channelKey; - std::string m_keyHVScaleCorr; - -#ifdef ASGTOOL_ATHENA - ToolHandle<ILArHVCorrTool> m_hvCorrTool; -#endif }; diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetOriginCorrectionTool.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetOriginCorrectionTool.h new file mode 100644 index 00000000000..90598f7134b --- /dev/null +++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetOriginCorrectionTool.h @@ -0,0 +1,48 @@ +// JetOriginCorrectionTool.h -*- C++ -*- + +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef JetMomentTools_JetOriginCorrectionTool_H +#define JetMomentTools_JetOriginCorrectionTool_H + +//////////////////////////////////////////////// +/// \file JetOriginCorrectionTool.h +/// \class JetOriginCorrectionTool +/// +/// This tool computes the constituent scale 4-vector of the jet +/// relative to a PV vertex and stores it in as an attribute +/// of the jet. +/// This tool DOES NOT change the scale of the jet, just adds an attribute. +/// +/// This tool is a wrapper : it only calls functions defined in JetUtils/JetOriginHelpers.h +/// +/// Properties: +/// - VertexContainer: Name of the vertex container. The 0th vertex is used to compute the correction +/// - OriginCorrectedName: Name of the output attribute +//////////////////////////////////////////////////////////// + +#include "AsgTools/AsgTool.h" +#include "JetInterface/IJetModifier.h" + +class JetOriginCorrectionTool : public asg::AsgTool, + virtual public IJetModifier { + ASG_TOOL_CLASS(JetOriginCorrectionTool, IJetModifier) + +public: + + /// Constructor from tool name. + JetOriginCorrectionTool(const std::string& myname); + + /// Inherited method to modify a jet. Compute the origin-corrected + /// momentum and put it in the jet + virtual int modify(xAOD::JetContainer& jet) const; + + protected: + + std::string m_vtxContainerName; + std::string m_correctionName; +}; + +#endif diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetVertexFractionTool.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetVertexFractionTool.h index 9ef955bdf85..b6ba71f2cc3 100644 --- a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetVertexFractionTool.h +++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetVertexFractionTool.h @@ -22,46 +22,70 @@ /// 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 +/// +/// Properties: +/// VertexContainer - name of the vertex container +/// AssociatedTracks - name for attribute holding the list of associated tracks +/// TrackVertexAssociation - name for the container holding the track-vertex associations +/// TrackSelector - tool to select tracks (none ==> no selection) +/// JVFName - name for the JVF array attribute (default is "JVF") + + +#include "AsgTools/ToolHandle.h" #include "xAODTracking/Vertex.h" #include "xAODTracking/VertexContainer.h" #include "xAODTracking/TrackParticle.h" #include "xAODTracking/TrackParticleContainer.h" +#include "JetInterface/IJetTrackSelector.h" #include "JetRec/JetModifierBase.h" #include "JetEDM/TrackVertexAssociation.h" #include <vector> #include <string> -class JetVertexFractionTool : public JetModifierBase -{ - ASG_TOOL_CLASS(JetVertexFractionTool,IJetModifier); +class JetVertexFractionTool : public JetModifierBase { + ASG_TOOL_CLASS(JetVertexFractionTool,IJetModifier); + +public: // methods + + // Constructor from tool name + JetVertexFractionTool(const std::string& name); + + // Initialization. + StatusCode initialize(); + + // Inherited methods to modify a jet + // Calls getJetVertexFraction and puts the result in the jet + virtual int modifyJet(xAOD::Jet& jet) const; - public: - // Constructor from tool name - JetVertexFractionTool(const std::string& name); + // 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; - // 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 the JVF for a given vertex + float getJetVertexFraction(const xAOD::Vertex*, + const std::vector<const xAOD::TrackParticle*>&, + const jet::TrackVertexAssociation*) 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 determine the highest JVF vertex and get an ElementLink to it + ElementLink<xAOD::VertexContainer> getMaxJetVertexFraction(const xAOD::VertexContainer*, + const std::vector<float>&) 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; +private: // data - // 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; + // Configurable parameters + std::string m_verticesName; + std::string m_assocTracksName; + std::string m_tvaName; + ToolHandle<IJetTrackSelector> m_htsel; + std::string m_jvfname; - private: - // Configurable parameters - std::string m_verticesName; - std::string m_assocTracksName; - std::string m_tvaName; +private: // methods - std::vector<float> getEmptyJetVertexFraction(const xAOD::VertexContainer*) const; + std::vector<float> getEmptyJetVertexFraction(const xAOD::VertexContainer*) const; }; diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetVertexTaggerTool.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetVertexTaggerTool.h new file mode 100644 index 00000000000..34f80a2456f --- /dev/null +++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetVertexTaggerTool.h @@ -0,0 +1,125 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +// JetVertexTaggerTool.h + +#ifndef JETMOMENTTOOLS_JETVERTEXTAGGERTOOL_H +#define JETMOMENTTOOLS_JETVERTEXTAGGERTOOL_H + +/// James Frost \n +/// November 2014 +/// +/// Tool to calculate the jet vertex tag (JVT) +/// JVT is a float per jet +/// +/// The JVT likelihood is parameterised in terms of: +/// RPT: ratio of the primary vertex track sum to the jet pT +/// JVFCorr: a corrected JVF calculation accounting for the number of PU tracks in the event. +/// +/// Calculation requires three main types of information +/// 1. Vertex container for the event (from evtStore), with respect to which the JVT track sums +/// 2. Tracks associated to each of the input jet (in the jet aux store) +/// 3. Track vertex association object (from evtStore) +/// 4. The track container needed for PU track counting +/// 5. The ID track selector tool to provide an input track selection. +/// 6. An input file with the JVT likelihood stored as a 2D histogram. +/// +/// Using this information, the procedure can be broken into three main steps: +/// 1. Count the total number of pileup tracks +/// 2. Per jet: Get track pT sums needed for JVFCorr and RpT +/// 3. Use JVT likelihood to get JVT value +/// +/// Attributes added (each is a float per jet): +/// 1. Jvt - the JVTLikelihood value for the jet. This is parameterised in terms of: +/// 2. JvtRpt - the ratio of the PV track pT sum to the jet pT +/// 3. JvtJvfcorr - a corrected Jvf quantity that corrects for the number of pileup tracks in an event. +/// +/// Properties: +/// VertexContainer - name of the vertex container +/// TrackParticleContainer - name of the track container +/// AssociatedTracks - name for attribute holding the list of associated tracks +/// TrackVertexAssociation - name for the container holding the track-vertex associations +/// TrackSelector - tool to select tracks (none ==> no selection) +/// JVTFileName - ROOT Filename containing JVT likelihood histogram +/// JVTLikelihoodHistName - JVT Likelihood histogram name +/// JVTName - name for the 3 JVT attributes (default is "JVT") +/// K_JVFCorrScale - the scale factor for pileup tracks in the JVFCorr calculation (default is 0.01) +/// Z0Cut - Z0 value (in mm) within which tracks associated to no vertex may be assigned to the primary vertex (default is 3.) +/// PUTrkPtCut - the track pT (in MeV) below which tracks associated to no vertex may be assigned to the primary vertex (default is 30000.) + +#include "AsgTools/ToolHandle.h" + +#include "xAODTracking/Vertex.h" +#include "xAODTracking/VertexContainer.h" +#include "xAODTracking/TrackParticle.h" +#include "xAODTracking/TrackParticleContainer.h" + +#include "JetInterface/IJetTrackSelector.h" +#include "JetRec/JetModifierBase.h" +#include "JetEDM/TrackVertexAssociation.h" + +#include "AsgTools/AsgTool.h" +#include "JetInterface/IJetModifier.h" + +#include <vector> +#include <string> + + +#include <TFile.h> +#include <TString.h> +#include <TH1D.h> +#include <TH2D.h> + +//class JetVertexTaggerTool : public JetModifierBase +class JetVertexTaggerTool : public asg::AsgTool,virtual public IJetModifier +{ + ASG_TOOL_CLASS(JetVertexTaggerTool,IJetModifier); + + public: + // Constructor from tool name + JetVertexTaggerTool(const std::string& name); + + // Initialization. + StatusCode initialize(); + + // Inherited methods to modify a jet + virtual int modify(xAOD::JetContainer& jetCont) const; + + // Finalization. + StatusCode finalize(); + + // Local method to return the primary and pileup track pT sums + // this method also allows the standard jvf to be calculated + std::pair<float,float> getJetVertexTrackSums(const xAOD::Vertex*, + const std::vector<const xAOD::TrackParticle*>&, + const jet::TrackVertexAssociation*) const; + + // Local method to count the number of pileup tracks in the event + int getPileupTrackCount(const xAOD::Vertex*, + const xAOD::TrackParticleContainer*&, + const jet::TrackVertexAssociation*) const; + + + private: // data + + // Configurable parameters + std::string m_verticesName; + std::string m_assocTracksName; + std::string m_tvaName; + std::string m_tracksName; + std::string m_jvtlikelihoodHistName; + std::string m_jvtfileName; + std::string m_jvtName; + ToolHandle<IJetTrackSelector> m_htsel; + TString fn; + TFile * m_jvtfile; + TH2F * m_jvthisto; + float m_kcorrJVF; + float m_z0cut; + float m_PUtrkptcut; +}; + + +#endif + diff --git a/Reconstruction/Jet/JetMomentTools/Root/JetCaloEnergies.cxx b/Reconstruction/Jet/JetMomentTools/Root/JetCaloEnergies.cxx new file mode 100644 index 00000000000..49112cededf --- /dev/null +++ b/Reconstruction/Jet/JetMomentTools/Root/JetCaloEnergies.cxx @@ -0,0 +1,117 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#include "JetMomentTools/JetCaloEnergies.h" + +#include "xAODJet/JetAccessorMap.h" + +#include "JetUtils/JetCaloQualityUtils.h" + +#include "xAODCaloEvent/CaloCluster.h" +#include "xAODPFlow/PFO.h" + +#include <vector> + +//********************************************************************** + +JetCaloEnergies::JetCaloEnergies(const std::string& name) +: JetModifierBase(name) { } + +//********************************************************************** + +int JetCaloEnergies::modifyJet(xAOD::Jet& jet) const { + ATH_MSG_VERBOSE("Begin modifying jet."); + 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.0; // re-initialize + + if ( numConstit == 0 ) { + ATH_MSG_VERBOSE("Jet has no constituents."); + return 1; + } + + // should find a more robust solution than using 1st constit type. + xAOD::Type::ObjectType ctype = jet.rawConstituent( 0 )->type(); + if ( ctype == xAOD::Type::CaloCluster ) { + ATH_MSG_VERBOSE(" Constituents are calo clusters."); + fillEperSamplingCluster(jet, ePerSampling); + + } else if (ctype == xAOD::Type::ParticleFlow) { + ATH_MSG_VERBOSE(" Constituents are pflow objects."); + fillEperSamplingPFO(jet, ePerSampling); + + }else { + ATH_MSG_VERBOSE("Constituents are not calo clusters nor pflow objects."); + } + + return 1; +} + + +void JetCaloEnergies::fillEperSamplingCluster(xAOD::Jet& jet, std::vector<float> & ePerSampling ) const { + // loop over raw constituents + size_t numConstit = jet.numConstituents(); + for ( size_t i=0; i<numConstit; i++ ) { + const xAOD::CaloCluster* constit = dynamic_cast<const xAOD::CaloCluster*>(jet.rawConstituent(i)); + for ( size_t s= CaloSampling::PreSamplerB; s< CaloSampling::Unknown; s++ ) { + ePerSampling[s] += constit->eSample( (xAOD::CaloCluster::CaloSample) s ); + } + } + 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 ); +} + +void JetCaloEnergies::fillEperSamplingPFO(xAOD::Jet & jet, std::vector<float> & ePerSampling ) const { + float emTot=0; + float hecTot=0; + float eTot =0; + size_t numConstit = jet.numConstituents(); + for ( size_t i=0; i<numConstit; i++ ) { + const xAOD::PFO* constit = dynamic_cast<const xAOD::PFO*>(jet.rawConstituent(i)); + + float att = 0.; + if( constit->attribute<float>(xAOD::PFODetails::eflowRec_LAYERENERGY_HEC0, att) ) + ePerSampling[CaloSampling::HEC0] += att; + // approximation : if there is energy in HEC0, assume nothing in EMB3 and all in EME3 + CaloSampling::CaloSample ems = att > 0 ? CaloSampling::EME3 : CaloSampling::EMB3 ; + // similarly either everything in TileExt0 or in TileBar0 + CaloSampling::CaloSample tiles = att > 0 ? CaloSampling::TileExt0 : CaloSampling::TileBar0 ; + + if( constit->attribute<float>(xAOD::PFODetails::eflowRec_LAYERENERGY_EM3, att) ) + ePerSampling[ems] += att; + if( constit->attribute<float>(xAOD::PFODetails::eflowRec_LAYERENERGY_Tile0, att) ) + ePerSampling[tiles]+= att; + + if(constit->ptEM() != 0. ){ + eTot += att*constit->eEM(); + if( constit->attribute<float>(xAOD::PFODetails::eflowRec_EM_FRAC_ENHANCED, att) ) + emTot += att*constit->eEM(); + if( constit->attribute<float>(xAOD::PFODetails::eflowRec_LAYERENERGY_HEC, att) ) + hecTot += att*constit->eEM(); + } + } + 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); + + if(eTot != 0.0){ + emFracAcc(jet) = emTot/eTot; + hecFracAcc(jet) = hecTot/eTot; + } else { + emFracAcc(jet) = 0.; + hecFracAcc(jet) = 0.; + } +} + + +//********************************************************************** diff --git a/Reconstruction/Jet/JetMomentTools/Root/JetCaloQualityTool.cxx b/Reconstruction/Jet/JetMomentTools/Root/JetCaloQualityTool.cxx new file mode 100644 index 00000000000..93139ea2861 --- /dev/null +++ b/Reconstruction/Jet/JetMomentTools/Root/JetCaloQualityTool.cxx @@ -0,0 +1,124 @@ +/* + 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 <iostream> +#include <iomanip> +using namespace std; +using namespace jet; + + +JetCaloQualityTool::JetCaloQualityTool(const std::string& name) + : JetModifierBase(name) +{ + + + declareProperty("TimingCuts", m_timingTimeCuts); + declareProperty("Calculations", m_calculationNames); +} + + +int JetCaloQualityTool::modifyJet( xAOD::Jet& jet ) const +{ + + if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Inside process() method" << endreq; + + if(m_doFracSamplingMax==true) + { + // Special case : we want to store more than 1 value (max sampling AND its index). + // AND there's no code available for direct cell access. + // So we just use a direct function instead of a calculator for now. + + 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 ); + } + + // Do all other calculations + std::vector<double> results = m_calcProcessor->process(&jet); + + // store them in the jet + for(size_t i=0;i < m_calcProcessor->numCalculators();i++) { + const JetCaloCalculator* calc = m_calcProcessor->at(i); + ATH_MSG_DEBUG( calc->name() << " -> "<<results[i] ); + jet.setAttribute<float>( calc->name(), results[i] ); + } + + return StatusCode::SUCCESS; +} + + +StatusCode JetCaloQualityTool::initialize() { + ATH_MSG_DEBUG( "Inside initialize() method" ); + + // In this tool we're using the cluster-based calculators + // (this is different in the cell-based calculation tool). + m_calcProcessor = & m_jetCalculations; + + // Define and use a macro to compactify the addition of calo calculators + // The macro defines a JetCalculator as 'c' in case addition setting are needed. +#define ADDCALCULATOR( klass ) klass *c = new klass(); m_jetCalculations.addCalculator(c) + + // convert names passed as property into calo calculators + for( std::string & calcN : m_calculationNames){ + + if ( calcN == "LArQuality") { + ADDCALCULATOR( jet::JetCalcQuality ); + } else if ( calcN == "TileQuality") { + ATH_MSG_ERROR( "TileQuality calculated from clusters is actually identical to LArQuality. "); + ATH_MSG_ERROR( "No meaningful TileQuality calculation possible from cluster yet"); + return StatusCode::FAILURE; + // ADDCALCULATOR( jet::JetCalcQuality ); + // c->includeTile = true; + // c->includeLAr = false; + } else if ( calcN == "Timing") { + ADDCALCULATOR( jet::JetCalcTimeCells ); + } else if ( calcN == "HECQuality") { + ADDCALCULATOR( jet::JetCalcQualityHEC ); + // the calculator class has its name backward (QualityHEC) + // but really sets the attribute has "HECQuality" + } else if ( calcN == "NegativeE") { + ADDCALCULATOR( jet::JetCalcNegativeEnergy ); + } else if ( calcN == "AverageLArQF") { + ADDCALCULATOR( jet::JetCalcAverageLArQualityF ); + } else if ( calcN == "Centroid") { + ADDCALCULATOR( jet::JetCalcCentroid ); + } else if ( calcN == "N90Constituents") { + ADDCALCULATOR( jet::JetCalcnLeadingCells ); + c->setName( xAOD::JetAttributeAccessor::name( xAOD::JetAttribute::N90Constituents )); + } else if ( calcN == "BchCorrCell") { + ADDCALCULATOR( jet::JetCalcBadCellsFrac ); + } else if (calcN == "FracSamplingMax") { + m_doFracSamplingMax = true; // no calculator, as this is a special case. + } + + }// end loop over m_calculationNames + + + // Define OOT calculators. + for( double & timeCut : m_timingTimeCuts){ + + // build the moment name from the base-name and the value of the timing cut + std::stringstream s; + s << std::setprecision(0) << std::fixed << "OotFracClusters" << timeCut; + + JetCalcOutOfTimeEnergyFraction* c = new jet::JetCalcOutOfTimeEnergyFraction(); + c->setName(s.str()); + c->timecut = timeCut; + m_jetCalculations.addCalculator( c ); + } + + return StatusCode::SUCCESS; +} + + + diff --git a/Reconstruction/Jet/JetMomentTools/src/JetIsolationTool.cxx b/Reconstruction/Jet/JetMomentTools/Root/JetIsolationTool.cxx similarity index 93% rename from Reconstruction/Jet/JetMomentTools/src/JetIsolationTool.cxx rename to Reconstruction/Jet/JetMomentTools/Root/JetIsolationTool.cxx index 99f862eb51b..3ae6c990bba 100644 --- a/Reconstruction/Jet/JetMomentTools/src/JetIsolationTool.cxx +++ b/Reconstruction/Jet/JetMomentTools/Root/JetIsolationTool.cxx @@ -73,8 +73,8 @@ namespace jet { virtual ~IsolationCalculator() {} - virtual std::string baseName() {return "";} - virtual IsolationCalculator * clone(const xAOD::Jet* ){return NULL;} + virtual std::string baseName() const {return "";} + virtual IsolationCalculator * clone(const xAOD::Jet* ) const {return NULL;} virtual void copyFrom( const IsolationCalculator* o, const xAOD::Jet*){ m_attNames = o->m_attNames; m_kinematics = o->m_kinematics; @@ -165,20 +165,21 @@ namespace jet { return result; } - virtual std::string baseName() {return m_iso.name();} + virtual std::string baseName() const {return m_iso.name();} - virtual IsolationCalculator * clone(const xAOD::Jet* j){ + virtual IsolationCalculator * clone(const xAOD::Jet* j) const { IsolationCalculator* iso; - m_iso.setup(j); if( j->getInputType() == xAOD::JetInput::EMTopo) { auto* isoT = new IsolationCalculatorT<ISOCRITERIA, EMClusterKinematicGetter >(); isoT->m_iso = m_iso; + isoT->m_iso.setup(j); iso = isoT; } else { auto* isoT= new IsolationCalculatorT(); isoT->m_iso = m_iso; + isoT->m_iso.setup(j); isoT->m_kine = m_kine; iso = isoT; } @@ -204,7 +205,7 @@ namespace jet { return dr2< m_deltaRmax2; } - std::string name(){ + std::string name() const { std::ostringstream oss; oss << m_name << int(10*m_parameter); return oss.str(); } @@ -225,7 +226,8 @@ namespace jet { /// 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 } + 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, ) ; @@ -248,7 +250,6 @@ namespace jet { 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; } @@ -363,14 +364,15 @@ int JetIsolationTool::modify(xAOD::JetContainer& jets) const { ATH_MSG_DEBUG("modify : retrieved map"); // adapt the calculators to these jets (radius, input type, etc...) + // Since we're in a const method, we must create a copy of the calculators we own. // 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] ) ); + for( const IsolationCalculator * calc : m_isoCalculators ){ + IsolationCalculator * cloned = calc->clone( jets[0] ); + calculators.push_back( cloned ); } - - /// Loop over jets in this collection. + // Loop over jets in this collection. for( xAOD::Jet* jet: jets ){ jet::ParticlePosition jetPos(jet); @@ -391,12 +393,17 @@ int JetIsolationTool::modify(xAOD::JetContainer& jets) const { 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 ){ + for( IsolationCalculator * calc : calculators ){ calc->setIsolationAttributes( jet, nearbyC ); } + + } - + // clear calculators : + for( IsolationCalculator * calc : calculators ){ + delete calc; } + ATH_MSG_DEBUG("done"); return 0; } diff --git a/Reconstruction/Jet/JetMomentTools/Root/JetLArHVTool.cxx b/Reconstruction/Jet/JetMomentTools/Root/JetLArHVTool.cxx new file mode 100644 index 00000000000..59707dfe6c1 --- /dev/null +++ b/Reconstruction/Jet/JetMomentTools/Root/JetLArHVTool.cxx @@ -0,0 +1,46 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#include "JetUtils/JetDistances.h" + +#include "JetMomentTools/JetLArHVTool.h" +#include "xAODCaloEvent/CaloCluster.h" + +JetLArHVTool::JetLArHVTool(const std::string& name) + : JetModifierBase(name) +{ +} + + +StatusCode JetLArHVTool::initialize() +{ + return StatusCode::SUCCESS; +} + + +int JetLArHVTool::modifyJet( xAOD::Jet& jet ) const +{ + double energyHVaff=0; + int numCellsHVaff=0; + + // loop over raw constituents. Look for clusters + size_t num = jet.numConstituents(); + for(size_t i=0;i<num;i++) { + const xAOD::CaloCluster * cl = dynamic_cast<const xAOD::CaloCluster *>(jet.rawConstituent(i)); + if( !cl) continue; + + double m; + if(cl->retrieveMoment(xAOD::CaloCluster::ENG_BAD_HV_CELLS,m) ) energyHVaff += m; + if(cl->retrieveMoment(xAOD::CaloCluster::N_BAD_HV_CELLS, m) ) numCellsHVaff += m; + } + + double emE = jet.jetP4(xAOD::JetEMScaleMomentum).E(); + // set the attributes + jet.setAttribute<float>("LArBadHVEnergyFrac", energyHVaff / emE); + jet.setAttribute<int>("LArBadHVNCell", numCellsHVaff); + + return 0; +} + + diff --git a/Reconstruction/Jet/JetMomentTools/Root/JetOriginCorrectionTool.cxx b/Reconstruction/Jet/JetMomentTools/Root/JetOriginCorrectionTool.cxx new file mode 100644 index 00000000000..a706edb180b --- /dev/null +++ b/Reconstruction/Jet/JetMomentTools/Root/JetOriginCorrectionTool.cxx @@ -0,0 +1,49 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +// JetOriginCorrectionTool.cxx + +#include "JetMomentTools/JetOriginCorrectionTool.h" + +#include "JetUtils/JetOriginHelpers.h" + +//********************************************************************** + +JetOriginCorrectionTool::JetOriginCorrectionTool(const std::string& myname) +: asg::AsgTool(myname) { + declareProperty("VertexContainer", m_vtxContainerName="PrimaryVertices"); + declareProperty("OriginCorrectedName", m_correctionName="JetOriginConstitScaleMomentum"); +} + +//********************************************************************** + +int JetOriginCorrectionTool::modify(xAOD::JetContainer& jetCont) const { + + const xAOD::VertexContainer *vxContainer = 0; + + // retrieve the VertexContainer. if fails, fill the jets with null vector + if( (evtStore()->retrieve(vxContainer, m_vtxContainerName).isFailure()) || vxContainer->empty() ){ + ATH_MSG_WARNING("Could not find origin Vertex. Filling jet with null "<< m_correctionName); + xAOD::JetFourMom_t null; + for(xAOD::Jet * j : jetCont) j->setAttribute<xAOD::JetFourMom_t>(m_correctionName, null); + return 0; + } + + // choose PV. + const xAOD::Vertex *vx = vxContainer->at(0); + + ATH_MSG_DEBUG(" correcting jets "); + for(xAOD::Jet * jet : jetCont){ + ATH_MSG_DEBUG(" ----> jet "<< jet); + ATH_MSG_DEBUG(" jet pT: "<< jet->pt()); + xAOD::JetFourMom_t fv = jet::clusterOriginCorrection(*jet,*vx); + ATH_MSG_DEBUG(" " << m_correctionName << " pT: " << fv.pt()); + jet->setAttribute<xAOD::JetFourMom_t>(m_correctionName, fv); + jet->setAssociatedObject("OriginVertex", vx); + + } + return 0; +} + +//********************************************************************** diff --git a/Reconstruction/Jet/JetMomentTools/Root/JetVertexFractionTool.cxx b/Reconstruction/Jet/JetMomentTools/Root/JetVertexFractionTool.cxx index 72948bd303f..2336fb7cdd7 100644 --- a/Reconstruction/Jet/JetMomentTools/Root/JetVertexFractionTool.cxx +++ b/Reconstruction/Jet/JetMomentTools/Root/JetVertexFractionTool.cxx @@ -2,20 +2,40 @@ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration */ +// JetVertexFractionTool.cxx #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); +, m_tvaName(""), + m_htsel("") { + declareProperty("VertexContainer", m_verticesName); + declareProperty("AssociatedTracks", m_assocTracksName); + declareProperty("TrackVertexAssociation", m_tvaName); + declareProperty("TrackSelector", m_htsel); + declareProperty("JVFName", m_jvfname ="JVF"); +} + +//********************************************************************** + +StatusCode JetVertexFractionTool::initialize() { + ATH_MSG_INFO("Initializing JetVertexFractionTool " << name()); + if ( m_htsel.empty() ) { + ATH_MSG_INFO(" No track selector."); + } else { + ATH_MSG_INFO(" Track selector: " << m_htsel->name()); + } + ATH_MSG_INFO(" Attribute name: " << m_jvfname); + return StatusCode::SUCCESS; } +//********************************************************************** + int JetVertexFractionTool::modifyJet(xAOD::Jet& jet) const { // Get the vertices container @@ -44,72 +64,87 @@ int JetVertexFractionTool::modifyJet(xAOD::Jet& jet) const { // Get and set the JVF vector const std::vector<float> jvf = getJetVertexFraction(vertices,tracks,tva); - jet.setAttribute("JVF",jvf); + jet.setAttribute(m_jvfname, jvf); // Get and set the highest JVF vertex - jet.setAttribute("HighestJVFVtx",getMaxJetVertexFraction(vertices,jvf)); + jet.setAttribute("Highest" + m_jvfname + "Vtx",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; +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; +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); +//********************************************************************** + +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; + bool notsel = m_htsel.empty(); + unsigned int nkeep = 0; + unsigned int nskip = 0; + for ( size_t itrk = 0; itrk < tracks.size(); ++itrk ) { + const xAOD::TrackParticle* ptrk = tracks.at(itrk); + if ( notsel || m_htsel->keep(*ptrk) ) { + const xAOD::Vertex* ptvtx = tva->associatedVertex(ptrk); if ( ptvtx != nullptr ) { - if ( ptvtx->index() == vertex->index() ) sumTrackPV += track->pt(); + if ( ptvtx->index() == vertex->index() ) sumTrackPV += ptrk->pt(); } - sumTrackAll += track->pt(); + sumTrackAll += ptrk->pt(); + ++nkeep; + } else { + ++nskip; } - - return sumTrackAll>0 ? sumTrackPV/sumTrackAll : -1; + } + double jvf = sumTrackAll>0 ? sumTrackPV/sumTrackAll : -1; + ATH_MSG_VERBOSE("JetVertexFractionTool " << name() + << ": nsel=" << nkeep + << ", nrej=" << nskip + << ", " << m_jvfname << "=" << jvf); + return jvf; } -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; +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/JetVertexTaggerTool.cxx b/Reconstruction/Jet/JetMomentTools/Root/JetVertexTaggerTool.cxx new file mode 100644 index 00000000000..b0a9e0bc77a --- /dev/null +++ b/Reconstruction/Jet/JetMomentTools/Root/JetVertexTaggerTool.cxx @@ -0,0 +1,263 @@ +///////////////////////// -*- C++ -*- //////////////////////////// + +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +// JetVertexTaggerTool.cxx +// Implementation file for class JetVertexTaggerTool +// Author: James Frost <james.frost@cern.ch> +/////////////////////////////////////////////////////////////////// + +#include "JetMomentTools/JetVertexTaggerTool.h" +#include "PathResolver/PathResolver.h" + +//********************************************************************** + +JetVertexTaggerTool::JetVertexTaggerTool(const std::string& name) +: asg::AsgTool(name) +, m_verticesName("") +, m_assocTracksName("") +, m_tvaName("") +, m_tracksName("") +, m_jvtlikelihoodHistName("") +, m_jvtfileName("") +, m_htsel("") +{ + declareProperty("VertexContainer",m_verticesName); + declareProperty("AssociatedTracks",m_assocTracksName); + declareProperty("TrackVertexAssociation",m_tvaName); + declareProperty("TrackParticleContainer",m_tracksName); + + declareProperty("JVTFileName",m_jvtfileName = "JVTlikelihood_20140805.root"); + declareProperty("JVTLikelihoodHistName",m_jvtlikelihoodHistName = "JVTRootCore_kNN100trim_pt20to50_Likelihood"); + declareProperty("TrackSelector", m_htsel); + declareProperty("JVTName", m_jvtName ="Jvt"); + declareProperty("K_JVFCorrScale",m_kcorrJVF = 0.01); + declareProperty("Z0Cut",m_z0cut = 3.); + declareProperty("PUTrkPtCut",m_PUtrkptcut = 30000.); +} + +//********************************************************************** + +StatusCode JetVertexTaggerTool::initialize() { + ATH_MSG_INFO("Initializing JetVertexTaggerTool " << name()); + + if ( m_htsel.empty() ) { + ATH_MSG_INFO(" No track selector."); + } else { + ATH_MSG_INFO(" Track selector: " << m_htsel->name()); + } + + // Use the Path Resolver to find the jvt file and retrieve the likelihood histogram + fn = PathResolverFindCalibFile(m_jvtfileName); + ATH_MSG_INFO(" Reading JVT file from:\n " << m_jvtfileName << "\n"); + ATH_MSG_INFO(" resolved in :\n " << fn << "\n\n"); + + m_jvtfile = TFile::Open(fn); + if ( !m_jvtfile ) { ATH_MSG_FATAL( "Cannot open JVTLikelihoodFile: " << fn ); return StatusCode::FAILURE; } + + ATH_MSG_VERBOSE("\n Reading JVT likelihood histogram from:\n " << fn << "\n\n"); + + m_jvthisto = (TH2F*)m_jvtfile->Get(m_jvtlikelihoodHistName.c_str() ); + if ( !m_jvthisto ) + { + ATH_MSG_FATAL( "\n Found JVT file, but JVT histogram missing. Aborting..." ); + return StatusCode::FAILURE; + } + + return StatusCode::SUCCESS; +} + +//********************************************************************** + +int JetVertexTaggerTool::modify(xAOD::JetContainer& jetCont) 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 container + const xAOD::TrackParticleContainer* tracksCont = NULL; + if ( evtStore()->retrieve(tracksCont,m_tracksName).isFailure() ) { + ATH_MSG_ERROR("Could not retrieve the TrackParticleContainer from evtStore: " << m_tracksName); + return 2; + } + ATH_MSG_DEBUG("Successfully retrieved TrackParticleContainer from evtStore: " << m_tracksName); + + // 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); + + if (vertices->size() == 0 ) { + ATH_MSG_WARNING("There are no vertices in the container. Exiting"); + return 4; + } + + // Count pileup tracks - currently done for each collection + const int n_putracks = getPileupTrackCount(vertices->at(0), tracksCont, tva); + + for(xAOD::Jet * jet : jetCont) + { + + // 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."); + } + + // Retrieve the Vertex associated to the jet. + const xAOD::Vertex* jetorigin; + if( ! jet->getAssociatedObject<xAOD::Vertex>("OriginVertex", jetorigin) ) + { ATH_MSG_VERBOSE("Jet has no vertex associated. Using leading vertex" ); + jetorigin = vertices->at(0); + } + else + { + ATH_MSG_VERBOSE("Jet has associated vertex. Using this vertex." ); + } + + // Get the track pTs sums associated to the primary (first key of pair) and PU (second key) vertices. + const std::pair<float,float> tracksums = getJetVertexTrackSums(jetorigin, tracks, tva); + + // Calculate RpT and JVFCorr + // Default JVFcorr to -1 when no tracks are associated. + float jvfcorr = -999.; + if(tracksums.first + tracksums.second > 0) { + jvfcorr = tracksums.first / (tracksums.first + ( tracksums.second * m_kcorrJVF) / std::max(n_putracks, 1) ); + } + else { jvfcorr = -1; } + + jet->setAttribute(m_jvtName+"Jvfcorr",jvfcorr); + const float rpt = tracksums.first/jet->pt(); + jet->setAttribute(m_jvtName+"Rpt",rpt); + + // Look up JVT value + float jvt = -999.; + if (jvfcorr == -1.) { jvt = -0.1; } + else { + float rpt_inputtojvt = std::min(rpt, (float) 1. ); + int bin = m_jvthisto->FindBin(jvfcorr, rpt_inputtojvt); + jvt = m_jvthisto->GetBinContent(bin); + jvt = m_jvthisto->Interpolate(jvfcorr, rpt_inputtojvt); + } + + jet->setAttribute(m_jvtName,jvt); + + ATH_MSG_VERBOSE("JetVertexTaggerTool " << name() + << ": Primary trk pT=" << tracksums.first + << ", Pileup trk pT=" << tracksums.second + << ", Old JVF=" << tracksums.first/(tracksums.first+tracksums.second) ); + + ATH_MSG_VERBOSE("JetVertexTaggerTool " << name() + << ": JVT=" << jvt + << ", RpT=" << rpt + << ", JVFCorr=" << jvfcorr ); + + // Done + + } + + return 0; +} + +//********************************************************************** + +StatusCode JetVertexTaggerTool::finalize() { + ATH_MSG_INFO("Finalizing JetVertexTaggerTool " << name()); + m_jvtfile->Close(); + + return StatusCode::SUCCESS; +} + +//********************************************************************** + +std::pair<float,float> JetVertexTaggerTool::getJetVertexTrackSums(const xAOD::Vertex* vertex, + const std::vector<const xAOD::TrackParticle*>& tracks, + const jet::TrackVertexAssociation* tva) const +{ + float sumTrackPV = 0; + float sumTracknotPV = 0; + bool notsel = m_htsel.empty(); + unsigned int nkeep = 0; + unsigned int nskip = 0; + for (size_t iTrack = 0; iTrack < tracks.size(); ++iTrack) { + const xAOD::TrackParticle* track = tracks.at(iTrack); + if ( notsel || m_htsel->keep(*track) ) { + const xAOD::Vertex* ptvtx = tva->associatedVertex(track); + + // Check track provenance + if ( ptvtx == nullptr ) { + + // No track associated, check if z0 within cut. + if( (fabs(track->z0()+track->vz()-vertex->z()) < m_z0cut) ) {sumTrackPV += track->pt(); } // if pass z0 cuts, assign track without vertex to PV track sum + } + else { + // Track has vertex, assign to appropriate pT sum + if ( ptvtx->index() == vertex->index() ) { sumTrackPV += track->pt(); } + else {sumTracknotPV += track->pt(); } + } + ++nkeep; + } + else { ++nskip; } + } + ATH_MSG_VERBOSE("JetVertexTaggerTool " << name() + << ": nsel=" << nkeep + << ", nrej=" << nskip ); + + return std::make_pair(sumTrackPV,sumTracknotPV); + +} + +//********************************************************************** + +int JetVertexTaggerTool::getPileupTrackCount(const xAOD::Vertex* vertex, + const xAOD::TrackParticleContainer*& tracksCont, + const jet::TrackVertexAssociation* tva) const +{ + int n_pileuptrackcount = 0; + bool notsel = m_htsel.empty(); + unsigned int nkeep = 0; + unsigned int nskip = 0; + int tot_count = 0; + for(size_t iTrack = 0; iTrack < tracksCont->size(); ++iTrack) + { + const xAOD::TrackParticle * track = tracksCont->at(iTrack); + if ( notsel || m_htsel->keep(*track) ) { + const xAOD::Vertex* ptvtx = tva->associatedVertex(track); + // Count track as PU if associated with non-primary vertex and within pT cut. + // N.B. tracks with no vertex associated may be added to PV track sums, but not PU sums, nor the PU vertex counting. + if ( ptvtx != nullptr ) { + if ( (ptvtx->index() != vertex->index() ) && (track->pt() < m_PUtrkptcut) ) ++n_pileuptrackcount; + } + tot_count++; + ++nkeep; + } + else { ++nskip; } + } + const int n_pileuptracks = n_pileuptrackcount; + + ATH_MSG_VERBOSE("JetVertexTaggerTool " << name() + << ": nsel=" << nkeep + << ", nrej=" << nskip + << ", total " << tracksCont->size() ); + ATH_MSG_VERBOSE("JetVertexTaggerTool " << name() + << ": n_PUtracks=" << n_pileuptracks + << ", total=" << tot_count ); + + return n_pileuptracks; +} + + + + diff --git a/Reconstruction/Jet/JetMomentTools/cmt/Makefile.RootCore b/Reconstruction/Jet/JetMomentTools/cmt/Makefile.RootCore index 08a084251e3..5e9463127b9 100644 --- a/Reconstruction/Jet/JetMomentTools/cmt/Makefile.RootCore +++ b/Reconstruction/Jet/JetMomentTools/cmt/Makefile.RootCore @@ -12,7 +12,7 @@ PACKAGE_OBJFLAGS = PACKAGE_LDFLAGS = PACKAGE_BINFLAGS = PACKAGE_LIBFLAGS = -PACKAGE_DEP = xAODJet JetEDM JetInterface JetRec +PACKAGE_DEP = xAODJet xAODPFlow JetEDM JetInterface JetRec PACKAGE_TRYDEP = PACKAGE_CLEAN = PACKAGE_NOGRID = diff --git a/Reconstruction/Jet/JetMomentTools/cmt/requirements b/Reconstruction/Jet/JetMomentTools/cmt/requirements index 031621bd8ed..ef443c29414 100755 --- a/Reconstruction/Jet/JetMomentTools/cmt/requirements +++ b/Reconstruction/Jet/JetMomentTools/cmt/requirements @@ -4,6 +4,8 @@ author P-A Delsart <delsart at in2p3 fr> private use xAODCaloEvent xAODCaloEvent-* Event/xAOD +use xAODPFlow xAODPFlow-* Event/xAOD + use CaloDetDescr CaloDetDescr-* Calorimeter use CaloEvent CaloEvent-* Calorimeter use PathResolver PathResolver-* Tools @@ -12,16 +14,20 @@ use AthenaKernel AthenaKernel-* Control public use AtlasPolicy AtlasPolicy-* use AtlasROOT AtlasROOT-* External -use GaudiInterface GaudiInterface-* External -use CaloIdentifier CaloIdentifier-* Calorimeter +use GaudiInterface GaudiInterface-* External +use AsgTools AsgTools-* Control/AthToolSupport + use xAODJet xAODJet-* Event/xAOD +use xAODTracking xAODTracking-* Event/xAOD + +use CaloIdentifier CaloIdentifier-* Calorimeter +use LArElecCalib LArElecCalib-* LArCalorimeter + 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 +use JetInterface JetInterface-* Reconstruction/Jet library JetMomentTools *.cxx ../Root/*.cxx -s=components *.cxx apply_pattern component_library diff --git a/Reconstruction/Jet/JetMomentTools/python/JetMomentsConfigHelpers.py b/Reconstruction/Jet/JetMomentTools/python/JetMomentsConfigHelpers.py index 616f98b1f95..f77f65d7e1c 100644 --- a/Reconstruction/Jet/JetMomentTools/python/JetMomentsConfigHelpers.py +++ b/Reconstruction/Jet/JetMomentTools/python/JetMomentsConfigHelpers.py @@ -108,6 +108,7 @@ defaultMoments = { 'ktDeltaR' : lambda *l : getsubJetsMomentTool(), 'isolation' : lambda *l : getJetIsolationTool("JetIsolation"), 'jvf' : lambda jetcollname,*l : getJetVertexAssociationTool(toolName=jetcollname+'JVAtool') , + 'jvt' : lambda jetcollname,*l : getJetVertexAssociationTool(toolName=jetcollname+'JVAtool') , 'trackMoments' : lambda *l : getJetTracksMomentTool(), 'truthMF' : lambda *l : scheduleMFTool("Truth",*l) , 'trackMF' : lambda *l : scheduleMFTool("Track",*l) , diff --git a/Reconstruction/Jet/JetMomentTools/share/JVTlikelihood_20140805.root b/Reconstruction/Jet/JetMomentTools/share/JVTlikelihood_20140805.root new file mode 100644 index 0000000000000000000000000000000000000000..f1bcd28953a479ee3f46f81f779d187f42f00329 GIT binary patch literal 19745 zcmbTdXH*m27d|R4BCjGuL5)g@iW(6G6@kza6{QMMP>~u10g)PdAOR5-5dmq^0#N~x z8hQ<(6Y0`hAdmo|Luer+ft&Ytzx+SkweDRvYt1@)?em;x?=v~G=ggkj?yjz$2M%=q zJ#gT_s{;oN&L24Nu#P9c<Q;sxBjVuyEvV@oIB@Lm0Y1Yqr5i%#$EV|g^{nu>#8qDJ z|EFKuf&Yj+49Gj@4w(FpnP=_50f8r44o<d;j)sOeZ`|;7cW_4NJ$qurb4k<H-4@~I zdFzI!tKtpBLkCA&Cx_RruCEX}4j!Jaa_;|m{r_{=0fGO#_@5^S4*cRd5y3nD=M8V} zXaC=(s4DT_O#bBkZ@m70GvSl_Z~Vb0I=8g0{-56eACvljn8*RMhyVUjIgt0~fG^*2 z@?rAegMmHQoAI$$?z+VfjoWWOeEjpzd(9*>yPtN^wl3)oNfq>@S1*%4WF^KrgeLx? z-n%My@!h8*SJ*^#Hu39&3{!5mPpSv^tm|rs!CdD)ox2rbKiiSGTs4bLtVB)R&9Gm+ z>a6~|WRcoy&Zz~&DM>!Lx+U*Cy}#i~GjiAfNO>>a)CP!q*WYdIQrl`$h=f>EpcpyG z>K8oDw@ZM%Rr%t1Efg#j(2|&$>eVve_z{ya?W`>C+)sNr>K|b2H&QX;VIJ#O;cYH% zweH&OVd-HwpSMTxnUQvOXS4!w$UQsFkQASZmn11_>8fl^#*@*xnO7V2k&&sj&2Jks z4PY;3nxecrOXjC*W`ifDvfx&34qHEL%4#~EIVIo9_a9gO{p9So7V?y$ja5TJvJ{4D zOPBEYkn9PYesAO;`H`b_);-?oPw(aUThsfH+*Pi!=h-S!dj5cv+f-9j^BaVkr}eLC zT;62uI7VDqIH_NWP(5k^n*x(Ee@@pfD!DwFh){RlSr&Kpt4cYszWG>euc(LVTITT+ zUBAXOq|6y>EK&a=pxe_a%1OJ&#ih1PQUa$dTruBq81Pa)(eA-kxnW=r5x5Pjpn=WN zi$3<liajeNvF&AcURn6^=<EHM;-nU|bv1CQYL33E-(rAXNK+mY?xDc(dlMWuddjo< z<}P<+3-KDW5FI&e(6j#=v242A)a^gqZG%OSH`dGG8cK~W5sU!VoQHpaHKWs8|BL^k zLfRqqt%ac_lbd_^q`<%-x%~|xwpRyjNY;Dvo5oUELbr1ut7`|b>QNOQfeH^Ob=E7b zd0kH81QauF7JI7{;kcJ)Q6zKRbN6*!<uS>d6bg}x5w3-nhj5i<F?8=I*q1&noYhfA zcOYd3EI(>l7(rRyo-r^aE=zA$Kkpv*H(c3T6m%`Cwy|Vxvi^z;#3+^dTVt33)0wMd zw%cdCU2{D;`1Xg{)1zmUE)V)s?q+W9Z{c9p%;u_9c%MqcO5_$LkY0^W&so~5n+lrH zCUffs0VDiDBz~Lxs=(Rmq@|#`nR*N;U@o${W=!fRT8e{dr<rd_IROqC+OSurr4hC^ z{P@C1kOu2_B<z(9e>nv}*kySD-V{*;18M^o3W0&{I}zKDCU^dV)n?V*I$ro+Swqi_ zSOUA$_pbUw+%P)|Bsn(CM5-u+pAv%DxjX*MUU5u1XhT|?`vXN0@+r~&XH|Qa(2wjt z8SHE^30)DY-k&EQHyNp@3nWwJ_m;shzhNn1DJhL9(!TFxX6N?nHOA|$bSVYbksc`_ zD=BsxLYYmI3|=mrq!<30a2<QoQNGNWoM$M&$$pd&P#Gy0Pg2XLOk}=-A9RId)oiu* zVg46({4ewbriTNhNC-mr3ZP@ZEF1tO?PFfACbomhU<dJmJG0%poWjD9y<;5cd{@Q5 zXXbEB{YDIhHAG@^$D+1ht!YeMSKHAODLGuvejiMl*q$>dz7l)4f0)8ua)s8gJ}8~v zu77VY9#~yKk95Ab6q30%QBHAyHxqpq)_mA)mq(2YVX3=4noP`oO4xGNFc@9~*{PTG zw%}m`5_h~+7i`6BC=b<l$0E^7T>s&wVVIb!ZzO*6PV#B?Ybseg6|$m<C<-maT@IQc z0!Pz_O=8^k1EigemhieT5M$FE=ACKq8KdOT+C^LsRg&2r>+{?Xy|nGQRGkUuko&a( zRt|8|p8u}NvOIdXz>3$=nCVDVvza>XKRi7g8N<j=1C+ED%7bycJLLFCjfu;oJuGLe zTYhH`D9`JFgO$DHCZ^s;Z9iNZvWM<)Gkj9q<2A(C+lj^AIlIFSV#|rcwcTVj;%p`U zdZE2F2-c^PzE@FfG8#MFZTy1Gu~NFUy%oC=Fznu{yj>zwSewb&<1V|_EZs>KWux2L zz-ro-D1wGKn>!LI#q3R`5dWlxWA>Mt19vylD3T7{Eqj66g&^)$AB@nRDh%`89sm#s zZqxVu<v3F+k>BhsYMWApYJFGQ<hdI&IBs*Q(7vynB(P?waxI4LZ_yr6Om9S_?oSIU z*DP)aFR|8=>0*|ov5hHMH(cdqJAI|{PbvVGRK=p~(JR>6qy1}{Ss|7rk)o_A6Jz{{ zV=^!1@MMH>CagF`)7iMO``4E0d<yY)@~J(g<v%79HT3E7v68iFhoKB%`@V+#U!5cD z8eD3?R|6p!b}!<7x&a^%3!RQ-;%+i8BufRTa4#RI8vGR!Xo&EXK03C!)(7-{$Ui1D zg_lROSeTV;#2fkd6r$Ur$u|4>24d;M8U>iW$C{$6PPoW_1l3k78r-mx{|Kv1=l~P} z8}SHEX9>R?YNmbLdWRVno;hN^kBaL27ltcz9vxPa=3umHyH|E2!eOoo4XgGNY&rYc z6UoQ-V6>ktOs-vVv9nfF&>5J%c$?jQ^@RF);tbPCt4RP$`&`cLx4Vx{;3L$y6ek`B zb8^p}8}0Z6RDkXMwZc(mk;!L}jW(>109J&(vt3hUqRG*T>Hv17k^=P;$xN`I#ZF5- zF_c+b-Djs|!FXOl*<o_V?X=K|Pj|ONlu@g*O0v|WA%*y1y9b0suStJwtEKa^rWDi6 z613<13VS&B6Rg+>vRl#l8FG7Zdv$4l3lPv|_dq?VdA@Xr&Vy!{3ahCUXYzpXVr6GZ zUpL%KB{gdQ65vZQXxS&RGs0JqZ7hALle@-r@@~5RF0}iPoi-<l8=R6dKgQD(j?_A` zr)$_0`(o+T&dX+!YwtTFDwL|17c>Qrq6h7?(a9TS-8bKLhC<@YHi8tDz|?W9^eAj2 z0`5JE{p9cQQNw#D9)KE{z0XPZaK11xdtYL+G5!<}H-(kE<80W}oPVP@vNL=K;v;_9 z?`>x|WZqGuCG(Ysv{e@uz3b%C)L98!_FIw<!CHGr2Xxg$y;bg*)sk3#$8t{=XL!o> zY_(I$E03FK>;yP)LSdz-;|O|T3Y^$E`v99#zW-Sas<sJD1~R;*1~Z{AH!iGnRdWQv z*SQNGO7hfIb=()Viq0^I`YK>me0NL1(aEZbh9WBwCmc2|$aS-aK!#jWfr)cK6P>ln zoiYJHSM!u$)Z6g`Aj>t8WDvtI;0UOCl=Eq91MfHA!L%Smc1G@K1a8LIr})kD@t`}L zOgwA1CO&X<R$Cp{?XcEq!(-bLNL~r(dOMXl{O;4LP_LI=!u3y(+R49%=?42<GiB$f z+<)vmwwTU$o45z*OqUXY<NH9)3O|rLD_B1A4kCC;eV^`Hg7Cmy|FrRfkF*?#D-b}+ z?Tk<M%I0IbPPMQK)-(}SPU~J9PA{=yGV0_<HbOMQ^SR=EPIyY$jV@uBc|K>x06*FQ zft{Lwx%+bEC*BCrB&Ld(KM*eVLNs;Lxy+#t<tuh!R14D*FL!nG6ZH0<#F+evSDSog zSGp1(kkW;n!>9yP&#CxYW*<e$F^o!{KQg50{1^)hliFJZPy(4Zc%%OQl>5i&F;+W_ z<2|P=VJ<<B2CW3#CG?yku}XaRCx)cdy-Th~cwDXD{);JtbmcThM65dSV1v};KE79B zG>%!vHX5DEgE%pYYFr3D=)CI@8%~9|GFhC(OMR+^{>#ZM6K7qltNtNN_;kS~WEYD@ z15U(U58ZIlJyKph5)EN~BX74|M+6Wfc}myTbeZzqC~JMHgoTk0e^lMKCR9a#N)F|V zjk5zbOy7s%zj?Cw<7h7g6fd>T90PZ&kiNC3RBq;UPcX_*@ngLZ`FPdc0L`Kdq|2`D zTo*#Fs^?n^t1^~`Y<+EhmOj_}&17}&ZxNeEm)Xwoc6&MCYkrFU=3VBj3a?$f<+E+D zFgsENv5!v~r|)$u^^X=oT;q``u4XUWsbouU$_7BrnH6$9Y{SQND9FB@DzW0uTy{cX zN$p&vRe+Sa6n%Z2yT|$i+v_KNBe4kTmSDl*bx3AB!c8poRXdGfwX)>bQ8a^e@sUJm zI-{_ZcH=v|S&&qS9DH9Yfr&7A{uE1pssX@P<=!62rdr^ssH&bPIiELbKGR=Bnm#*8 zN5=_<;9UrH*;E4F{YDwO{Z`IrkJ`_36}mM~n^_32?cu6!=g~J$k$~S->y|4j@}^J8 zg4u4}o%?oAk@CN*6b<^bsR54!M$KKUpCaMElh)YyfO7h`XD68-vOUiak7iRR9yR0n z_csi4KC?BSR!ywDdWwwvy@G5SDhp(@^Hf*F2Fx>qhqGx`k88%N$BVa%7b*}>HXIw2 zxV+9Y1A>YvTVu|GW?p3F>^`J-!BbLATGDIcxP7kb3TOaG9Mf2SofF}&v-K(GBW5!N zQLs=E`D8yfn0}w9)!mW|?yg+7%;XJ+XZLyRJme|$y!<9o><5sDtd%1>&m@`O-neY2 zSvm|L1GDzVCSl>@sP_X+8o%r|LgY`cPk_t}$o6l%&;7NsS485o&P%NL#SAoYf3?WR z2rhf)@Q|6a1=25j-%yFh8_g#sb1$!GC}N@2u-A&n{jAERu|Qgvkr?y)*;Jsvx`QGT zlwGMJa4IL%|8b{Aq{Z&EU?$H%<QvW{r4r^fp4uoJ1WjHpt>f`MiIMQD9B8V^`Z%;R zvGq)hnR%A%fP0n`?*GIYpF;LP6pdEE2Uf5Wt#4>nPpjdJaqUcY{6Gsg+xu_NIr~y% zceVjc+u5NMN#O}5&d#OC)$A=&t(a3OazEQ(%-)_y^?cduFeE3O{roz9hRftJ#t%f; z-ea3Ole^#aG5RQPXsOSV<zCyBA`v;w)rX{bV_kC$@c70&r4?7N0MxSD-0Pe$e~WRo z-sV3GRH#XOst?#qvin5s&;VxIqoaR;8ern{$Hx5O?zHJMAd{d<Q4cs@SXHq)FcltR z;pBx}&7pb1HuZAC*a*w2{P2RR3c<mQHLvi2W=uf@f+SMiwLrCcPNFdslua$kfW_3s zxu^a+c|5MBs(k+81uF45sasQ<Y~wX=pMc4YKvW-sz1~6w=2F)E&DqaNb}O_8X%X5Q z4qHfEE`7^w^Holyzb%(9lHx|C|1!mp<rh&C5(@1Lv{ch><UUK_8Hm~GN8k?6@PsZu z7irCfJm?B;3l?C_Uab9Vw}k}c`Ow#MVh2c=l2n8;ai?e?6L(3j+%I(qu^T@?;+7B> z@7iu5;dyvJD&5#{j%OK>#{{*jzu;Lj7w3=$tia1Jgpq;H!&e3Vas1@(q**0ZOF2u0 zj@3EZ&shx~V~&0d_O5P@yyf_b?W{Q}0xJGkL039zdG~$U$BI>@V+KUO0>!49^7C9Q zszR}e%M;MPC5lZN74`DLl$}u|k7DXjz_QxLeI7lsQ9-k!&`gsYDc<{mb&aQS*Yq<? zwHHPrTy%zV^mEUzbcS>Myu0ECPc!K)KEJ}U6<ai_kdZLXwmsd_CLD=#*#w0_1NOXl zfQ1(L{#w@z@Jmrk1wwg!qvb=nnGLyEymni_rG(eNto3t1d8M$W!bjPCta|6J;}^CM zO}-*Q<4rgc;7VX$x`0(|!4RSlX>$H8XK7Xwe;8=iZ5|3`GxXulSz7AECx48`4<4o_ zYR-Gc-SM-(zgVryQP+ZB6}8F?kNNV^T$vSy0f*x=pg(wm$+~P6LkSp3qFPa!kc^Z- zQ)BXKvf=Je%3I~qUlnxalo*j0j-T1QGy~Dq*+4ywGeL>XFIIbQOx{oCkX};m4W3wi z`%6&ehIjcdU%+h^o*@#Y*j)WRYca^tfB}7-%|-5yf&x2|s40K@w>F{`oBe-`q}i`J z0t#n+m2b208+*ZcMH2gmN{@euB1!G%2uOSk5r|bJ`Trb2)UGfFJGw>o@?RTv*Y@SA zOaMDG4J0ZU;r6B;1cjv9D5bAFV4BJY-rt}yA=w#_><uWJ^$omDK@lk5@jNl)_C`0L zwxT78rqPStV@>)wM)>G@_|<(Wn)Suq_TWso=kTm%u;OT3ElfyhlYd5Y+D2k+70HY` z-rLZ9saDGIIwQbgisd$zHYmE9@GWv*08sqRx4UG=k-J(Bv$@X5clg1kJml#pB|DKP z1;=lWnNd^v+?z|ugxpILYS%dgAG*7vF3ecC^i(Z()x#vXw27T*+tC-oauwlmem_N? zS%I1A6GF_e*Y=x-%&3<Ay<T@qzxjen)9Ks4c#!nL-adIJ??AX?IMXOcd%qeq++s!> z?@wyj7a2UenmLH+6Zm6BlN|7yJRBfR88@S83?%KbUdTJn*!u;Vl>WjsvdtU`c(ZP} zU|DI4T2ASrs!bka;>}i{^+H<6b4Nn88Ml<G{iC#noukNA-%@c5_j&H<moI6b$B=R~ zfI+aY@)&ZC&X6A0rA$*tLbsWir9H3h8|AEf%)Q^pyFTIQ98gu^_*@*O&#un9&Lyal zN2%PwL6}NGX|r1XvT^s?!O}9{?$YA{ml3YcjC{upkWG!OV;IwN?2=Oj!ql9wF}j1I zji}{bpYV4MILi`i8>Nm9&Ej*|sUSyZ#-1Y$lCt}oM^@mxAltQekT>*&`A~50DAjVf z1r}IyqjQvs;YrFhkg-u(%J7QeLhW~Re*$o<QmQMrXxmr0%m$iMKT6{cL)S*|SXg!E zC{1#tt2@YJ@6|FR)rpQ$XWBYO`q;7yYL&5&TAq^V1n*8U$kUieTy|hHOSXMo%h<39 zj9rh?pww!YccAd@hY23xNNT`D0KSx@R#LMlc+WAC`MP>vwbpkeBv)f1V_62vEZz3S zdHb$05Sw|@6IqgDZiV8nBpFuAIxw)LwWRtYcYJiJc6&RUo9orYE~QOE33a*B6NQrV ze(qSWCjYX=K<8cmM%}^VtD_NNRU$hPz7+5DKtpzwSCd+~G<V}Ki`O+S0yE(1xk@rt zmUV`fr+6Al#KHXmUOO4F$^SQ!$Ji*z;fx|{ywY*?WyK-^5pxo@DGcb!CI}ajzZ2)x zqE;D0r-4=rx_sq*L=CN~hE!?e<oSKo8SMnG7Ix)^f-=&`)Si&0#-e|S6U$YY=u6pK zObEHZ_{GPIKI8Ove(=odqG1H7wzI6u*WPD0Qh&vor&>0I<z*6$Y*Xg#X+H*AyL|#( zc_spVh&Ps##|Wb?DGZU-CBI2?djd(Va)rBba&gHoJ~L1m#7k8+^<KWjm3h$1qC9Zg zoDXFAbVPx1UJC+1?C#-xn*;leq2W-Ag`2z-<>S%eeg79(q*`;-$o}u2i-D<3SoScd zAUxUmK`qM$zSkAX$OM2#KyVb(JEX>0jC=UKB3uH+k&*?EOovLt&-gq^F5RAC7nH&! zSQ$v|(r<bk<b$&rJHhXyGO^{M@Wr3Tn1sg}2Ao=Mb@C@=%Nbr~)X#rN^w5%0U|8?B zRJy#P6_YU}v1=}vd@Rtq$F{kIU@!vGORox7Qp*$~$hAM+KYvikKXuHO6}XQnKJRv5 zj08(xe~{RzxegQ7+tPt^rXK1|V^#rcFp;x2)^{`nf4e9TRlB^v)0h&zlYeM}RTw?B zFS{yLd=vKBF?+Ad$Nk6JjS&6LMST^oG{^R|ONBT0qBT%Xy1=TG-b86{shwaxy2z%= zQHdFwD=_O!B!7t0&2u9!gKOHR&)}dvGd8Rphle=hY?2{OL_Y3N^jtCn-O0~|ryYt| z9`kw8=#VMAnE|jUA|Z0=CG%z6Q~XY0eDq)HGerlx#xW|u;prJ{SW|d6pX85omvF9B zewSY;D^U$4S#BgwUHCm<a8>mNyn>x%B)NMu<{kYFeU$xPe!+M@Ump_`6&dCBEANKZ zyOV^Y0M1`wWNz$a$bAjfU2qR>cDS8=?uWLvWEu5a+4LFuH>5U@+~Fk22M!gd|B--? zjoWq)d0a?I&g}!IPCk%4u}6cXHD-S0NAFDwLQaJO)7P!yZ-n5HUDAFPcMTV*zz*Y- z$T-g)=*wxyC38=zAhD+QJ*Xn8uGjepxm{%qbz-LhKbFb83wr~vEetdYbc4f&W10X2 zvYUg~`s1C%Id9^(rm*gtjCZxZfrE{U&vpq{FBX``z%6pao!74;%u(X=GYA=xDoZNh zj0{qBl(lwCrWuaV{2cz+5QT818mku@y1fci`xXi>2qJyWlDbxdd#?;VuretesUkSx zIwidsCX5#L0Uei@CS0Ep!Yn;@;a9$BJ+kBB0)p6fl<{Ly`4L%rvCc|%<{On?MG?}L z4fdv2NzRzK53Um+)waVn1(M@W4J(XABt@;7v>PaNT-WslG&jpC;IhAymfB$Mi1E6j z?>o>C*s4wvU1QK7jy;+DQLSw4oK^jqrKo<?zQEqP&_%zCO1+H$D9z(x0Kc-VR<|-P z%(u{4KCettqpNgRDyTLT8h=gY$>Nr=g>h!q($1$4tip^KkhnLk?DkVVfnC+phuhH$ znBHwL`P0~U+jX&QJtj6pY+H7RCUiXKh@{07tY7x=%8ak{Nc;NX>0?YFgD=&FP!pNr zdki$GpL$<>LB)0LjJ0Vu#KwCqp)8C9nH7MMPZnJ~I;8c^)<yn{ha`$|9^;Y;ofaQ% z3yl8aPOT^DZ_Q!Lsp|Itz+P#|kZHu)jKw8fRK%?WwNvxgRMgS$p((ZLLbZC%h0U%6 zN-S9pc>w%k#pa^M9+KXsE5kZr#JQ(EHnpPPq_|?Ip(Tam9187rT4k0@NPCw-SKuN> zpHXDsg~R1uQHLFI<)cSj+uDG)mg=5Yy^4<udrtSYk%sB{K@Mn=jBTAneI~TC)1(M! z*k#<0@C`nMjG-`+{2TU{Svr>9Pg$J&r0R;h5ikD|?D9GUR2LL!(j?dRqWd)T*>T)o zfDd8nm_Mh51t-@DIb<TWlOt`})v)7Z>a#VEYx(*(?S<c<FBd(2=Axq`On!{zwS&@u zu0<rnfHVkrU(b(pyLIy742Ii^M=wol?i?LHrmRy##-T-I0$4jA(nSnEs-^A-Aj5pM z64iJ501`;=w>6th=dDGDFL)t7%A)J&O^LV_d#K@O2e(GaQRYXMR6(#MJP`s6-NG-% z-`Ly=ir+m(d@!;Yr`>y5*HshB;$vL_Zbcy0#hr8_0@s1v?XIuQo0LgMlTY9!Vt$I` z6Yx!dVVaO8S$j=<=UmFXZ#1FyJuKm!hWguD0PbRdrX=~{W&k$CQ}B>KBuE@CAZT71 zp@R(bPn55kJZMFRM|U4}uZeUdseq^1zW(S|F-YEGd`0Vx>uwJO%O>YgTe!9O)T-Zb zBin@<gU^ri#e5DgUTJ(v6^8zu5|1#A(qiwe%y`Y!135DOlvd@*KTJpJ+$&PtZ8Ddi zat<JGL;k41ulo)Fi3ZgLGvlp!$<~!d>=tpj=`MgN0Qq5v8s|rM9o`R^*t<gSCY2PY z04U0##maNa0;v(2RK9te&w}X%-5b53#XqkZ5+a*^Fu?($0HF6Wx9RZ<D=f<@#9&IH zOZzHKoM`2%LwNM;?23hN^%M9~15685!%Pw7d<8?Jz%KPdJAxCfu710zBteK887_(- z0(vn=<^Hmhc#S#Q4zPcn`GWqCfL)pCQ%{wo6N9B9)Pu%G(TCfQ!ls_D(Gv7Bxd!_J zf)m?{wxp5*ABv$fA#xhPi$bpPzDTn18I{guNd!x)3XckuDN_N33cALp)8f$Yh6P}a zDdM8G9Dx8dxdE|urb`mUkI=a(q~_UlYX`X`2CUv8(sTb3!K<`i3VV;3DK+C@KmTMR zhHVC8p2Q{pAO#d`K$^NQpaml)g;Kz7X34HzWP;rGj-Ny0>FBmFI7`r33_ht=<hTl! zs-v)FJl$ruvMlFf04(TG%{5jiEq+45Z4Pr;>XWIo1&qx22!ourGw|GRa9snfUqm;D z7%*jVw*y{kq~usl0lSc@<@J&jKL0ME6-)lEN)4pj7xrQGM^e+6y7reV@U?qq<ksSU zbX_9yYnL#lHNDv<)gfgMCJJ>DcFjsDGPQQD!Ixx4Z{B$W7k3H{?mvXpbxGbVsc#|4 ztYykq+@Z5hRR59v;X&|SY&1|co)nM_-~+|hKwrYhs2{=-p@3D{fhlnX4W*t)pZG53 z8&~+X;F#1dVWb|BZ_teoHdk+<evY4>o@q$<8WR7b{xChfTQNm#Co&V8@J}E@vmAxc ztnX}-!-@fjgG0Kj;Fl2x%|*TCnjcdITy$o(hxOnS+zX6*Fxi^i+nvZw4J^9Ft9ey! zLWx1Wd5ID2!6)dmtKzPSRiIT<O{!K^m1V3{Ks$!2?d9&{Wkeh%Vt(bp8PG^0@;wO6 zPdN_&?=&2Xd<qR|lEg}^Zq3#09i&I8N8ch~R|7&;k8a&m4^4KdT&q0cm8JpzM{EdK zal{8NyL}5qx;7q{?&B1;D^vqOKCt5ckdBP-{-fit%5+DZQ#rH%YragsoU%v+#ZV?? zCq+OL&=H*w{4vs`pq2S}`Q)Mf=TtiPL~YCzcudWEM;xV#GoMn}Ibv~ieR!G!7FiSI zy#7@FK^}5uY44CAG=LgHs1Q@8L%MfF9A2MMj?g@w_AE;N#tFL@hiCOBlwKZ@*X#eG zy}SFtp<WpCNM6QVcWtb5(P0sDF+w72T2@xl$4$DPkFnrZz8{2>j+KmT0;vCB-ZUD& z9x&J_91R7lX0FH!047|k_-V3ZYF#@d!-%YI>6d^#F8Pq5NQBzZ*ERl&!vJr?k5w8d zQ55?9LBTa-7E`loSf$1FFYiiLPBW<w?=D{BE&6R_P<gol>LU8z!1ar+xGrmw`~a$E z?V(o?ss}XlI&H52UNt9Q3L_V`z;)dC<bkT7biy_0nlR3KJ?dsx9y4xl{Jpv#$NmAq z+;3+u1Q;Ee0;z|ixveI<A+<XaqJlb{SHFh<NEbtI(0Tm9S$jW!Q54Rm{=>rLXTJsY znvcnGmJ-wS*2OLX*O=|{)zbILczIHF9RO|WkN<SnpS_FGkXd2x8j*i2IqIMmeAb=^ z8=~4uOcDv6n|$DY6}eC#n%_{X0#F~bKDWAAEWoLk3XSyTT%GtJ@_x$`UG8HsWurh! zRHuKAOk8;ZEq}lNIDexsOl|N56VQ9CoKAprEl(Y4J;M*%`6qyH8|V+SSA)$>1*sFh zPOIgMDf6TAWLT<U9`DsHu7qz6Yn}xbI@rd;>a<5mH4<t!jE&9^>W^nlFJ9TM=q|5- z2NgTD+Jtm$tm(EnPw{aW!qx*<=mYem>DRcv+pB4#V<o%Ku`<(`*47D@o&==3*%Zzh z?Lm$=-dtWxXlSlqfpI0;p)T<C$_&ph9$C{HdZz(i<odI6I$+<xH7$<2x)^BffPH^j zCN|@pID+m7`LTpK(K;0G@iXCx+dlweZ7=kYxtnwx?co^<4Adbo803$grd^U;o_v4H zdpM>agc|sfb`<EQ(8Jb_@Cul0XCAlL(aW0vqUaYGG<`ir+ubWG)@yv)Y_(@Kk_zik zt$Jb3B%n!t)9tP6*FTLx7s+pHAm7o2?1v~(5b@U&L3rv_cbXD@>OSG2hlKj}DdZ73 zZy^{%@uF1SZBk6pT8z|TV3Sm{Wcucd@W39hv+l6To5+YVw}<b9D$Y-7Lg}6)vw8J$ z)j{aKC(MdUq{(p0iqhJl5o6ygT%$mf%oo0MbYB@>eB$2>*7!IjyEf@@#rgypo6>v* zp`RTBEC><2$T_D=EV~#iHV{XZMYlK`lh?gY*%#v@6*4Y>RX$Ri**cKf=X!!G)AbvN zV1n;uZqW}<6lB2Dm;SQNjLg@V7~{2b&^pWm>^za^If3^)s)T%4?xn7)p{}*%IRSxR zq473c?CXsF?qz?qS)EN%Mdu4T|57W!U43d{^ONq{`TaZEQqWj;^3;d-B5KnoiNLRO zt^P;e=;5nt8?@71FalPe9*3}r9|doK55ryx=PNg<7^Lp-|FY3o`GTaJ3tGG$fV=b< z$A0HI{775PdI(K-{)L;Lj&f^|a<}@d%-Ke!O4nn+q*I~zitbFxxz@%YJU9b75#UNX z8`@pRKM!#QM?^4{LRe!n2Uh`dlXp;@e><jx2$C_ef1gf%pZw^zTR)+4wc-q&t|mjf zI9`m#UkYSwrEeUZcf3jkvXgF;nP6SsBD_+5Dd@Vzx%wfJuKk8@q{_7Jq=BjW`DTzG zvS-2n07oN(at;x?(MVbC&7jPeyOI4^P$uh?1Auf)N@qt%lrlWt{kmdqarmg4+?pU- zptX|9mS>#|&@{!?@zK@DN4-*88wF3NX-3zHABO9YWTqma-C^;QMJFg{8UKKZA3aai z3GQ96T0806TFKj0KYRytOQo*5i2Bk|jQ8nX%+M~UDXXSTsfZ!F2(cpjp+qAc$jBf~ zgz#1^E0T6Vxf+s5jG#Kh&jkDqiPx-wBaBr|A~v>zB}n|3qdRZKds`<C)z-r&`L8v{ zpbth^?`J~1kWDvr7g`ZhHVr^jMcPqB>Si2UzY54c8&*uj0&$4XHfuL;S%2f)Q?F|t z{seP7FxY@5dYCU`bmZn$lRe~&Fv09#wiflJ5J^x9ZFSgZ&3?a;Pj%r+8fAjfEU)8M zRBpp3xOR_NYr{AD=WfIU7zub-c}IwMo^ki3x+WJ-uYz1A1q|z0o}ge1pv&#v85J`7 zpwIg7VgK}<w_=w(`%zabkl0>v_Y54c>1M?m3fMXn#kr521S{Vsze|2BRo&EF@yS>R zk~nW41bMxckbZFGqHhp*UAC%g&2z}hP*i@B7-z#jvGyHP8@7sQn}lm%zk-*ieLqW& z)&t0!9$(OwGFtVGQ(a(&KBq2p)9vm%ame|45@TZj5&bzKDnzItL{{{<vDW91FpP!K z3q=3pmI&36YP7KO8r9K&uq0X7KEya{hz6w(3&@u@soWYxb#VEiIveB4-ykbNDhx2! zJQk2-Ev8;4M)T973-&TUu2i)G>p)IEJ{A|J1Vp~N9Z-&)3c96>Z4mP&{yI@TY5kTK z^j@6NH^fx;+qwJS(;sJQ`qipnM4cok3TH95!Mm&3z)f{+-AG!yUqTvcZ>3bobg{>a z6dn+%n3Jx@(Fd(2X@-%e={LJ*zA)`ys18xkz>46ch&+A=tF)o9pu6|oLp=rPcDDvH zKVR|ako}GE#-Y(M=7mso36^@88j$c7rn)c{?X5Sr4&7C8!APXtm2#b=UvUlm9a&!6 z9kYJ3OO|j&ShyC_M-nt|1N-l<y^T1i-G9C(_$CT{*LJ>1SQOV9v&wCXSv=x`wrn2Z zErw<G<&+J#9|ytfacaRag~-}>@e7k_1M~&S+Cx5FN7lZb0Ik*S+2LYhNrH97cs$@* z6cm3UNKK+XloPr;0WI_W{xH?ARA3(~hOnm|3TSGdpO`(k3XJ$og8xmK<<{|gxk($o zS0rDNtLI0FZpCI8OPO{F&HN6OR3?6GYX27ocQ@2R<%&xZFF;#t>AJyeQrjm3scKPq zH@?9o$8_Ig+HUnyA%e?ZZYO|Tj3qxhKYE(U+<J6sMF)LHBKymAiAlipd`y{?P;^9v zK+Ddd!Xt;b?fnDS?liy(&IhSkywFcnPY45DO!%orVfMIz*w#4uNvM+QOw_P$v`NFW z*0Du>d*W)8C(|AC*<85-_thoTr}Z>F{PqM5=%9y{$&;$K3iS}Tnig}r_bzxDvpgAj zKjhVF4B1|6aB_@7C)RfR+<77NE%+as^|`G0>j`~Y%fB0?K#(?45dY5IS%L)JgU|jF z=YZtt#wgA~4W+r?!D<D9=Fx3)f5MVw(W!HCTdYswq^FGRuftjL>!w$m;{$+<vU!W+ z;B~Q`bKxGaGppEh;N40nrQ6TfbbiMdCQ~nj)F;%1K<=(g`-{lMZm)3@X+i-P3j~a< zkDv-3ME_omfGxe2hW-I~^8&ihmsPQCjf9wwEhBt;y))!GAm+DJK%JTTZ(-Gepi-N} z8zG=6Vfb`ZzjR?U*fu!5D?>AhxcLh-g>-td9RUQZ$kQ4TcDpsc4>)!MjfJC4H>I~W zy^hsRgVlDErfh(um%P{77k7ron9mV7<^f#$Ih=l12;?P&Uh(d__e>d8pYYevxhD!l z4LBnIA{2fl$i53;!MlpSHXUBBGbQ{@^&4%B@#V93`Y5OzJwy1e4NiTcrA0h}P&Ic6 zRb!-#a2=hnTRUyHBT)*pyxs+m^|I#5f57G{hp6|_OH+<>%VpG`k=-@5t0zTIhY2n_ zzPxH!7lV-{Pu+!{Pi<q@rv|*1JsdG9ORg8otv7<SXHJH^PxV`ht+P~JMzt;ZxP9Ay z1YT-21|<OfJoSdHdEcYHOEOZWYm-c2m#q_g6XzI?cH03K^Bap(6)sc$qM(t%o5NRq z>)|@pYNysN81Aww_|bBcM~IVxdkuz=H-W~_qIuc3(14~W-!!YT%`{?bqyCD*k_RK2 zZ2-NDwhq32IKK%@ZV0~7sgfpadqeAphb|*ZfA_f0c_Eice!18d5(?Idu?^r5iyCZ( zWhVp+W$Ctp6q_;j?$j$F=?-CY16AL8jv_QVHfidB5E65{cwB$DdM;}275^UU9Y5o% zhnxoHaNycE+*_Y~Fub8$$4$oYZ{I1G^kC51DIir4TPU_Y-SZg!kN;>?(_z(&*CwCk zk+-^Nrw&`$eZ*Z0_2KHRiW%x?NC!ACR@8@5f!NLMEvhUAy}Gozac0_v3SM%KUK8zl zA+oJEdBx%@RfqE%(*;u9jd>5+f~)O*rJSoj6Hs@G5O+l%Ebb=IjTLvA<hQ6dQEE7= z`;X1Y+v)`mXu%J|va8jG2qg1DKxZVyESr0Mo@hj_<weEXINF_e@sjBAYL;rKCGS2G z6F~gN2mJ+kjOD}YS$qOi9FuP_@|ux&{0mQjk4&rX@|*062}Z9410G!lw7f7hvO?=x z!=?NSTI5I)Qd!!wUN0gH{Yb>u+EBV|n}D)0C-2sE>X0ZiZm>e!BY&D$R72wftEnng zvy&zdR+e{X;v&Z>bhSUEUeg2eFPitBihis-0c?0|%eZ8Y6qKi$Dn)0f39Ig&07Op; zc3l<nxVZBeP)7Y1;}V=&<oh0UX`#*x{m19GpwGgYU_VjSfhaI4*$DmMsN*NfI|PF` zC4hK*R1B=f@@p2vZf+GG+BqFefOnm2ma8}xF#W-`>$oc8+~L%=9IFcbpWrMDxq*lP z!;Kng3lF(<)rG<|DXTgOj?Rm$BDV*#s5R25!)5(J<>|h>6w1XylhoZg>EH9a`}~$u z_BEL&8YkaxSTP*6NiQ=e)64N`kz-q10Xl3pdh)~RUkSBW4Vz*LY!U5<y61#xzX53v zsjDMDjm=*XYfaS-OL|fHvy`ZZoa)V^z1z~_T|@bIq~VmwoQq3Qb4~oVOCFZbc#DU7 zYWs6}uvj!?WWq)eyF^_aZ)Bbb=<TLw+`a-)=^v{LqaGGrifal<ha8hqQR2O&fAURu zIDT<1Qxv?;?<N4hc>0~J*Fwm4tnP4Bok2xT$_}52TT{edQ*g_pGvZV7i9?}b7Fv|M z21BdY-R^PT0eVHGLMHG3XQGEf*FQUHBU<Y@vB){8OucX9FUe^uJPFmAGS>a9y-<*W zEq_|$Ib0ahEHwNTTyK#2DokVN2yFe0pkTIUh_LOhvDIl~xTNj+#Tz7%Gdr0@qvx_B z-Rr}<P9TeHnhKy$TmbZtL!hhjIod&liuy{9+dKFl`K!q2RfCr}+Le@xUZ~VtzVC3G zD5;cPO?VsXxE56~w_VKlYLHlHk#Bn5?p?gPkwlb9%GxPVgVl_Mb#k0O<D`!@<oxmq zq)tB&BxtC!qq`zoB_1Z$>TXq1`+=UFUhv@CSWxK=v$%t<Q;{bIUmQ|uC?spyF4%@? zEYB-Y1Y<wI$|kdmZv48Aie7I>En2U9Vjh3jAPZvby#Xl<S%vw|(>=~kaz#EXeAs!= zVA6()qD&K}Z~l^?6R;aU2n!zsjk46Nnhu{noGD?mc1)`QQI@Tl&3no}1DH6CQy10P zI;-vf9gA2PD2Ipw&j9dKHx`6kAST|-n$a=Rhj=a~V(kJfNUZrrjcHTJMZ&Z_I-^Fu zquo{o+#eR^ZrUXbiFzMdV7si>ASh!g_Y_~Ivuj+tAg;C$TPG#n#*2_~L%!7DfZO?= zD3x{~`BMyW!aERDs~#{VOKx0OyM>*zJS>l8d`65jjC2uQ`YOPg8EmrqQk{<P>y)@x zDsn}~TVAw7mX5#b686%YvC)MGehIn{mMb(7U+8K7cFU_L{{+G>b8`cF2!-hq!hQn! z{>6MLez5lrpL-Zgy+ktGJD*i0v9AsQ&Y)i}v#LIL`9yEg537v^pX4e4keX5{ZhesG zk|4EGaHIz2vslZ|FW2ar8-FDlc*My3D<2tUs}b>J#N-3E)o`$a%utTuu>KL|deLTd zu8FR{adhbs@L=h$!w=ofcWsU$e{%1g(Xja)&ifX2>fUm52;{|=byg2Rb)mqxHrk4J z6{RqExv~%G9uJ8BfT<7YpVxq(I7q$LFg<$f>(wcFsjIrkcl}AFT@s5=XFdeE@q>14 z<T`F5LMo1b{F0$EjqS;tjbg&pI*mQr7Cj_)buxB7v*XbjYF3|YS(QbqhF2LE3@^$~ z)n1jXcw&YhthgM#Wg}>}`QR#D%Q0Ad;k~?vCW+G;2Hr}mOdC3#B<AdHx%SmF$#_nW zagw3$L*qErn|KQ5c*=iUeG*-Xd#$n97CLzC%U83FDm|n$Iw}b_E6e#%;~Q(`{545> zr!BSe@hR~1)2@x8uqs~c;~0_KSatBLEWJx|>A_-yztXoZDDIcW<)&=q0}S8CthFO; ztvX^Vi$&U4!!!PxGvv!_*CBtgQQ5{;wy(a)H=B?vtNhj8sX+6eUki-4vK@Kd9U4<Q zbH02iX$UB91{?mAj3%rTT1cts+`R4HEFvLrX(Qe5VyJgpj%z&cY$U|t7%CIfvW%}| zR{}yKl^3KN?(aA(vs8yJ$xF<^F$$-!=VUJ}hoME@g!QXG!K=3Q-@F`6v<(ZNy<ll4 z-XHrp5OiWs=Vi@h!gPQUarK$Wiipr`e(CHj1Zi~k!-s;au6U%oQ;&mOQM_eklcCbI zVyb7XLrlDBpo_RO?<~4%v+nX9mEsrN6LH`i+BrPFDP;CJae|Pmm>2d%5$Nyhg_)+! z`evsr$ixgKt@{l;r9ShX+4$`K$k@!5(1sF~Ut9nu1y=T}nt)4FOXFC(t$0T^A!YVz zN?Y`oG<Bu3CVuk;=jADfF)~DI-le5WCx<Q$34%Lb8ldklM!)J>jIR4>WY$MqFGGpH z{QRGX2+!QOi^&$zCZ@{@rk_{6nDQ0J04Q>=3XYVij#vSpOgH+;L(kq^ca$zJJ+(7a z^>iz4EYBbnTB@}ECV4SEKGxuwT(b$5-2jcB88+^^iFmpfzW!tDTq1MyVjeS4^Hw(* zaa(Y6G7cZC+E%$qAX2aA`r|8P=$Gh0g&Vg}$n%LdY^`;c*|T`2*GYdb3|&#lv32N* zM)b|7qPuUTiekc-p5j`K&;rWi$RDU_Bf_-my6gMkmC-i_IWZ%nN7M)74Q+%r^cGt_ z`sq6TQ4rir3iF(?-F*Dtw#Ha0r;TLwv{LQiP>9urWlF!Qa_lQ~9+-5`KUQ_VmA;X8 z8@y$b)#bJOOg#CvU1h(8a(c(mbD33MR*hCQExVPt*q`ZXCj0zTyt~B}8~B*JX?|mw z1>F{tBazlykwp3UtsyOXlZ6)p$L*SlrL^=VeF*C^Wu)K2^lBBU<M*n1;;H)N*SUH| z8!ULcUc!u{Su#3<IQW2+@W$k@?0ntT@(hCvmUkL&fj8;ya^W<}xUX65s|>F^!_2Zp z<6DTWp>cA%v%7z%scPGF27vb`)qVuzDdK!au<b)@RYjvNS4^R#b$1_Hk%&{N)KZ&Q zY)`eN3~^u2C4oMgzH(#anBNVnn1Jf!TF6g;Q(MO3Cf5W>+t|fV9XwjC-=Kh#S6%W_ z`g15teH*1Rlp{YBM_u&htjv}Yw~`z)P;ZEFuO?hD(?hQ;S*DQ~sD0;M?z9(+lhxgR za!HTD#k;3%JxQQA&LW842;;kdgu3;N&^phGdN{X%ahDvV125pW^sGamOZq&lTbiOy z3-)%oOUBYrNbdVWPwJ28SOdAl?xfrxrAo^Qxg4Q4dtBE$tmlhrwW`-PSkKz<`P|X| z;xfI2w06qZxb^X=pA6Z~8Z}AhC-HR4OxNI@@Y$-%y>z3FLiO6q%O>{55%j9{=yKUr zV~;DiK-3wv`Tk$;@=s9wW-fzV4c3$FD=lKm#z$gl-{L*c8|X`*R<BCuR1->v?Ni6^ zBH~N3BQvVyvdzqo74srE#UA>GOgCPa?Jho5ay5~bjqbXeL0GDa=}1HGU3HAOs-OSM zK!#X4t3|7MoQ3+<uDI;;RN_hYqR5yYT4%^oeDkc$QwaWz^BZ~bw&|K$BXslWqP&Yk zC-91;eQ5W}U(T|ft_EJ?FR#m8tS*j8WvsJi7yoRq2(RE8s8<tkUfdhHuDsxg(VvaY z9N5C)SHq`2I9`%5{2I5=<E0E=SVy%+C`Oy43_d(lvyjpKN{09|Z}$;J{GJ^R=#g_f zVz!9*%Ewn9(WMJ@@!2pQb9d?g<~_6cRAp5xrR_1|GyZQ*fM)CI@&wlmWJ$j8d5>!w z=#0Xvox+l3DDxTFt)P5#Cn|C8bz}sZx7AyPH5e$ZEVeAVSc>1Py?U*hripzGUC=bz z)w_jIMuk7lH&KCHZ8jOs8UOf|HM~&UYo=6UXoC;GhKqqjUBjE0I?Ba=2{*ij2k=&? za94PXokb@bCUb@v{%f6;{}~vb;g#c@5011_s@&w=HC%ls&K5`aCnQk<^`0fQi}yTD z^6}w1bY^#t8HHFGDcSNc;$OtWaox_Evm0J6bA-WP$SWp#+3`=*+`bx$T5dk(#WbZN zOOU&B#b4jJ(cgwey`xXcf7nifC930QHom8LE(8~?XZvK~3rwn&FL-pkFWi-t*Wb5` zju>++Sjs0<DQ2ZCJrObe>{{ctvCqHKkfmn&k#bc<RH`x^vHNI$e2czyx3kuHz|O!H zD#ucEX4(}miZtr&qMcQsy?vLv>6q`4?lQ!NXU@;S8#m(+mdVOCn252YS(Eavfh|q{ z@}BZ%<WFH1`7Th`Uw)(q{yr(w*@M|c+cb9pp4IS**2+H{XX%}CvXjMz<Xl}?!<&<3 zJc~5k==3f6`EI(B>n!K$+Q<uX?wgJ}_fqn-MJ%1tfxf|fI#sOXmDBLjc{j?rBX2X1 z5~%pkKH%5l-`)Qe$nt7KHTeEl8+zxz+R*>XW)Hj)c=C`}wR+Up)7|!^v#q<Xi=C_J zeOE7QC)7D72Ny@%SG>l<lUEQLa_j(~#(xE&=H^ZkCTE~(+oOz2x-k&t%-b=W(M`H% z{mtc^9jb`^2CwCu2@Xr?1`5?`Ht-^~;>5)8<va7ojy-vM=iISJcXZx935C`P{d4j6 z*Q+i7a_X^9z^{*<zQ2Px(0}0xaUWJpOY9t;K&I%NYe`B>8CRgI?ay$^J+Q#lz4AYw zztsKsD?pkh{5dhv^f+Nn!ZK%}>v4O508!&w)%uQPC$Cz}LQM!92Hc7U8*;t<{*|}? z>m*zIWi>=nJ2?|r{PVi+hXZamo>_eYG`VTr9=djN)=_oh_{T?XHRnc=8Yk`@IbZdq z?;ENzA0#<ur%+}_w;#B+_|5+8ua)luX4qv2P(#%aR7%LKAsp)?9+Ri-l{|7wv2Ln) z+?0FM+HYm{=X};&jnOkp_{KuMs-v3i`6fSX-ksy%S<LyasenY&%1ih8WQR_C)L2b$ zn5Q!EQhT>|FIq`GV0Y7?<c@l!+vk)l@9>}Hq(5@gXbz@Gddt3_Q~KTgxY}B6{3O~) zfAcu(>(M!j%|BHNWtvCd9Jgsm#?Og1W<<{U4&3ZYuyxh4?&u{A-mJ^*FhsZqE?yvm zlx0sc{zU4<%$*_w12z5K5QN?i%&d68-=4P>`>{Mbjbw@#9lF&@Mz~JKWWUtQyL0gR zKZXBZxFa)jK&o_^nW*htt`Mp8vFX+K=?};K4be9J{?G1|pnecS%O)e&yY2QofAXn{ zZOyVejoA)8A1Q_k<ts<{*0dUvoL=aEPgO+~eQ#5V0BO6)p-Li;!qvZ=x~S3!Tt+l< z)iqC^Ug)~=1Tr?*xA{o1+oQvJyl?!{`FEf+o%ny!MWuBYZ@Ni}I@|?WrvH;GA?F~s z@};5lOpNB@<)>qhLM-5UBfn-|_B+1wc)T$9Wn|0xfqn4sKP~-#tFpwyZzjfgPUo(N zr%8k7pPxDOM7?+{UMU0SK$K_=4@(VB)~fe5Vg-J`J`z)_)D096{+sywVB3?~S$k7w zCMn$<61~_hDWQ<F@MYd_yCm)Dz|Yx+T$9_%aqH1G1N?-{)2{CrE1#%dr-nl#@pN%r z^6%lK_#;>N+hS#*#kz$F2IBX=9EI>FH6g4B>qmUw-j#{`DPtPfV<%_Z^Zluu@D?EN zNWAKLE8o2im#wp6m-2p!o3$3FSncJu2_A(E^am@x3^5YbOC=vUEMsS&8vf24-__+R zH=O$HF?ivrXj`s}VYe^8@r6s+8U2da1ztB|g||-C&6HiHjzgL%S@)Y-V?RBLPldnL zX|m7&89u2w&sg}A0SkJfQ<eGt*Q1(LkIw9JVPT!ERN>1V*Wd|P+v#qxFC{;S{hQd% z|KUu^Y{q={$K~HiZ+gmu6pWd@eS8B7HJ^z9$r#M@^-ro2(x)5tv5wWJk3D-eVR*}x z`2aO$aC{xvS(&aTKJEM=z`Rw65rD|vJ9!tbs1>Y20!1icQ(8i!t+FcEd48q8)oyow z{x!2Ad3Ma_&J3czrm_&Tu-6x2Xw@S*+-!O)Z1as;*MmqW=~XwI`o~|jXMW_|tAF}Y zyD%Z6F`<9<p-*GFo^ok!^(9hn{q(=*rJ7i>*f6P)h8tFGLVS`4wkA3LZK8ETmbkFT zx0hYSO9}R9nU<<YcVarN(KoKPxY*qHFg*J)$KQwNfrNxa0qZ3Q-tlwxcgC4chQ-Og zlArLgi<*}pNPVa(_=+O+#e5A>;;n2f;aKkl^<$EEfWmJ#zNuY=g5`1qU}ANM7rtk1 zHr@KF*wT2)+I)Q6?LxoJbn9lg&iwU`1Nt_+8so+1Kgs879w6S~b%dV*>ave2b-rL~ zvCjH7kh@MN`Bz4q5bAY;H22r#RZUjj2R_WQxf)5y7E0Cku8-l*Ixa7WRZJ|tG`sX7 zw%G9tEg^BhuM}31V^X4LNk43yGn{H~*5}abB+y@S@izhg&s-5<`KIBpp%<1`2HpMr zv?C|Uesu5%t(-3Ev|I4t!HtAN0w!61TqBG9`^(g;SC(S$mf7dn_T@b-L7h0VBY66k ztlx**tQ&?A58UAD@7`jr$sfG8Vt}h37d?`%aE<PBKS=H9`^e(T@t01Kf1)oNQYG|0 zNYMUi`;>m;bR$YCbaScjKvf^aR65`Fxy~D>y?3vrS}KorivII2;`4cZ-0tcBtCchV zXTtI0xZB8@t64L`99t;L$W@;l6CsoPNbWJB+(*SSVJXR#+srU$xh1ABw_%~A!<9R= zNF*PRe)|6Y`Thg1$Lp8Z>xb9-{dg|W5d5=WyyNhzN_B2J{Yf|6+is*|0UI;CAl`to z6a#L1@V74P@5&97g@6@Gwj{(h3=tg#P#F4U_~_v))!?%!-d_g3?URO^=#X;p_zlx6 z+pQ0u>Q-ji^TIy)DcA3~JpekHwLCm>YnE<Ok;i`hprsd)Z`tQOSYaqB7O#>JWPIxd zv>}ci8k{CWj?_=lU(6nnj#xZ-fBtrntY9=Zd;QL21!!n)OTa;l0mVct+67seKLaO* zXeGb`^jyo=i0xn4e;pR}Z9EE@q3R7ZVN>TSS++Zjs*DRDkQCyKnt#Gp@PMNP^paih zi||gXr-1NldNeW~^c>Pin9&z_Z0FeajpckoZ<7B*)+=$@B)XUgH$L2EQwGJGiHOw| z3{{F-A_Nry@^K{bdn)<3F{)4cmZlm8FVZv*uQbiPod`}x57w$NnpNJ37nW2eEfNqs z=e!Ac1$*-c#xE5LlPLaIIMKrnR0kzpTnm{a!lPrhPUm@ndxRhVTb_vCP3N||lb2;( zD_?s^sKq0mi{%)O5_!}V|KM&}o0v|FdlHsm%aFw*tIomgT@x3z9*Cs1S2-r&mxXeD zbzL3R{|rH&SdYJZU~2JZU|eetI3CN%E}CN5O)wTgW&lO9)SF7KRY+_4(~1DiGHPTY zn4-aW0)=rWb?1$Sn}Yz4{pgD1Vx<<}*q8w?T8fHQW5l7%omS66+u%m9ID>L7o1-aU z6v%!K0%<JoY)GF~$~hCf$?tZ9SOF&~>PxgyPzKp`9@}8uiQH?KG%u`NcJ7cr%;;Bl z@jWb7iR<kY-u?~ULr~Mh@qWs8bx4AR0HKt6>NQS`BUsm2->2UKv&4|7GkD?#H@C~{ ze)Ldix;l_!jW^%d91$T6lLYF-Md!)ef*yJ0l5hbdS&PeIs^?eSc4lL!BQrHeKZIcs zVN30M%%b=EKi@|kND&|&UmSpYg_b#cF_PBh)5Y#fafL)c(IiF+U2n4^Bp8lc#!_Ac z8U%Tti_xF+>p1QI1Kl_}S@S7XLPUfU6b;*drP6S4eBG5YeK6%miDJYsYKl4snq7|Z zS4*-2l(a%+k)EuiWy`?UyPeU5k2kGd18J>BGc_cd>h*=qPAI9jK4jYF+XA9ruTSJe zALWCDgkjUgv}cnaWM)v!OGW*C#C|9yuO~>~%{0WyO%=|7%k{*p>eDRzOk(04zYFyM zT2DXDuq*bQff!Xce4r6n?~k#`*{|8Hx7()55;0A^XRJ)cGJ(gm=ssJ7QfFaK_aYGZ zk1-u7V{gRQmGqd(&KXiP?~}<@j`n*TBOkszraF_o?X`O{8m<Czp;~iAk&pXax@<D- z-g;>p_^igm<0A*$_=u;quH5W%(rQLR*I+WxtNNsa5?VeUe|}=$tZ0M2=`eXfE|Y7d zKp5!V4kYs+J%$g_6Lwqu(X6ANWK%qtYnZI4Ms_RyKKoT#Ch{{poP8UkGwpZc8;?c| z>GaFAf`zq6xI!Mw<KZdS62kY;mFq3Vx>IT?CC`3OW!}4sME;QnNiE)uT;~pcYB+bg zP!;(5>fx*LDM+|kebLiS)?%Ktew^SczsRm9pal=KQrXbnw3@6+t#*~p-%jBASn(QO zjD5K_Dn)47zU-3vT||CiLJ;iak3ZHCacBMJicSYqZA%mkX!ZYGOg-j2axu@JKz!p5 zab6rqr3`vMrL1GJNwu?2I?msb{1LNyJbOyek_TskHFp`~=KaXhXb6r;>@B&6{ndsf zm5z0S9bg`Nd|z4%xgji-VH-G!9J#71O=n;xq}<V}F15`iwm(ju)D;!#aZR+otWeuV zGTaRs%yiJ`5%ztDgfQ6+d!`>srnd<Ouogt$y;qXx<elG!I^@3Ix8A6|Aj|fAO<5-% z6*e+{AGy7B;rIO4%~X8j1)8=jD?sCCzrj345GYT%cv6^&a676@tgczw4{mXkupw^e zYT+W!QnicI*X`~1!2=S>%a<Tu?xC6#;_#*PDbcRbbq~$fl`7+;a6{8g(lsyR*f3E^ zY5nB$k7+l#qAw!CO{O$S5Tnq|8&V#`OQpVcau>@)CLFm(4ocvkL#@8@F!r9r%7`7R ziI@)I2|{E?a<6mOg=;%Kw#uugn24xe?{e>-walCqDQ3kKibbp`w%t^(C{*p5`Gc-& z7H3zFA2>>;&E5fl!i>s{mk2)_LR<qJcg@sHByn^IZJ>Dsz&$DR^MqoB+f3+E*VNL( zKIxW6iC&O=)q&brW_^kfCN7Pt#o$;JbvK`}8tj9OKYdXid}fik;yPv{m+BonGw+PK zebG=H{{*$0r@~vy=Lb*vt$hD}O5jlL6P&v8yDL)iRdRH2U0n}m;0Mq4`EzC2fLe*6 zTb|J5%~&s#Pr7H+O4~-c=k)B5T-EEVhK~Uq2vMhj%HC&kwrfFdM*tl+wD+U^&rzXX z3z8bODSEq;;9Q<S3yKHRzUb5Pnc|*CW+e*kPan&>>hQMmMbiN9o6icNHud|!U-3?= zdd^f?%VA0>=e^MQ`<v%L(q8%Et5;|I#-LksuAh{A^0j@A=)UJ(Y8_P%8j>vVU(C+) zJDZbXHFY`R>#w_|S`zrl;mJ`bx{lT0RcuT+@qNCVhRydZM(R#fFUb~`O^;0i2Y!Ei zQ;G~P%*@n=-of-WwVE{*jzp|3KYz6w$aBefE7=a5h8ivmeX*XUTT!ehE7s;?DGh*& z9idttUm}w=Gv<Wy)tmMK>Kj%3*T?3E47n%Vb!M|XpOyq<HucQs^?+ipa^(3~ZrQt- zUvoU{YSVnpXRuNtjuYnlotss^m{rC<c$&#;!}C(9^VzKdb$cbLq?Y>PVJ!;<<dhnh zLLUmOjqX2{BQ84NHU@;BwG>E9=~Kij%ANIT?P6+J3NGVPkV2k;@Udmg+sw#%gOdED z!W`lLf`5Exy3KgcjUH?zE?ggdKBqkPdX#6l7B;1oCizc;RT%m34!$v^#URk43^UuA zOM_ok(B7Y5QAiAbo9x><@z3p8bji_QtX!#^H(+KdH;CgjhYv$+l~^S5AHQw@abui# z>gpP0OC*&}t+LP&qwVRjx7fb5kW3VAR(~I>*D~IT7AmB^@vZk5pSmjp$B8y4v3z&X zP#ao{;>RDFHWt?|_1E-P5x*X+fle3Mawpdmsp@er_y(0RVtGVQUN^saCqyLF1;Wlu z{<!zRxEtdWY}u{%s`22Uz%2VfEVAOjyO<#_Wra_3udJJ_d4Y*Ww4dE=<trI1oRPp@ z@r}aOy-RHZ$v9;M246Aqqq8}|c43W_Z!XDacFB*zHO<;m$VRGj57_t+^Lj5JUEQT= zkgQ{-#=*+#)n)#Ci1MuQ2`yv3+X5V++}Pk#kp}A83ew#+>$1XoLynMggK`;}f_m+! zjII3CfO!o$!5XJtUwMRm2BJj=@rRQdZ!dkgTu}Z}yRQ7*GyYP$4*6Jr_l*BQ?K%`V f97g_c7_CD&*S{|L|Hmcmzq>s2XeJ$oN@M#MC*69U literal 0 HcmV?d00001 diff --git a/Reconstruction/Jet/JetMomentTools/src/JetCaloCellQualityTool.cxx b/Reconstruction/Jet/JetMomentTools/src/JetCaloCellQualityTool.cxx new file mode 100644 index 00000000000..0041895fede --- /dev/null +++ b/Reconstruction/Jet/JetMomentTools/src/JetCaloCellQualityTool.cxx @@ -0,0 +1,87 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#include "JetCaloCellQualityTool.h" +#include "xAODJet/JetAccessorMap.h" +#include <iostream> + +JetCaloCellQualityTool::JetCaloCellQualityTool(const std::string& name) + : JetCaloQualityTool(name) + , m_attSuffix("FromCells") + +{ + declareProperty("CellSuffix", m_attSuffix); + + declareProperty("LArQualityCut", m_LArQualityCut=4000); + declareProperty("TileQualityCut", m_TileQualityCut=254); + +} + +StatusCode JetCaloCellQualityTool::initialize(){ + ATH_MSG_DEBUG(" initialize() "); + + // In this tool we're using the cell-based calculators + m_calcProcessor = & m_jetCalculations; + + // Define and use a macro to compactify the addition of calo calculators + // The macro defines a JetCalculator as 'c' in case addition setting are needed. +#define ADDCALCULATOR( klass ) klass *c = new klass(); c->setName( c->name() + m_attSuffix); m_cellCalculators.addCellCalculator(c) + + // convert names passed as property into calo calculators + for( std::string & calcN : m_calculationNames){ + + if ( calcN == "LArQuality") { + ADDCALCULATOR( jet::JetCalcQuality_fromCells ); + c->LArQualityCut = m_LArQualityCut; // c is defined in the macro. + } else if ( calcN == "TileQuality") { + ADDCALCULATOR( jet::JetCalcQuality_fromCells ); + c->setName( "TileQuality" ); + c->TileQualityCut = m_TileQualityCut; // c is defined in the macro. + c->includeTile = true; + c->includeLAr = false; + } else if ( calcN == "Timing") { + ADDCALCULATOR( jet::JetCalcTimeCells_fromCells ); + } else if ( calcN == "QualityHEC") { + ADDCALCULATOR( jet::JetCalcQualityHEC_fromCells ); + } else if ( calcN == "NegativeE") { + ADDCALCULATOR( jet::JetCalcNegativeEnergy_fromCells ); + } else if ( calcN == "AverageLArQF") { + ADDCALCULATOR( jet::JetCalcAverageLArQualityF_fromCells ); + } else if ( calcN == "Centroid") { + ADDCALCULATOR( jet::JetCalcCentroid_fromCells ); + } else if ( calcN == "N90Cells") { + ADDCALCULATOR( jet::JetCalcnLeadingCells_fromCells ); + c->setName( xAOD::JetAttributeAccessor::name( xAOD::JetAttribute::N90Cells ) ); + } else if ( calcN == "BchCorrCell") { + ATH_MSG_ERROR(" No BchCorrCell implemented yet using CaloCell direct access"); + return StatusCode::FAILURE; + } else if (calcN == "FracSamplingMax") { + m_doFracSamplingMax = true; // no caculator, as this is a special case. + } + + }// end loop over m_calculationNames + + + // Define OOT calculators. + for( double & timeCut : m_timingTimeCuts){ + + // build the moment name from the base-name and the value of the timing cut + std::stringstream s; + s << std::setprecision(0) << std::fixed << "OotFracCells" << timeCut; + + jet::JetCalcOutOfTimeEnergyFraction_fromCells* c = new jet::JetCalcOutOfTimeEnergyFraction_fromCells(); + c->setName(s.str()); + c->timecut = timeCut; + m_cellCalculators.addCalculator( c ); + } + + + for(size_t i=0;i<m_cellCalculators.numCalculators();i++) { + const jet::JetCaloCalculator* calc = m_cellCalculators.at(i); + ATH_MSG_INFO( "Will calculate cell calo attribute : "<< calc->name() ); + } + + return StatusCode::SUCCESS; + +} diff --git a/Reconstruction/Jet/JetMomentTools/src/JetCaloCellQualityTool.h b/Reconstruction/Jet/JetMomentTools/src/JetCaloCellQualityTool.h new file mode 100644 index 00000000000..33e3174ec11 --- /dev/null +++ b/Reconstruction/Jet/JetMomentTools/src/JetCaloCellQualityTool.h @@ -0,0 +1,53 @@ +// this file is -*- C++ -*- + +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +/** + \class JetCaloCellQualityTool + \brief Calculates calorimeter based variables for jet quality + using CaloCell contained in jet object directly. + + This class is identical (same properties) as JetCaloQualityTool except it uses CaloCell based calculations instead of the cluster moments based calculations. + + *IMPORANT* : cell based calculations often differ from cluster based calculations (because of how cluster + moments are calculated/stored). + + This tool therefore give different jet attribute names to these calculations. + Ex : "NegativeE" -> "NegativeEFromCells" + The suffix can be changed by the CellSuffix property (default to "FromCells"). + + + This class is available only in Athena. + + \author P-A Delsart + \date (run 2 implementation) November , 2014 +*/ +#ifndef JETREC_JETCALOCELLQUALITYTOOL_H +#define JETREC_JETCALOCELLQUALITYTOOL_H + +#include "JetMomentTools/JetCaloQualityTool.h" + +#include "JetUtils/JetCaloCellQualityUtils.h" + +class JetCaloCellQualityTool : public JetCaloQualityTool { + ASG_TOOL_CLASS0(JetCaloCellQualityTool); +public: + JetCaloCellQualityTool(const std::string & name); + virtual StatusCode initialize(); + + protected: + + /// This objects holds a list of cell-based calculators + jet::JetCaloCellCalculations m_cellCalculators; + + /// Suffix added to the name of the jet attribute + std::string m_attSuffix; + + // parameters for Quality cuts + int m_LArQualityCut; + int m_TileQualityCut; + +}; +#endif diff --git a/Reconstruction/Jet/JetMomentTools/src/JetCaloEnergies.cxx b/Reconstruction/Jet/JetMomentTools/src/JetCaloEnergies.cxx deleted file mode 100644 index 06b160c4747..00000000000 --- a/Reconstruction/Jet/JetMomentTools/src/JetCaloEnergies.cxx +++ /dev/null @@ -1,57 +0,0 @@ -/* - 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 deleted file mode 100644 index ed7b98a679d..00000000000 --- a/Reconstruction/Jet/JetMomentTools/src/JetCaloQualityTool.cxx +++ /dev/null @@ -1,177 +0,0 @@ -/* - 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/JetLArHVTool.cxx b/Reconstruction/Jet/JetMomentTools/src/JetLArHVTool.cxx deleted file mode 100644 index 2697e970b07..00000000000 --- a/Reconstruction/Jet/JetMomentTools/src/JetLArHVTool.cxx +++ /dev/null @@ -1,105 +0,0 @@ -/* - 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 index a5fefd200f9..7bf1c891334 100644 --- a/Reconstruction/Jet/JetMomentTools/src/components/JetMomentTools_entries.cxx +++ b/Reconstruction/Jet/JetMomentTools/src/components/JetMomentTools_entries.cxx @@ -3,36 +3,45 @@ #include "JetMomentTools/JetCaloEnergies.h" #include "JetMomentTools/JetCaloQualityTool.h" +#include "../JetCaloCellQualityTool.h" #include "JetMomentTools/JetWidthTool.h" #include "JetMomentTools/JetBadChanCorrTool.h" #include "JetMomentTools/JetVertexFractionTool.h" +#include "JetMomentTools/JetVertexTaggerTool.h" #include "JetMomentTools/JetTrackMomentsTool.h" #include "JetMomentTools/JetMuonSegmentMomentsTool.h" #include "JetMomentTools/JetPtAssociationTool.h" #include "JetMomentTools/JetIsolationTool.h" #include "JetMomentTools/JetLArHVTool.h" +#include "JetMomentTools/JetOriginCorrectionTool.h" DECLARE_TOOL_FACTORY(JetCaloEnergies) DECLARE_TOOL_FACTORY(JetCaloQualityTool) +DECLARE_TOOL_FACTORY(JetCaloCellQualityTool) DECLARE_TOOL_FACTORY(JetWidthTool) DECLARE_TOOL_FACTORY(JetBadChanCorrTool) DECLARE_TOOL_FACTORY(JetVertexFractionTool) +DECLARE_TOOL_FACTORY(JetVertexTaggerTool) DECLARE_TOOL_FACTORY(JetTrackMomentsTool) DECLARE_TOOL_FACTORY(JetMuonSegmentMomentsTool) DECLARE_TOOL_FACTORY(JetPtAssociationTool) DECLARE_TOOL_FACTORY(JetIsolationTool) DECLARE_TOOL_FACTORY(JetLArHVTool) +DECLARE_TOOL_FACTORY(JetOriginCorrectionTool) DECLARE_FACTORY_ENTRIES(JetRec) { DECLARE_TOOL(JetCaloEnergies) DECLARE_TOOL(JetCaloQualityTool) + DECLARE_TOOL(JetCaloCellQualityTool) DECLARE_TOOL(JetWidthTool) DECLARE_TOOL(JetBadChanCorrTool) DECLARE_TOOL(JetMuonSegmentMomentsTool) DECLARE_TOOL(JetTrackMomentsTool) DECLARE_TOOL(JetVertexFractionTool) + DECLARE_TOOL(JetVertexTaggerTool) DECLARE_TOOL(JetPtAssociationTool) DECLARE_TOOL(JetIsolationTool) DECLARE_TOOL(JetLArHVTool) + DECLARE_TOOL(JetOriginCorrectionTool) } -- GitLab