JetPlusTracks/plugins/JetPlusTrackCorrector.h

0001 #ifndef RecoJets_JetPlusTrack_JetPlusTrackCorrector_h
0002 #define RecoJets_JetPlusTrack_JetPlusTrackCorrector_h
0003
0004 #include "FWCore/Framework/interface/ConsumesCollector.h"
0005 #include "FWCore/Utilities/interface/InputTag.h"
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/Event.h"
0008
0009 #include "DataFormats/Common/interface/Handle.h"
0010 #include "DataFormats/Common/interface/ValueMap.h"
0011 #include "DataFormats/Common/interface/View.h"
0012 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0013 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
0014 #include "DataFormats/JetReco/interface/JetTracksAssociation.h"
0015 #include "DataFormats/Math/interface/LorentzVector.h"
0016 #include "DataFormats/MuonReco/interface/Muon.h"
0017 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0018 #include "DataFormats/TrackReco/interface/Track.h"
0019 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0020 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0021 #include "DataFormats/VertexReco/interface/Vertex.h"
0022 #include "DataFormats/PatCandidates/interface/Electron.h"
0023 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0024 #include "DataFormats/PatCandidates/interface/Muon.h"
0025
0026 #include "boost/range/iterator_range.hpp"
0027 #include <sstream>
0028 #include <string>
0029
0030 namespace edm {
0031   class Event;
0032   class EventSetup;
0033   class ParameterSet;
0034 }  // namespace edm
0035
0036 // --------------------------------------------------------
0037 // -------------------- Helper classes --------------------
0038 // --------------------------------------------------------
0039
0040 namespace jpt {
0041
0042   /// Container class for response & efficiency maps
0043   class Map {
0044   public:
0045     Map(std::string, bool verbose = false);
0046     Map();
0047     ~Map();
0048
0049     uint32_t nEtaBins() const;
0050     uint32_t nPtBins() const;
0051
0052     double eta(uint32_t) const;
0053     double pt(uint32_t) const;
0054
0055     uint32_t etaBin(double eta) const;
0056     uint32_t ptBin(double pt) const;
0057
0058     double value(uint32_t eta_bin, uint32_t pt_bin) const;
0059
0060     double binCenterEta(uint32_t) const;
0061     double binCenterPt(uint32_t) const;
0062
0063     void clear();
0064     void print(std::stringstream& ss) const;
0065
0066   private:
0067     class Element {
0068     public:
0069       Element() : ieta_(0), ipt_(0), eta_(0.), pt_(0.), val_(0.) { ; }
0070       uint32_t ieta_;
0071       uint32_t ipt_;
0072       double eta_;
0073       double pt_;
0074       double val_;
0075     };
0076
0077     typedef std::vector<double> VDouble;
0078     typedef std::vector<VDouble> VVDouble;
0079
0080     std::vector<double> eta_;
0081     std::vector<double> pt_;
0082     VVDouble data_;
0083   };
0084
0085   inline uint32_t Map::nEtaBins() const { return eta_.size(); }
0086   inline uint32_t Map::nPtBins() const { return pt_.size(); }
0087
0088   /// Generic container class
0089   class Efficiency {
0090   public:
0091     Efficiency(const jpt::Map& response, const jpt::Map& efficiency, const jpt::Map& leakage);
0092
0093     Efficiency() = delete;
0094
0095     typedef std::pair<uint16_t, double> Pair;
0096
0097     uint16_t nTrks(uint32_t eta_bin, uint32_t pt_bin) const;
0098
0099     double inConeCorr(uint32_t eta_bin, uint32_t pt_bin) const;
0100     double outOfConeCorr(uint32_t eta_bin, uint32_t pt_bin) const;
0101
0102     uint32_t nEtaBins() const;
0103     uint32_t nPtBins() const;
0104
0105     uint32_t size() const;
0106     bool empty() const;
0107
0108     void addE(uint32_t eta_bin, uint32_t pt_bin, double energy);
0109     void reset();
0110
0111     void print() const;
0112
0113   private:
0114     double sumE(uint32_t eta_bin, uint32_t pt_bin) const;
0115     double meanE(uint32_t eta_bin, uint32_t pt_bin) const;
0116
0117     bool check(uint32_t eta_bin, uint32_t pt_bin, std::string name = "check") const;
0118
0119     typedef std::vector<Pair> VPair;
0120     typedef std::vector<VPair> VVPair;
0121     VVPair data_;
0122
0123     const jpt::Map& response_;
0124     const jpt::Map& efficiency_;
0125     const jpt::Map& leakage_;
0126   };
0127
0128   inline uint32_t Efficiency::nEtaBins() const { return response_.nEtaBins(); }
0129   inline uint32_t Efficiency::nPtBins() const { return response_.nPtBins(); }
0130   inline uint32_t Efficiency::size() const { return data_.size(); }
0131   inline bool Efficiency::empty() const { return data_.empty(); }
0132
0133   /// Tracks associated to jets that are in-cone at Vertex and CaloFace
0134   class JetTracks {
0135   public:
0136     JetTracks();
0137     ~JetTracks();
0138     void clear();
0139     reco::TrackRefVector vertex_;
0140     reco::TrackRefVector caloFace_;
0141   };
0142
0143   /// Particles matched to tracks that are in/in, in/out, out/in at Vertex and CaloFace
0144   class MatchedTracks {
0145   public:
0146     MatchedTracks();
0147     ~MatchedTracks();
0148     void clear();
0149     reco::TrackRefVector inVertexInCalo_;
0150     reco::TrackRefVector outOfVertexInCalo_;
0151     reco::TrackRefVector inVertexOutOfCalo_;
0152   };
0153
0154   //class Sum {
0155   //public:
0156   //  Sum() { theResponseOfChargedWithEff = 0.; theResponseOfChargedFull = 0.;}
0157   //  void set( double a, double b) const {theResponseOfChargedWithEff=a;theResponseOfChargedFull=b;}
0158   //  double theResponseOfChargedWithEff;
0159   //  double theResponseOfChargedFull;
0160   //};
0161
0162 }  // namespace jpt
0163
0164 // -------------------------------------------------------
0165 // -------------------- JPT algorithm --------------------
0166 // -------------------------------------------------------
0167
0168 /**
0169    \brief Jet energy correction algorithm using tracks
0170 */
0171 class JetPlusTrackCorrector {
0172   // ---------- Public interface ----------
0173
0174 public:
0175   /// Constructor
0176   JetPlusTrackCorrector(const edm::ParameterSet& iPS, edm::ConsumesCollector&& iC);
0177
0178   /// Destructor
0179   virtual ~JetPlusTrackCorrector();
0180
0181   // Typedefs for 4-momentum
0182   typedef math::XYZTLorentzVector P4;
0183   typedef math::PtEtaPhiELorentzVectorD PtEtaPhiE;
0184   typedef math::PtEtaPhiMLorentzVectorD PtEtaPhiM;
0185
0186   /// Vectorial correction method (corrected 4-momentum passed by reference)
0187   double correction(const reco::Jet&,
0188                     const reco::Jet&,
0189                     const edm::Event&,
0190                     const edm::EventSetup&,
0191                     P4&,
0192                     jpt::MatchedTracks& pions,
0193                     jpt::MatchedTracks& muons,
0194                     jpt::MatchedTracks& elecs,
0195                     bool& validMatches);
0196   double correction(const reco::Jet&,
0197                     const reco::Jet&,
0198                     const edm::Event&,
0199                     const edm::EventSetup&,
0200                     const reco::TrackRefVector& tracksinvert,
0201                     const reco::TrackRefVector& tracksincalo,
0202                     P4&,
0203                     jpt::MatchedTracks& pions,
0204                     jpt::MatchedTracks& muons,
0205                     jpt::MatchedTracks& elecs);
0206
0207   /// Scalar correction method
0208   double correction(const reco::Jet&,
0209                     const reco::Jet&,
0210                     const edm::Event&,
0211                     const edm::EventSetup&,
0212                     jpt::MatchedTracks& pions,
0213                     jpt::MatchedTracks& muons,
0214                     jpt::MatchedTracks& elecs,
0215                     bool& validMatches);
0216
0217   /// Correction method (not used)
0218   double correction(const reco::Jet&) const;
0219
0220   /// Correction method (not used)
0221   double correction(const P4&) const;
0222
0223   /// For AA - correct in tracker
0224
0225   //  double correctAA(  const reco::Jet&, const reco::TrackRefVector&, double&, const reco::TrackRefVector&,const reco::TrackRefVector&, double&) const;
0226   double correctAA(const reco::Jet&,
0227                    const reco::TrackRefVector&,
0228                    double&,
0229                    const reco::TrackRefVector&,
0230                    const reco::TrackRefVector&,
0231                    double,
0232                    const reco::TrackRefVector&) const;
0233
0234   /// Returns true
0235   bool eventRequired() const;
0236
0237   /// Returns value of configurable
0238   bool vectorialCorrection() const;
0239
0240   // ---------- Extended interface ----------
0241
0242   /// Get responses/sumPT/SumEnergy with and without Efficiency correction
0243
0244   double getResponseOfChargedWithEff() { return theResponseOfChargedWithEff; }
0245   double getResponseOfChargedWithoutEff() { return theResponseOfChargedWithoutEff; }
0246   double getSumPtWithEff() { return theSumPtWithEff; }
0247   double getSumPtWithoutEff() { return theSumPtWithoutEff; }
0248   double getSumEnergyWithEff() { return theSumEnergyWithEff; }
0249   double getSumEnergyWithoutEff() { return theSumEnergyWithoutEff; }
0250   double getSumPtForBeta() { return theSumPtForBeta; }
0251
0252   /// Can jet be JPT-corrected?
0253   bool canCorrect(const reco::Jet&) const;
0254
0255   /// Matches tracks to different particle types
0256   bool matchTracks(const reco::Jet&,
0257                    const edm::Event&,
0258                    const edm::EventSetup&,
0259                    jpt::MatchedTracks& pions,
0260                    jpt::MatchedTracks& muons,
0261                    jpt::MatchedTracks& elecs);
0262
0263   void matchTracks(const reco::Jet&,
0264                    const edm::Event&,
0265                    const edm::EventSetup&,
0266                    const reco::TrackRefVector& tracksinvert,
0267                    const reco::TrackRefVector& tracksincalo,
0268                    jpt::MatchedTracks& pions,
0269                    jpt::MatchedTracks& muons,
0270                    jpt::MatchedTracks& elecs);
0271
0272   /// Calculates corrections to be applied using pions
0273   P4 pionCorrection(const P4& jet, const jpt::MatchedTracks& pions);
0274
0275   /// Calculates correction to be applied using muons
0276   P4 muonCorrection(const P4& jet, const jpt::MatchedTracks& muons);
0277
0278   /// Calculates correction to be applied using electrons
0279   P4 elecCorrection(const P4& jet, const jpt::MatchedTracks& elecs) const;
0280
0281   /// Get reponses
0282
0283   //  double theResponseOfChargedWithEff;
0284   //  double theResponseOfChargedFull;
0285   //  double setResp(double a) const {return a;}
0286   //  double getChargedResponsesWithEff() { return theResponseOfChargedWithEff;}
0287   //  double getChargedResponsesFull() { return theResponseOfChargedFull;}
0288
0289   // ---------- Protected interface ----------
0290
0291 protected:
0292   // Some useful typedefs
0293   typedef reco::MuonCollection RecoMuons;
0294   typedef reco::GsfElectronCollection RecoElectrons;
0295   typedef edm::ValueMap<float> RecoElectronIds;
0296   typedef reco::JetTracksAssociation::Container JetTracksAssociations;
0297   typedef reco::TrackRefVector TrackRefs;
0298
0299   /// Associates tracks to jets (overriden in derived class)
0300   virtual bool jetTrackAssociation(const reco::Jet&, const edm::Event&, const edm::EventSetup&, jpt::JetTracks&) const;
0301
0302   /// JTA using collections from event
0303   bool jtaUsingEventData(const reco::Jet&, const edm::Event&, jpt::JetTracks&) const;
0304
0305   /// Matches tracks to different particle types (overriden in derived class)
0306   virtual void matchTracks(const jpt::JetTracks&,
0307                            const edm::Event&,
0308                            jpt::MatchedTracks& pions,
0309                            jpt::MatchedTracks& muons,
0310                            jpt::MatchedTracks& elecs);
0311
0312   /// Calculates individual pion corrections
0313   P4 pionCorrection(
0314       const P4& jet, const TrackRefs& pions, jpt::Efficiency&, bool in_cone_at_vertex, bool in_cone_at_calo_face);
0315
0316   /// Calculates individual muons corrections
0317   P4 muonCorrection(const P4& jet, const TrackRefs& muons, bool in_cone_at_vertex, bool in_cone_at_calo_face);
0318
0319   /// Calculates individual electron corrections
0320   P4 elecCorrection(const P4& jet, const TrackRefs& elecs, bool in_cone_at_vertex, bool in_cone_at_calo_face);
0321
0322   /// Calculates vectorial correction using total track 3-momentum
0323   P4 jetDirFromTracks(const P4& jet,
0324                       const jpt::MatchedTracks& pions,
0325                       const jpt::MatchedTracks& muons,
0326                       const jpt::MatchedTracks& elecs) const;
0327
0328   /// Generic method to calculates 4-momentum correction to be applied
0329   P4 calculateCorr(const P4& jet,
0330                    const TrackRefs&,
0331                    jpt::Efficiency&,
0332                    bool in_cone_at_vertex,
0333                    bool in_cone_at_calo_face,
0334                    double mass,
0335                    bool is_pion,
0336                    double mip);
0337
0338   /// Correction to be applied using tracking efficiency
0339   P4 pionEfficiency(const P4& jet, const jpt::Efficiency&, bool in_cone_at_calo_face);
0340
0341   /// Check corrected 4-momentum does not give negative scale
0342   double checkScale(const P4& jet, P4& corrected) const;
0343
0344   /// Get RECO muons
0345   bool getMuons(const edm::Event&, edm::Handle<RecoMuons>&) const;
0346
0347   bool getMuons(const edm::Event&, edm::Handle<pat::MuonCollection>&) const;
0348
0349   /// Get RECO electrons
0350   bool getElectrons(const edm::Event&, edm::Handle<RecoElectrons>&, edm::Handle<RecoElectronIds>&) const;
0351
0352   bool getElectrons(const edm::Event&, edm::Handle<pat::ElectronCollection>&) const;
0353
0354   /// Matches tracks to RECO muons
0355   bool matchMuons(TrackRefs::const_iterator&, const edm::Handle<RecoMuons>&) const;
0356   bool matchMuons(TrackRefs::const_iterator&, const edm::Handle<pat::MuonCollection>&) const;
0357
0358   /// Matches tracks to RECO electrons
0359   bool matchElectrons(TrackRefs::const_iterator&,
0360                       const edm::Handle<RecoElectrons>&,
0361                       const edm::Handle<RecoElectronIds>&) const;
0362   bool matchElectrons(TrackRefs::const_iterator&, const edm::Handle<pat::ElectronCollection>&) const;
0363
0364   /// Check on track quality
0365   bool failTrackQuality(TrackRefs::const_iterator&) const;
0366
0367   /// Find track in JetTracks collection
0368   bool findTrack(const jpt::JetTracks&, TrackRefs::const_iterator&, TrackRefs::iterator&) const;
0369
0370   /// Find track in MatchedTracks collections
0371   bool findTrack(const jpt::MatchedTracks& pions,
0372                  const jpt::MatchedTracks& muons,
0373                  const jpt::MatchedTracks& electrons,
0374                  TrackRefs::const_iterator&) const;
0375
0376   /// Determines if any tracks in cone at CaloFace
0377   bool tracksInCalo(const jpt::MatchedTracks& pions,
0378                     const jpt::MatchedTracks& muons,
0379                     const jpt::MatchedTracks& elecs) const;
0380
0381   /// Rebuild jet-track association
0382   void rebuildJta(const reco::Jet&, const JetTracksAssociations&, TrackRefs& included, TrackRefs& excluded) const;
0383
0384   /// Exclude jet-track association
0385   void excludeJta(const reco::Jet&, const JetTracksAssociations&, TrackRefs& included, const TrackRefs& excluded) const;
0386
0387   const jpt::Map& responseMap() const;
0388   const jpt::Map& efficiencyMap() const;
0389   const jpt::Map& leakageMap() const;
0390
0391   /// Default constructor
0392   JetPlusTrackCorrector() { ; }
0393
0394   // ---------- Protected member data ----------
0395
0396 protected:
0397   // Some general configuration
0398   bool verbose_;
0399   bool usePAT_;
0400   bool vectorial_;
0401   bool vecResponse_;
0402   bool useInConeTracks_;
0403   bool useOutOfConeTracks_;
0404   bool useOutOfVertexTracks_;
0405   bool usePions_;
0406   bool useEff_;
0407   bool useMuons_;
0408   bool useElecs_;
0409   bool useTrackQuality_;
0410
0411   // Jet-track association
0412   edm::InputTag jetTracksAtVertex_;
0413   edm::InputTag jetTracksAtCalo_;
0414   int jetSplitMerge_;
0415   edm::InputTag srcPVs_;
0416   double ptErrorQuality_;
0417   double dzVertexCut_;
0418   reco::Particle::Point vertex_;
0419
0420   // Muons and electrons
0421   edm::InputTag muons_;
0422   edm::InputTag electrons_;
0423   edm::InputTag electronIds_;
0424   edm::InputTag patmuons_;
0425   edm::InputTag patelectrons_;
0426
0427   // Filter tracks by quality
0428   reco::TrackBase::TrackQuality trackQuality_;
0429
0430   // Response and efficiency maps
0431   const jpt::Map response_;
0432   const jpt::Map efficiency_;
0433   const jpt::Map leakage_;
0434
0435   // Muon/Electron
0436   double muonPtmatch_;
0437   double muonEtamatch_;
0438   double muonPhimatch_;
0439   double electronDRmatch_;
0440
0441   // Mass
0442   double pionMass_;
0443   double muonMass_;
0444   double elecMass_;
0445
0446   // Jet-related
0447   double maxEta_;
0448   float theResponseOfChargedWithEff;
0449   float theResponseOfChargedWithoutEff;
0450   float theSumPtWithEff;
0451   float theSumPtWithoutEff;
0452   float theSumEnergyWithEff;
0453   float theSumEnergyWithoutEff;
0454   float theSumPtForBeta;
0455   jpt::Efficiency not_used{response_, efficiency_, leakage_};
0456
0457   // ---------- Private data members ----------
0458
0459 private:
0460   edm::EDGetTokenT<reco::JetTracksAssociation::Container> input_jetTracksAtVertex_token_;
0461   edm::EDGetTokenT<reco::JetTracksAssociation::Container> input_jetTracksAtCalo_token_;
0462   edm::EDGetTokenT<RecoMuons> input_reco_muons_token_;
0463   edm::EDGetTokenT<reco::VertexCollection> input_pvCollection_token_;
0464   edm::EDGetTokenT<RecoElectrons> input_reco_elecs_token_;
0465   edm::EDGetTokenT<RecoElectronIds> input_reco_elec_ids_token_;
0466   edm::EDGetTokenT<pat::MuonCollection> input_pat_muons_token_;
0467   edm::EDGetTokenT<pat::ElectronCollection> input_pat_elecs_token_;
0468 };
0469
0470 // ---------- Inline methods ----------
0471
0472 inline double JetPlusTrackCorrector::correction(const reco::Jet& fJet,
0473                                                 const reco::Jet& fJetcalo,
0474                                                 const edm::Event& event,
0475                                                 const edm::EventSetup& setup,
0476                                                 jpt::MatchedTracks& pions,
0477                                                 jpt::MatchedTracks& muons,
0478                                                 jpt::MatchedTracks& elecs,
0479                                                 bool& validMatches) {
0480   P4 not_used_for_scalar_correction;
0481   return correction(fJet, fJetcalo, event, setup, not_used_for_scalar_correction, pions, muons, elecs, validMatches);
0482 }
0483
0484 inline bool JetPlusTrackCorrector::eventRequired() const { return true; }
0485 inline bool JetPlusTrackCorrector::vectorialCorrection() const { return vectorial_; }
0486 inline bool JetPlusTrackCorrector::canCorrect(const reco::Jet& jet) const { return (fabs(jet.eta()) <= maxEta_); }
0487
0488 inline JetPlusTrackCorrector::P4 JetPlusTrackCorrector::pionCorrection(
0489     const P4& jet, const TrackRefs& pions, jpt::Efficiency& eff, bool in_cone_at_vertex, bool in_cone_at_calo_face) {
0490   return calculateCorr(jet, pions, eff, in_cone_at_vertex, in_cone_at_calo_face, pionMass_, true, -1.);
0491 }
0492
0493 inline JetPlusTrackCorrector::P4 JetPlusTrackCorrector::muonCorrection(const P4& jet,
0494                                                                        const TrackRefs& muons,
0495                                                                        bool in_cone_at_vertex,
0496                                                                        bool in_cone_at_calo_face) {
0497   return calculateCorr(jet, muons, not_used, in_cone_at_vertex, in_cone_at_calo_face, muonMass_, false, 2.);
0498 }
0499
0500 inline JetPlusTrackCorrector::P4 JetPlusTrackCorrector::elecCorrection(const P4& jet,
0501                                                                        const TrackRefs& elecs,
0502                                                                        bool in_cone_at_vertex,
0503                                                                        bool in_cone_at_calo_face) {
0504   return calculateCorr(jet, elecs, not_used, in_cone_at_vertex, in_cone_at_calo_face, elecMass_, false, 0.);
0505 }
0506
0507 inline double JetPlusTrackCorrector::checkScale(const P4& jet, P4& corrected) const {
0508   if (jet.energy() > 0. && (corrected.energy() / jet.energy()) < 0.) {
0509     corrected = jet;
0510   }
0511   return corrected.energy() / jet.energy();
0512 }
0513
0514 inline const jpt::Map& JetPlusTrackCorrector::responseMap() const { return response_; }
0515 inline const jpt::Map& JetPlusTrackCorrector::efficiencyMap() const { return efficiency_; }
0516 inline const jpt::Map& JetPlusTrackCorrector::leakageMap() const { return leakage_; }
0517
0518 #endif  // RecoJets_JetPlusTracks_JetPlusTrackCorrector_h