Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:35:01

0001 #include "JetPlusTrackCorrector.h"
0002 #include "DataFormats/Common/interface/Handle.h"
0003 #include "DataFormats/Common/interface/View.h"
0004 #include "DataFormats/Math/interface/deltaR.h"
0005 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0006 #include "FWCore/Framework/interface/ESHandle.h"
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "FWCore/Framework/interface/EventSetup.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0012 #include "DataFormats/Math/interface/deltaPhi.h"
0013 #include "DataFormats/Math/interface/deltaR.h"
0014 #include <fstream>
0015 #include <vector>
0016 #include <iomanip>
0017 #include <cmath>
0018 
0019 using namespace std;
0020 using namespace jpt;
0021 
0022 // -----------------------------------------------------------------------------
0023 //
0024 JetPlusTrackCorrector::JetPlusTrackCorrector(const edm::ParameterSet& pset, edm::ConsumesCollector&& iC)
0025     : verbose_(pset.getParameter<bool>("Verbose")),
0026       usePAT_(pset.getParameter<bool>("UsePAT")),
0027       vectorial_(pset.getParameter<bool>("VectorialCorrection")),
0028       vecResponse_(pset.getParameter<bool>("UseResponseInVecCorr")),
0029       useInConeTracks_(pset.getParameter<bool>("UseInConeTracks")),
0030       useOutOfConeTracks_(pset.getParameter<bool>("UseOutOfConeTracks")),
0031       useOutOfVertexTracks_(pset.getParameter<bool>("UseOutOfVertexTracks")),
0032       usePions_(pset.getParameter<bool>("UsePions")),
0033       useEff_(pset.getParameter<bool>("UseEfficiency")),
0034       useMuons_(pset.getParameter<bool>("UseMuons")),
0035       useElecs_(pset.getParameter<bool>("UseElectrons")),
0036       useTrackQuality_(pset.getParameter<bool>("UseTrackQuality")),
0037       jetTracksAtVertex_(pset.getParameter<edm::InputTag>("JetTracksAssociationAtVertex")),
0038       jetTracksAtCalo_(pset.getParameter<edm::InputTag>("JetTracksAssociationAtCaloFace")),
0039       jetSplitMerge_(pset.getParameter<int>("JetSplitMerge")),
0040       srcPVs_(pset.getParameter<edm::InputTag>("srcPVs")),
0041       ptErrorQuality_(pset.getParameter<double>("PtErrorQuality")),
0042       dzVertexCut_(pset.getParameter<double>("DzVertexCut")),
0043       muons_(pset.getParameter<edm::InputTag>("Muons")),
0044       electrons_(pset.getParameter<edm::InputTag>("Electrons")),
0045       electronIds_(pset.getParameter<edm::InputTag>("ElectronIds")),
0046       patmuons_(pset.getParameter<edm::InputTag>("PatMuons")),
0047       patelectrons_(pset.getParameter<edm::InputTag>("PatElectrons")),
0048       trackQuality_(reco::TrackBase::qualityByName(pset.getParameter<std::string>("TrackQuality"))),
0049       response_(Map(pset.getParameter<std::string>("ResponseMap"), verbose_)),
0050       efficiency_(Map(pset.getParameter<std::string>("EfficiencyMap"), verbose_)),
0051       leakage_(Map(pset.getParameter<std::string>("LeakageMap"), verbose_)),
0052       muonPtmatch_(pset.getParameter<double>("muonPtmatch")),
0053       muonEtamatch_(pset.getParameter<double>("muonEtamatch")),
0054       muonPhimatch_(pset.getParameter<double>("muonPhimatch")),
0055       electronDRmatch_(pset.getParameter<double>("electronDRmatch")),
0056       pionMass_(0.140),
0057       muonMass_(0.105),
0058       elecMass_(0.000511),
0059       maxEta_(pset.getParameter<double>("MaxJetEta")) {
0060   if (verbose_) {
0061     std::stringstream ss;
0062     ss << "[JetPlusTrackCorrector::" << __func__ << "] Configuration for JPT corrector: " << std::endl
0063        << " Particles" << std::endl
0064        << "  UsePions             : " << (usePions_ ? "true" : "false") << std::endl
0065        << "  UseMuons             : " << (useMuons_ ? "true" : "false") << std::endl
0066        << "  UseElecs             : " << (useElecs_ ? "true" : "false") << std::endl
0067        << " Corrections" << std::endl
0068        << "  UseInConeTracks      : " << (useInConeTracks_ ? "true" : "false") << std::endl
0069        << "  UseOutOfConeTracks   : " << (useOutOfConeTracks_ ? "true" : "false") << std::endl
0070        << "  UseOutOfVertexTracks : " << (useOutOfVertexTracks_ ? "true" : "false") << std::endl
0071        << "  ResponseMap          : " << pset.getParameter<std::string>("ResponseMap") << std::endl
0072        << " Efficiency" << std::endl
0073        << "  UsePionEfficiency    : " << (useEff_ ? "true" : "false") << std::endl
0074        << "  EfficiencyMap        : " << pset.getParameter<std::string>("EfficiencyMap") << std::endl
0075        << "  LeakageMap           : " << pset.getParameter<std::string>("LeakageMap") << std::endl
0076        << " Tracks" << std::endl
0077        << "  JetTracksAtVertex    : " << jetTracksAtVertex_ << std::endl
0078        << "  JetTracksAtCalo      : " << jetTracksAtCalo_ << std::endl
0079        << "  JetSplitMerge        : " << jetSplitMerge_ << std::endl
0080        << "  UseTrackQuality      : " << (useTrackQuality_ ? "true" : "false") << std::endl
0081        << " Collections" << std::endl
0082        << "  Muons                : " << muons_ << std::endl
0083        << "  Electrons            : " << electrons_ << std::endl
0084        << " Vectorial" << std::endl
0085        << "  UseTracksAndResponse : " << ((vectorial_ && vecResponse_) ? "true" : "false") << std::endl
0086        << "  UseTracksOnly        : " << ((vectorial_ && !vecResponse_) ? "true" : "false");
0087     edm::LogInfo("JetPlusTrackCorrector") << ss.str();
0088   }
0089 
0090   if (!useInConeTracks_ || !useOutOfConeTracks_ || !useOutOfVertexTracks_) {
0091     std::stringstream ss;
0092     ss << "[JetPlusTrackCorrector::" << __func__ << "]"
0093        << " You are using JPT algorithm in a non-standard way!" << std::endl
0094        << " UseInConeTracks      : " << (useInConeTracks_ ? "true" : "false") << std::endl
0095        << " UseOutOfConeTracks   : " << (useOutOfConeTracks_ ? "true" : "false") << std::endl
0096        << " UseOutOfVertexTracks : " << (useOutOfVertexTracks_ ? "true" : "false");
0097     edm::LogWarning("JetPlusTrackCorrector") << ss.str();
0098   }
0099 
0100   input_jetTracksAtVertex_token_ = iC.consumes<reco::JetTracksAssociation::Container>(jetTracksAtVertex_);
0101   input_jetTracksAtCalo_token_ = iC.consumes<reco::JetTracksAssociation::Container>(jetTracksAtCalo_);
0102   input_pvCollection_token_ = iC.consumes<reco::VertexCollection>(srcPVs_);
0103 
0104   input_reco_muons_token_ = iC.consumes<RecoMuons>(muons_);
0105   input_reco_elecs_token_ = iC.consumes<RecoElectrons>(electrons_);
0106   input_reco_elec_ids_token_ = iC.consumes<RecoElectronIds>(electronIds_);
0107   if (usePAT_) {
0108     input_pat_muons_token_ = iC.consumes<pat::MuonCollection>(patmuons_);
0109     input_pat_elecs_token_ = iC.consumes<pat::ElectronCollection>(patelectrons_);
0110   }
0111 }
0112 
0113 // -----------------------------------------------------------------------------
0114 //
0115 JetPlusTrackCorrector::~JetPlusTrackCorrector() { ; }
0116 // -----------------------------------------------------------------------------
0117 //
0118 
0119 double JetPlusTrackCorrector::correction(const reco::Jet& fJet,
0120                                          const reco::Jet& fJetcalo,
0121                                          const edm::Event& event,
0122                                          const edm::EventSetup& setup,
0123                                          const reco::TrackRefVector& tracksinvert,
0124                                          const reco::TrackRefVector& tracksincalo,
0125                                          P4& corrected,
0126                                          MatchedTracks& pions,
0127                                          MatchedTracks& muons,
0128                                          MatchedTracks& elecs) {
0129   double scale = 1.;
0130   corrected = fJet.p4();
0131   matchTracks(fJetcalo, event, setup, tracksinvert, tracksincalo, pions, muons, elecs);
0132   if (!usePAT_) {
0133     if (usePions_) {
0134       corrected += pionCorrection(fJet.p4(), pions);
0135     }
0136     if (useMuons_) {
0137       corrected += muonCorrection(fJet.p4(), muons);
0138     }
0139     if (useElecs_) {
0140       corrected += elecCorrection(fJet.p4(), elecs);
0141     }
0142   } else {
0143     corrected += fJetcalo.p4();
0144   }
0145   scale = checkScale(fJet.p4(), corrected);
0146   return scale;
0147 }
0148 
0149 // -----------------------------------------------------------------------------
0150 //
0151 double JetPlusTrackCorrector::correction(const reco::Jet& fJet,
0152                                          const reco::Jet& fJetcalo,
0153                                          const edm::Event& event,
0154                                          const edm::EventSetup& setup,
0155                                          P4& corrected,
0156                                          MatchedTracks& pions,
0157                                          MatchedTracks& muons,
0158                                          MatchedTracks& elecs,
0159                                          bool& validMatches) {
0160   theResponseOfChargedWithEff = 0.;
0161   theResponseOfChargedWithoutEff = 0.;
0162   theSumPtWithEff = 0.;
0163   theSumPtWithoutEff = 0.;
0164   theSumEnergyWithEff = 0.;
0165   theSumEnergyWithoutEff = 0.;
0166   theSumPtForBeta = 0.;
0167 
0168   // Corrected 4-momentum for jet
0169   corrected = fJet.p4();
0170 
0171   // Match tracks to different particle types
0172   validMatches = matchTracks(fJetcalo, event, setup, pions, muons, elecs);
0173   if (!validMatches) {
0174     return 1.;
0175   }
0176 
0177   // Check if jet can be JPT-corrected
0178   if (!canCorrect(fJet)) {
0179     return 1.;
0180   }
0181 
0182   // Debug
0183   if (verbose_) {
0184     edm::LogInfo("JetPlusTrackCorrector") << "[JetPlusTrackCorrector::" << __func__ << "]"
0185                                           << " Applying JPT corrections...";
0186   }
0187 
0188   // Pion corrections (both scalar and vectorial)
0189   if (usePions_) {
0190     corrected += pionCorrection(fJet.p4(), pions);
0191   }
0192 
0193   // Muon corrections (both scalar and vectorial)
0194   if (useMuons_) {
0195     corrected += muonCorrection(fJet.p4(), muons);
0196   }
0197 
0198   // Electrons corrections (both scalar and vectorial)
0199   if (useElecs_) {
0200     corrected += elecCorrection(fJet.p4(), elecs);
0201   }
0202 
0203   // Define jet direction using total 3-momentum of tracks (overrides above)
0204   if (vectorial_ && !vecResponse_) {
0205     if (fabs(corrected.eta()) < 2.) {
0206       corrected = jetDirFromTracks(corrected, pions, muons, elecs);
0207     }
0208   }
0209 
0210   // Check if corrected 4-momentum gives negative scale
0211   double scale = checkScale(fJet.p4(), corrected);
0212 
0213   // Debug
0214   if (verbose_) {
0215     std::stringstream ss;
0216     ss << "Total correction:" << std::endl
0217        << std::fixed << std::setprecision(6) << " Uncorrected (Px,Py,Pz,E)   : "
0218        << "(" << fJet.px() << "," << fJet.py() << "," << fJet.pz() << "," << fJet.energy() << ")" << std::endl
0219        << " Corrected (Px,Py,Pz,E)     : "
0220        << "(" << corrected.px() << "," << corrected.py() << "," << corrected.pz() << "," << corrected.energy() << ")"
0221        << std::endl
0222        << " Uncorrected (Pt,Eta,Phi,M) : "
0223        << "(" << fJet.pt() << "," << fJet.eta() << "," << fJet.phi() << "," << fJet.mass() << ")" << std::endl
0224        << " Corrected (Pt,Eta,Phi,M)   : "
0225        << "(" << corrected.pt() << "," << corrected.eta() << "," << corrected.phi() << "," << corrected.mass() << ")"
0226        << std::endl
0227        << " Scalar correction to E     : " << scale << std::endl
0228        << " Scalar correction to Et    : " << (fJet.et() > 0. ? corrected.Et() / fJet.et() : 1.);  // << std::endl
0229     edm::LogVerbatim("JetPlusTrackCorrector") << ss.str();
0230   }
0231 
0232   // Return energy correction
0233   return scale;
0234 }
0235 
0236 // -----------------------------------------------------------------------------
0237 //
0238 double JetPlusTrackCorrector::correction(const reco::Jet& jet) const {
0239   edm::LogError("JetPlusTrackCorrector") << "JetPlusTrackCorrector can be run on entire event only";
0240   return 1.;
0241 }
0242 
0243 // -----------------------------------------------------------------------------
0244 //
0245 double JetPlusTrackCorrector::correction(const reco::Particle::LorentzVector& jet) const {
0246   edm::LogError("JetPlusTrackCorrector") << "JetPlusTrackCorrector can be run on entire event only";
0247   return 1.;
0248 }
0249 
0250 // -----------------------------------------------------------------------------
0251 //
0252 void JetPlusTrackCorrector::matchTracks(const reco::Jet& fJet,
0253                                         const edm::Event& event,
0254                                         const edm::EventSetup& setup,
0255                                         const reco::TrackRefVector& tracksinvert,
0256                                         const reco::TrackRefVector& tracksincalo,
0257                                         jpt::MatchedTracks& pions,
0258                                         jpt::MatchedTracks& muons,
0259                                         jpt::MatchedTracks& elecs) {
0260   JetTracks jet_tracks;
0261   jet_tracks.vertex_ = tracksinvert;
0262   jet_tracks.caloFace_ = tracksincalo;
0263   matchTracks(jet_tracks, event, pions, muons, elecs);
0264 
0265   return;
0266 }
0267 
0268 bool JetPlusTrackCorrector::matchTracks(const reco::Jet& fJet,
0269                                         const edm::Event& event,
0270                                         const edm::EventSetup& setup,
0271                                         jpt::MatchedTracks& pions,
0272                                         jpt::MatchedTracks& muons,
0273                                         jpt::MatchedTracks& elecs) {
0274   // Associate tracks to jet at both the Vertex and CaloFace
0275   JetTracks jet_tracks;
0276   bool ok = jetTrackAssociation(fJet, event, setup, jet_tracks);
0277 
0278   if (!ok) {
0279     return false;
0280   }
0281 
0282   // Track collections propagated to Vertex and CaloFace for "pions", muons and electrons
0283   matchTracks(jet_tracks, event, pions, muons, elecs);
0284 
0285   // Debug
0286   if (verbose_) {
0287     std::stringstream ss;
0288     ss << "Number of tracks:" << std::endl
0289        << " In-cone at Vertex   : " << jet_tracks.vertex_.size() << std::endl
0290        << " In-cone at CaloFace : " << jet_tracks.caloFace_.size();
0291     edm::LogVerbatim("JetPlusTrackCorrector") << ss.str();
0292   }
0293 
0294   return true;
0295 }
0296 
0297 // -----------------------------------------------------------------------------
0298 //
0299 bool JetPlusTrackCorrector::jetTrackAssociation(const reco::Jet& fJet,
0300                                                 const edm::Event& event,
0301                                                 const edm::EventSetup& setup,
0302                                                 JetTracks& trks) const {
0303   // Some init
0304   trks.clear();
0305 
0306   // Check if labels are given
0307   if (!jetTracksAtVertex_.label().empty() && !jetTracksAtCalo_.label().empty()) {
0308     return jtaUsingEventData(fJet, event, trks);
0309   } else {
0310     edm::LogWarning("PatJPTCorrector") << "[JetPlusTrackCorrector::" << __func__ << "]"
0311                                        << " Empty label for the reco::JetTracksAssociation::Containers" << std::endl
0312                                        << " InputTag for JTA \"at vertex\"    (label:instance:process) \""
0313                                        << jetTracksAtVertex_.label() << ":" << jetTracksAtVertex_.instance() << ":"
0314                                        << jetTracksAtVertex_.process() << "\"" << std::endl
0315                                        << " InputTag for JTA \"at calo face\" (label:instance:process) \""
0316                                        << jetTracksAtCalo_.label() << ":" << jetTracksAtCalo_.instance() << ":"
0317                                        << jetTracksAtCalo_.process() << "\"";
0318     return false;
0319   }
0320 }
0321 
0322 // -----------------------------------------------------------------------------
0323 //
0324 bool JetPlusTrackCorrector::jtaUsingEventData(const reco::Jet& fJet, const edm::Event& event, JetTracks& trks) const {
0325   // Get Jet-track association at Vertex
0326   edm::Handle<reco::JetTracksAssociation::Container> jetTracksAtVertex;
0327   event.getByToken(input_jetTracksAtVertex_token_, jetTracksAtVertex);
0328 
0329   if (!jetTracksAtVertex.isValid() || jetTracksAtVertex.failedToGet()) {
0330     if (verbose_ && edm::isDebugEnabled()) {
0331       edm::LogWarning("JetPlusTrackCorrector")
0332           << "[JetPlusTrackCorrector::" << __func__ << "]"
0333           << " Invalid handle to reco::JetTracksAssociation::Container (for Vertex)"
0334           << " with InputTag (label:instance:process) \"" << jetTracksAtVertex_.label() << ":"
0335           << jetTracksAtVertex_.instance() << ":" << jetTracksAtVertex_.process() << "\"";
0336     }
0337     return false;
0338   }
0339 
0340   // Retrieve jet-tracks association for given jet
0341   const reco::JetTracksAssociation::Container jtV = *(jetTracksAtVertex.product());
0342   TrackRefs excluded;
0343   TrackRefs relocate;
0344   if (jetSplitMerge_ < 0) {
0345     trks.vertex_ = reco::JetTracksAssociation::getValue(jtV, fJet);
0346   } else {
0347     rebuildJta(fJet, jtV, trks.vertex_, excluded);
0348     rebuildJta(fJet, jtV, relocate, excluded);
0349     trks.vertex_ = relocate;
0350   }
0351 
0352   // Check if any tracks are associated to jet at vertex
0353   if (trks.vertex_.empty()) {
0354     return false;
0355   }
0356 
0357   // Get Jet-track association at Calo
0358   edm::Handle<reco::JetTracksAssociation::Container> jetTracksAtCalo;
0359   event.getByToken(input_jetTracksAtCalo_token_, jetTracksAtCalo);
0360 
0361   if (!jetTracksAtCalo.isValid() || jetTracksAtCalo.failedToGet()) {
0362     if (verbose_ && edm::isDebugEnabled()) {
0363       edm::LogWarning("JetPlusTrackCorrector")
0364           << "[JetPlusTrackCorrector::" << __func__ << "]"
0365           << " Invalid handle to reco::JetTracksAssociation::Container (for CaloFace)"
0366           << " with InputTag (label:instance:process) \"" << jetTracksAtCalo_.label() << ":"
0367           << jetTracksAtCalo_.instance() << ":" << jetTracksAtCalo_.process() << "\"";
0368     }
0369     return false;
0370   }
0371 
0372   // Retrieve jet-tracks association for given jet
0373   const reco::JetTracksAssociation::Container jtC = *(jetTracksAtCalo.product());
0374   if (jetSplitMerge_ < 0) {
0375     trks.caloFace_ = reco::JetTracksAssociation::getValue(jtC, fJet);
0376   } else {
0377     excludeJta(fJet, jtC, trks.caloFace_, excluded);
0378   }
0379 
0380   return true;
0381 }
0382 
0383 // -----------------------------------------------------------------------------
0384 //
0385 bool JetPlusTrackCorrector::getMuons(const edm::Event& event, edm::Handle<RecoMuons>& reco_muons) const {
0386   event.getByToken(input_reco_muons_token_, reco_muons);
0387   return true;
0388 }
0389 
0390 bool JetPlusTrackCorrector::getMuons(const edm::Event& event, edm::Handle<pat::MuonCollection>& pat_muons) const {
0391   event.getByToken(input_pat_muons_token_, pat_muons);
0392   return true;
0393 }
0394 
0395 // -----------------------------------------------------------------------------
0396 //
0397 void JetPlusTrackCorrector::matchTracks(const JetTracks& jet_tracks,
0398                                         const edm::Event& event,
0399                                         MatchedTracks& pions,
0400                                         MatchedTracks& muons,
0401                                         MatchedTracks& elecs) {
0402   // Some init
0403   pions.clear();
0404   muons.clear();
0405   elecs.clear();
0406 
0407   // Need vertex for track cleaning
0408 
0409   vertex_ = reco::Particle::Point(0, 0, 0);
0410   edm::Handle<reco::VertexCollection> pvCollection;
0411   event.getByToken(input_pvCollection_token_, pvCollection);
0412   if (pvCollection.isValid() && !pvCollection->empty())
0413     vertex_ = pvCollection->begin()->position();
0414 
0415   // Get RECO muons
0416   edm::Handle<RecoMuons> reco_muons;
0417   edm::Handle<pat::MuonCollection> pat_muons;
0418   bool found_reco_muons = true;
0419   bool found_pat_muons = true;
0420   if (useMuons_) {
0421     if (!usePAT_) {
0422       getMuons(event, reco_muons);
0423     } else {
0424       getMuons(event, pat_muons);
0425       found_reco_muons = false;
0426     }
0427   }
0428 
0429   // Get RECO electrons and their ids
0430   edm::Handle<RecoElectrons> reco_elecs;
0431   edm::Handle<pat::ElectronCollection> pat_elecs;
0432   edm::Handle<RecoElectronIds> reco_elec_ids;
0433   bool found_reco_elecs = true;
0434   bool found_pat_elecs = true;
0435   if (useElecs_) {
0436     if (!usePAT_) {
0437       getElectrons(event, reco_elecs, reco_elec_ids);
0438     } else {
0439       getElectrons(event, pat_elecs);
0440       found_reco_elecs = false;
0441     }
0442   }
0443 
0444   // Check RECO products found
0445   if (!found_reco_muons || !found_reco_elecs) {
0446     if (!found_pat_muons || !found_pat_elecs) {
0447       edm::LogError("JetPlusTrackCorrector") << "[JetPlusTrackCorrector::" << __func__ << "]"
0448                                              << " Unable to access RECO collections for muons and electrons";
0449       return;
0450     }
0451   }
0452 
0453   // Identify pions/muons/electrons that are "in/in" and "in/out"
0454   {
0455     TrackRefs::const_iterator itrk = jet_tracks.vertex_.begin();
0456     TrackRefs::const_iterator jtrk = jet_tracks.vertex_.end();
0457 
0458     double theSumPtForBetaOld = theSumPtForBeta;
0459 
0460     for (; itrk != jtrk; ++itrk) {
0461       if (useTrackQuality_ && (*itrk)->quality(trackQuality_) && theSumPtForBetaOld <= 0.)
0462         theSumPtForBeta += (**itrk).pt();
0463       //
0464       // Track either belongs to PV or do not belong to any vertex
0465 
0466       const reco::TrackBaseRef ttr1(*itrk);
0467 
0468       int numpv = 0;
0469 
0470       int itrack_belong = -1;
0471 
0472       for (reco::VertexCollection::const_iterator iv = pvCollection->begin(); iv != pvCollection->end(); iv++) {
0473         numpv++;
0474         std::vector<reco::TrackBaseRef>::const_iterator rr = find((*iv).tracks_begin(), (*iv).tracks_end(), ttr1);
0475         if (rr != (*iv).tracks_end()) {
0476           itrack_belong++;
0477           break;
0478         }
0479 
0480       }  // all vertices
0481 
0482       if (numpv > 1 && itrack_belong == 0) {
0483         continue;
0484       }
0485 
0486       if (failTrackQuality(itrk)) {
0487         continue;
0488       }
0489 
0490       TrackRefs::iterator it = jet_tracks.caloFace_.end();
0491       bool found = findTrack(jet_tracks, itrk, it);
0492       bool muaccept = false;
0493       bool eleaccept = false;
0494       if (found_reco_muons)
0495         muaccept = matchMuons(itrk, reco_muons);
0496       else if (found_pat_muons) {
0497         muaccept = matchMuons(itrk, pat_muons);
0498       }
0499       if (found_reco_elecs)
0500         eleaccept = matchElectrons(itrk, reco_elecs, reco_elec_ids);
0501       else if (found_pat_elecs) {
0502         eleaccept = matchElectrons(itrk, pat_elecs);
0503       }
0504 
0505       bool is_muon = useMuons_ && muaccept;
0506       bool is_ele = useElecs_ && eleaccept;
0507 
0508       if (found) {
0509         if (is_muon) {
0510           muons.inVertexInCalo_.push_back(*it);
0511         } else if (is_ele) {
0512           elecs.inVertexInCalo_.push_back(*it);
0513         } else {
0514           pions.inVertexInCalo_.push_back(*it);
0515         }
0516       } else {
0517         if (is_muon) {
0518           muons.inVertexOutOfCalo_.push_back(*itrk);
0519         } else if (is_ele) {
0520           elecs.inVertexOutOfCalo_.push_back(*itrk);
0521         } else {
0522           pions.inVertexOutOfCalo_.push_back(*itrk);
0523         }
0524       }
0525     }
0526   }
0527 
0528   // Identify pions/muons/electrons that are "out/in"
0529   {
0530     TrackRefs::iterator itrk = jet_tracks.caloFace_.begin();
0531     TrackRefs::iterator jtrk = jet_tracks.caloFace_.end();
0532     for (; itrk != jtrk; ++itrk) {
0533       if (failTrackQuality(itrk)) {
0534         continue;
0535       }
0536 
0537       if (!tracksInCalo(pions, muons, elecs)) {
0538         continue;
0539       }
0540 
0541       bool found = findTrack(pions, muons, elecs, itrk);
0542 
0543       if (!found) {
0544         bool muaccept = false;
0545         bool eleaccept = false;
0546         if (found_reco_muons)
0547           muaccept = matchMuons(itrk, reco_muons);
0548         else if (found_pat_muons) {
0549           muaccept = matchMuons(itrk, pat_muons);
0550         }
0551         if (found_reco_elecs)
0552           eleaccept = matchElectrons(itrk, reco_elecs, reco_elec_ids);
0553         else if (found_pat_elecs) {
0554           eleaccept = matchElectrons(itrk, pat_elecs);
0555         }
0556 
0557         bool is_muon = useMuons_ && muaccept;
0558         bool is_ele = useElecs_ && eleaccept;
0559 
0560         if (is_muon) {
0561           muons.outOfVertexInCalo_.push_back(*itrk);
0562         } else if (is_ele) {
0563           elecs.outOfVertexInCalo_.push_back(*itrk);
0564         } else {
0565           pions.outOfVertexInCalo_.push_back(*itrk);
0566         }
0567       }
0568     }
0569   }
0570 
0571   if (verbose_ && edm::isDebugEnabled()) {
0572     std::stringstream ss;
0573     ss << "[JetPlusTrackCorrector::" << __func__ << "] Number of tracks:" << std::endl
0574        << " In-cone at Vertex and in-cone at CaloFace:" << std::endl
0575        << "  Pions      : " << pions.inVertexInCalo_.size() << std::endl
0576        << "  Muons      : " << muons.inVertexInCalo_.size() << std::endl
0577        << "  Electrons  : " << elecs.inVertexInCalo_.size() << std::endl
0578        << " In-cone at Vertex and out-of-cone at CaloFace:" << std::endl
0579        << "  Pions      : " << pions.inVertexOutOfCalo_.size() << std::endl
0580        << "  Muons      : " << muons.inVertexOutOfCalo_.size() << std::endl
0581        << "  Electrons  : " << elecs.inVertexOutOfCalo_.size() << std::endl
0582        << " Out-of-cone at Vertex and in-cone at CaloFace:" << std::endl
0583        << "  Pions      : " << pions.outOfVertexInCalo_.size() << std::endl
0584        << "  Muons      : " << muons.outOfVertexInCalo_.size() << std::endl
0585        << "  Electrons  : " << elecs.outOfVertexInCalo_.size() << std::endl;
0586     LogTrace("JetPlusTrackCorrector") << ss.str();
0587   }
0588 }
0589 
0590 // -----------------------------------------------------------------------------
0591 //
0592 bool JetPlusTrackCorrector::getElectrons(const edm::Event& event,
0593                                          edm::Handle<RecoElectrons>& reco_elecs,
0594                                          edm::Handle<RecoElectronIds>& reco_elec_ids) const {
0595   event.getByToken(input_reco_elecs_token_, reco_elecs);
0596   event.getByToken(input_reco_elec_ids_token_, reco_elec_ids);
0597   return true;
0598 }
0599 
0600 bool JetPlusTrackCorrector::getElectrons(const edm::Event& event,
0601                                          edm::Handle<pat::ElectronCollection>& pat_elecs) const {
0602   event.getByToken(input_pat_elecs_token_, pat_elecs);
0603   return true;
0604 }
0605 
0606 // -----------------------------------------------------------------------------
0607 //
0608 bool JetPlusTrackCorrector::failTrackQuality(TrackRefs::const_iterator& itrk) const {
0609   bool retcode = false;
0610 
0611   if (useTrackQuality_ && !(*itrk)->quality(trackQuality_)) {
0612     retcode = true;
0613     return retcode;
0614   }
0615   if (((*itrk)->ptError() / (*itrk)->pt()) > ptErrorQuality_) {
0616     retcode = true;
0617     return retcode;
0618   }
0619   if (fabs((*itrk)->dz(vertex_)) > dzVertexCut_) {
0620     retcode = true;
0621     return retcode;
0622   }
0623 
0624   return retcode;
0625 }
0626 
0627 // -----------------------------------------------------------------------------
0628 //
0629 bool JetPlusTrackCorrector::findTrack(const JetTracks& jet_tracks,
0630                                       TrackRefs::const_iterator& itrk,
0631                                       TrackRefs::iterator& it) const {
0632   it = find(jet_tracks.caloFace_.begin(), jet_tracks.caloFace_.end(), *itrk);
0633   if (it != jet_tracks.caloFace_.end()) {
0634     return true;
0635   } else {
0636     return false;
0637   }
0638 }
0639 
0640 // -----------------------------------------------------------------------------
0641 //
0642 bool JetPlusTrackCorrector::findTrack(const MatchedTracks& pions,
0643                                       const MatchedTracks& muons,
0644                                       const MatchedTracks& elecs,
0645                                       TrackRefs::const_iterator& itrk) const {
0646   TrackRefs::iterator ip = find(pions.inVertexInCalo_.begin(), pions.inVertexInCalo_.end(), *itrk);
0647   TrackRefs::iterator im = find(muons.inVertexInCalo_.begin(), muons.inVertexInCalo_.end(), *itrk);
0648   TrackRefs::iterator ie = find(elecs.inVertexInCalo_.begin(), elecs.inVertexInCalo_.end(), *itrk);
0649   if (ip == pions.inVertexInCalo_.end() && im == muons.inVertexInCalo_.end() && ie == elecs.inVertexInCalo_.end()) {
0650     return false;
0651   } else {
0652     return true;
0653   }
0654 }
0655 
0656 // -----------------------------------------------------------------------------
0657 //
0658 bool JetPlusTrackCorrector::tracksInCalo(const MatchedTracks& pions,
0659                                          const MatchedTracks& muons,
0660                                          const MatchedTracks& elecs) const {
0661   if (!pions.inVertexInCalo_.empty() || !muons.inVertexInCalo_.empty() || !elecs.inVertexInCalo_.empty()) {
0662     return true;
0663   } else {
0664     return false;
0665   }
0666 }
0667 
0668 // -----------------------------------------------------------------------------
0669 //
0670 JetPlusTrackCorrector::P4 JetPlusTrackCorrector::pionCorrection(const P4& jet, const MatchedTracks& pions) {
0671   P4 corr_pions;
0672 
0673   // In-cone
0674 
0675   P4 corr_pions_in_cone;
0676   P4 corr_pions_eff_in_cone;
0677   Efficiency in_cone(responseMap(), efficiencyMap(), leakageMap());
0678 
0679   if (useInConeTracks_) {
0680     corr_pions_in_cone = pionCorrection(jet, pions.inVertexInCalo_, in_cone, true, true);
0681     corr_pions += corr_pions_in_cone;
0682     if (useEff_) {
0683       corr_pions_eff_in_cone = pionEfficiency(jet, in_cone, true);
0684       corr_pions += corr_pions_eff_in_cone;
0685     }
0686   }
0687 
0688   // Out-of-cone
0689 
0690   P4 corr_pions_out_of_cone;
0691   P4 corr_pions_eff_out_of_cone;
0692   Efficiency out_of_cone(responseMap(), efficiencyMap(), leakageMap());
0693 
0694   if (useOutOfConeTracks_) {
0695     corr_pions_out_of_cone = pionCorrection(jet, pions.inVertexOutOfCalo_, out_of_cone, true, false);
0696     corr_pions += corr_pions_out_of_cone;
0697     if (useEff_) {
0698       corr_pions_eff_out_of_cone = pionEfficiency(jet, out_of_cone, false);
0699       corr_pions += corr_pions_eff_out_of_cone;
0700     }
0701   }
0702 
0703   // Out-of-vertex
0704 
0705   P4 corr_pions_out_of_vertex;
0706   Efficiency not_used(responseMap(), efficiencyMap(), leakageMap());
0707 
0708   if (useOutOfVertexTracks_) {
0709     corr_pions_out_of_vertex = pionCorrection(jet, pions.outOfVertexInCalo_, not_used, false, true);
0710     corr_pions += corr_pions_out_of_vertex;
0711   }
0712 
0713   if (verbose_) {
0714     std::stringstream ss;
0715     ss << " Pion corrections:" << std::endl
0716        << "  In/In      : "
0717        << "(" << pions.inVertexInCalo_.size() << ") " << corr_pions_in_cone.energy() << std::endl
0718        << "  In/Out     : "
0719        << "(" << pions.inVertexOutOfCalo_.size() << ") " << corr_pions_out_of_cone.energy() << std::endl
0720        << "  Out/In     : "
0721        << "(" << pions.outOfVertexInCalo_.size() << ") " << corr_pions_out_of_vertex.energy() << std::endl;
0722     if (useEff_) {
0723       ss << " Pion efficiency corrections:" << std::endl
0724          << "  In/In      : "
0725          << "(" << pions.inVertexInCalo_.size() << ") " << corr_pions_eff_in_cone.energy() << std::endl
0726          << "  In/Out     : "
0727          << "(" << pions.inVertexOutOfCalo_.size() << ") " << corr_pions_eff_out_of_cone.energy();
0728     }
0729     edm::LogVerbatim("JetPlusTrackCorrector") << ss.str();
0730   }
0731 
0732   return corr_pions;
0733 }
0734 
0735 // -----------------------------------------------------------------------------
0736 //
0737 JetPlusTrackCorrector::P4 JetPlusTrackCorrector::muonCorrection(const P4& jet, const MatchedTracks& muons) {
0738   P4 corr_muons;
0739 
0740   P4 corr_muons_in_cone;
0741   P4 corr_muons_out_of_cone;
0742   P4 corr_muons_out_of_vertex;
0743 
0744   if (useInConeTracks_) {
0745     corr_muons_in_cone = muonCorrection(jet, muons.inVertexInCalo_, true, true);
0746     corr_muons += corr_muons_in_cone;
0747   }
0748 
0749   if (useOutOfConeTracks_) {
0750     corr_muons_out_of_cone = muonCorrection(jet, muons.inVertexOutOfCalo_, true, false);
0751     corr_muons += corr_muons_out_of_cone;
0752   }
0753 
0754   if (useOutOfVertexTracks_) {
0755     corr_muons_out_of_vertex = muonCorrection(jet, muons.outOfVertexInCalo_, false, true);
0756     corr_muons += corr_muons_out_of_vertex;
0757   }
0758 
0759   if (verbose_) {
0760     std::stringstream ss;
0761     ss << " Muon corrections:" << std::endl
0762        << "  In/In      : "
0763        << "(" << muons.inVertexInCalo_.size() << ") " << corr_muons_in_cone.energy() << std::endl
0764        << "  In/Out     : "
0765        << "(" << muons.inVertexOutOfCalo_.size() << ") " << corr_muons_out_of_cone.energy() << std::endl
0766        << "  Out/In     : "
0767        << "(" << muons.outOfVertexInCalo_.size() << ") " << corr_muons_out_of_vertex.energy();
0768     edm::LogVerbatim("JetPlusTrackCorrector") << ss.str();
0769   }
0770 
0771   return corr_muons;
0772 }
0773 
0774 // -----------------------------------------------------------------------------
0775 //
0776 JetPlusTrackCorrector::P4 JetPlusTrackCorrector::elecCorrection(const P4& jet, const MatchedTracks& elecs) const {
0777   P4 null;  //@@ null 4-momentum
0778 
0779   if (verbose_) {
0780     std::stringstream ss;
0781     ss << " Electron corrections:" << std::endl
0782        << "  In/In      : "
0783        << "(" << elecs.inVertexInCalo_.size() << ") " << null.energy() << std::endl
0784        << "  In/Out     : "
0785        << "(" << elecs.inVertexOutOfCalo_.size() << ") " << null.energy() << std::endl
0786        << "  Out/In     : "
0787        << "(" << elecs.outOfVertexInCalo_.size() << ") " << null.energy();
0788     edm::LogVerbatim("JetPlusTrackCorrector") << ss.str();
0789   }
0790 
0791   return null;  //@@ to be implemented
0792 }
0793 
0794 // -----------------------------------------------------------------------------
0795 //
0796 JetPlusTrackCorrector::P4 JetPlusTrackCorrector::jetDirFromTracks(const P4& corrected,
0797                                                                   const MatchedTracks& pions,
0798                                                                   const MatchedTracks& muons,
0799                                                                   const MatchedTracks& elecs) const {
0800   // Correction to be applied to jet 4-momentum
0801   P4 corr;
0802 
0803   //  bool tracks_present = false;
0804   bool tracks_present_inin = false;
0805 
0806   // Correct using pions in-cone at vertex
0807 
0808   if (!pions.inVertexInCalo_.empty()) {
0809     TrackRefs::iterator itrk = pions.inVertexInCalo_.begin();
0810     TrackRefs::iterator jtrk = pions.inVertexInCalo_.end();
0811     for (; itrk != jtrk; ++itrk) {
0812       corr += PtEtaPhiM((*itrk)->pt(), (*itrk)->eta(), (*itrk)->phi(), 0.);
0813       tracks_present_inin = true;
0814     }
0815   }
0816 
0817   if (!pions.inVertexOutOfCalo_.empty()) {
0818     TrackRefs::iterator itrk = pions.inVertexOutOfCalo_.begin();
0819     TrackRefs::iterator jtrk = pions.inVertexOutOfCalo_.end();
0820     for (; itrk != jtrk; ++itrk) {
0821       corr += PtEtaPhiM((*itrk)->pt(), (*itrk)->eta(), (*itrk)->phi(), 0.);
0822     }
0823   }
0824 
0825   // Correct using muons in-cone at vertex
0826 
0827   if (!muons.inVertexInCalo_.empty()) {
0828     TrackRefs::iterator itrk = muons.inVertexInCalo_.begin();
0829     TrackRefs::iterator jtrk = muons.inVertexInCalo_.end();
0830     for (; itrk != jtrk; ++itrk) {
0831       corr += PtEtaPhiM((*itrk)->pt(), (*itrk)->eta(), (*itrk)->phi(), 0.);
0832       //      tracks_present = true;
0833     }
0834   }
0835 
0836   if (!muons.inVertexOutOfCalo_.empty()) {
0837     TrackRefs::iterator itrk = muons.inVertexOutOfCalo_.begin();
0838     TrackRefs::iterator jtrk = muons.inVertexOutOfCalo_.end();
0839     for (; itrk != jtrk; ++itrk) {
0840       corr += PtEtaPhiM((*itrk)->pt(), (*itrk)->eta(), (*itrk)->phi(), 0.);
0841       //      tracks_present = true;
0842     }
0843   }
0844 
0845   // Correct using electrons in-cone at vertex
0846 
0847   if (!elecs.inVertexInCalo_.empty()) {
0848     TrackRefs::iterator itrk = elecs.inVertexInCalo_.begin();
0849     TrackRefs::iterator jtrk = elecs.inVertexInCalo_.end();
0850     for (; itrk != jtrk; ++itrk) {
0851       corr += PtEtaPhiM((*itrk)->pt(), (*itrk)->eta(), (*itrk)->phi(), 0.);
0852       //      tracks_present = true;
0853     }
0854   }
0855 
0856   if (!elecs.inVertexOutOfCalo_.empty()) {
0857     TrackRefs::iterator itrk = elecs.inVertexOutOfCalo_.begin();
0858     TrackRefs::iterator jtrk = elecs.inVertexOutOfCalo_.end();
0859     for (; itrk != jtrk; ++itrk) {
0860       corr += PtEtaPhiM((*itrk)->pt(), (*itrk)->eta(), (*itrk)->phi(), 0.);
0861       //      tracks_present = true;
0862     }
0863   }
0864 
0865   // Adjust direction if in cone tracks are present
0866 
0867   if (!tracks_present_inin) {
0868     corr = corrected;
0869   } else {
0870     corr *= (corr.P() > 0. ? corrected.P() / corr.P() : 1.);
0871     corr = P4(corr.px(), corr.py(), corr.pz(), corrected.energy());
0872   }
0873 
0874   return corr;
0875 }
0876 
0877 // -----------------------------------------------------------------------------
0878 //
0879 JetPlusTrackCorrector::P4 JetPlusTrackCorrector::calculateCorr(const P4& jet,
0880                                                                const TrackRefs& tracks,
0881                                                                jpt::Efficiency& eff,
0882                                                                bool in_cone_at_vertex,
0883                                                                bool in_cone_at_calo_face,
0884                                                                double mass,
0885                                                                bool is_pion,
0886                                                                double mip) {
0887   // Correction to be applied to jet 4-momentum
0888   P4 correction;
0889 
0890   // Reset efficiency container
0891   eff.reset();
0892 
0893   double theSumResp = 0;
0894   double theSumPt = 0;
0895   double theSumEnergy = 0;
0896 
0897   // Iterate through tracks
0898   if (!tracks.empty()) {
0899     TrackRefs::iterator itrk = tracks.begin();
0900     TrackRefs::iterator jtrk = tracks.end();
0901 
0902     for (; itrk != jtrk; ++itrk) {
0903       // Ignore high-pt tracks (only when in-cone and not a mip)
0904       if (in_cone_at_calo_face && is_pion && (*itrk)->pt() >= 50.) {
0905         continue;
0906       }
0907 
0908       // Inner track 4-momentum
0909       P4 inner;
0910       if (vectorial_ && vecResponse_) {
0911         inner = PtEtaPhiM((*itrk)->pt(), (*itrk)->eta(), (*itrk)->phi(), mass);
0912       } else {
0913         double energy = sqrt((*itrk)->px() * (*itrk)->px() + (*itrk)->py() * (*itrk)->py() +
0914                              (*itrk)->pz() * (*itrk)->pz() + mass * mass);
0915         inner = (jet.energy() > 0. ? energy / jet.energy() : 1.) * jet;
0916       }
0917 
0918       // Add track momentum (if in-cone at vertex)
0919       if (in_cone_at_vertex) {
0920         correction += inner;
0921       }
0922 
0923       // Find appropriate eta/pt bin for given track
0924       double eta = fabs((*itrk)->eta());
0925       double pt = fabs((*itrk)->pt());
0926       uint32_t ieta = response_.etaBin(eta);
0927       uint32_t ipt = response_.ptBin(pt);
0928 
0929       // Check bins (not for mips)
0930       if (is_pion && (ieta == response_.nEtaBins() || ipt == response_.nPtBins())) {
0931         continue;
0932       }
0933 
0934       // Outer track 4-momentum
0935       P4 outer;
0936       if (in_cone_at_calo_face) {
0937         if (vectorial_ && vecResponse_) {
0938           // Build 4-momentum from outer track (SHOULD USE IMPACT POINT?!)
0939           double outer_pt = (*itrk)->pt();
0940           double outer_eta = (*itrk)->eta();
0941           double outer_phi = (*itrk)->phi();
0942           if ((*itrk)->extra().isNonnull()) {
0943             outer_pt = (*itrk)->pt();
0944             outer_eta = (*itrk)->outerPosition().eta();  //@@ outerMomentum().eta()
0945             outer_phi = (*itrk)->outerPosition().phi();  //@@ outerMomentum().phi()
0946           }
0947           outer = PtEtaPhiM(outer_pt, outer_eta, outer_phi, mass);
0948           // Check if mip or not
0949           if (!is_pion) {
0950             outer *= (outer.energy() > 0. ? mip / outer.energy() : 1.);
0951           }  //@@ Scale to mip energy
0952           else {
0953             outer *= (outer.energy() > 0. ? inner.energy() / outer.energy() : 1.);
0954           }  //@@ Scale to inner track energy
0955         } else {
0956           // Check if mip or not
0957           if (!is_pion) {
0958             outer = (jet.energy() > 0. ? mip / jet.energy() : 1.) * jet;
0959           }  //@@ Jet 4-mom scaled by mip energy
0960           else {
0961             outer = inner;
0962           }  //@@ Set to inner track 4-momentum
0963         }
0964         if (is_pion) {
0965           outer *= response_.value(ieta, ipt);
0966         }                     //@@ Scale by pion response
0967         correction -= outer;  //@@ Subtract
0968 
0969         // Calculate the sum of responses
0970         theSumResp += response_.value(ieta, ipt);
0971       }
0972 
0973       // Calculate the sum of pt and energies
0974       theSumPt += inner.pt();
0975       theSumEnergy += inner.energy();
0976 
0977       // Record inner track energy for pion efficiency correction
0978       if (is_pion) {
0979         eff.addE(ieta, ipt, inner.energy());
0980       }
0981 
0982       // Debug
0983       if (verbose_ && edm::isDebugEnabled()) {
0984         std::stringstream temp;
0985         temp << " Response[" << ieta << "," << ipt << "]";
0986         std::stringstream ss;
0987         ss << "[JetPlusTrackCorrector::" << __func__ << "]" << std::endl
0988            << " Track eta / pt    : " << eta << " / " << pt << std::endl
0989            << temp.str() << std::setw(21 - temp.str().size()) << " : " << response_.value(ieta, ipt) << std::endl
0990            << " Track momentum added : " << inner.energy() << std::endl
0991            << " Response subtracted  : " << outer.energy() << std::endl
0992            << " Energy correction    : " << correction.energy();
0993         LogDebug("JetPlusTrackCorrector") << ss.str();
0994       }
0995 
0996     }  // loop through tracks
0997   }    // ntracks != 0
0998 
0999   if (in_cone_at_vertex) {
1000     theResponseOfChargedWithEff += theSumResp;
1001     theResponseOfChargedWithoutEff += theSumResp;
1002     theSumPtWithEff += theSumPt;
1003     theSumPtWithoutEff += theSumPt;
1004     theSumEnergyWithEff += theSumEnergy;
1005     theSumEnergyWithoutEff += theSumEnergy;
1006   }
1007   return correction;
1008 }
1009 
1010 // -----------------------------------------------------------------------------
1011 //
1012 JetPlusTrackCorrector::P4 JetPlusTrackCorrector::pionEfficiency(const P4& jet,
1013                                                                 const Efficiency& eff,
1014                                                                 bool in_cone_at_calo_face) {
1015   // Total correction to be applied
1016   P4 correction;
1017 
1018   double theSumResp = 0;
1019   double theSumPt = 0;
1020   double theSumEnergy = 0;
1021 
1022   // Iterate through eta/pt bins
1023   for (uint32_t ieta = 0; ieta < response_.nEtaBins() - 1; ++ieta) {
1024     for (uint32_t ipt = 0; ipt < response_.nPtBins() - 1; ++ipt) {
1025       // Check tracks are found in this eta/pt bin
1026       if (!eff.nTrks(ieta, ipt)) {
1027         continue;
1028       }
1029 
1030       for (uint16_t ii = 0; ii < 2; ++ii) {
1031         // Check which correction should be applied
1032         double corr = 0.;
1033         if (ii == 0) {
1034           corr = eff.outOfConeCorr(ieta, ipt);
1035         } else if (ii == 1 && in_cone_at_calo_face) {
1036           corr = eff.inConeCorr(ieta, ipt);
1037         } else {
1038           continue;
1039         }
1040 
1041         // Calculate correction to be applied
1042         P4 corr_p4;
1043         if (vectorial_ && vecResponse_) {
1044           double corr_eta = response_.binCenterEta(ieta);
1045           double corr_phi = jet.phi();  //@@ jet phi!
1046           double corr_pt = response_.binCenterPt(ipt);
1047           corr_p4 =
1048               PtEtaPhiM(corr_pt, corr_eta, corr_phi, pionMass_);  //@@ E^2 = p^2 + m_pion^2, |p| calc'ed from pt bin
1049           corr_p4 *= (corr_p4.energy() > 0. ? corr / corr_p4.energy() : 1.);  //@@ p4 scaled up by mean energy for bin
1050         } else {
1051           corr_p4 = (jet.energy() > 0. ? corr / jet.energy() : 1.) * jet;
1052         }
1053 
1054         // Apply correction
1055         if (ii == 0) {
1056           correction += corr_p4;
1057           theSumPt += corr_p4.pt();
1058           theSumEnergy += corr_p4.energy();
1059         }  //@@ Add out-of-cone
1060         else if (ii == 1) {
1061           correction -= corr_p4;
1062           theSumResp += corr_p4.energy();
1063         }  //@@ Subtract in-cone
1064       }
1065     }
1066   }
1067 
1068   theResponseOfChargedWithEff += theSumResp;
1069   theSumPtWithEff += theSumPt;
1070   theSumEnergyWithEff += theSumEnergy;
1071   return correction;
1072 }
1073 
1074 // -----------------------------------------------------------------------------
1075 //
1076 bool JetPlusTrackCorrector::matchMuons(TrackRefs::const_iterator& itrk, const edm::Handle<RecoMuons>& muons) const {
1077   if (muons->empty()) {
1078     return false;
1079   }
1080 
1081   RecoMuons::const_iterator imuon = muons->begin();
1082   RecoMuons::const_iterator jmuon = muons->end();
1083   for (; imuon != jmuon; ++imuon) {
1084     if (imuon->innerTrack().isNull() || !muon::isGoodMuon(*imuon, muon::TMLastStationTight) ||
1085         imuon->innerTrack()->pt() < 3.0) {
1086       continue;
1087     }
1088 
1089     if (itrk->id() != imuon->innerTrack().id()) {
1090       edm::LogError("JetPlusTrackCorrector") << "[JetPlusTrackCorrector::" << __func__ << "]"
1091                                              << "Product id of the tracks associated to the jet " << itrk->id()
1092                                              << " is different from the product id of the inner track used for muons "
1093                                              << imuon->innerTrack().id() << "!" << std::endl
1094                                              << "Cannot compare tracks from different collection. Configuration Error!";
1095       return false;
1096     }
1097 
1098     if (*itrk == imuon->innerTrack())
1099       return true;
1100   }
1101 
1102   return false;
1103 }
1104 
1105 bool JetPlusTrackCorrector::matchMuons(TrackRefs::const_iterator& itrk,
1106                                        const edm::Handle<pat::MuonCollection>& muons) const {
1107   if (muons->empty()) {
1108     return false;
1109   }
1110 
1111   for (auto const& muon : *muons) {
1112     if (muon.innerTrack().isNull() == 1)
1113       continue;
1114     if (std::abs((**itrk).pt() - muon.innerTrack()->pt()) < muonPtmatch_ &&
1115         std::abs((**itrk).eta() - muon.innerTrack()->eta()) < muonEtamatch_ &&
1116         std::abs(reco::deltaPhi((**itrk).phi(), muon.innerTrack()->phi())) < muonPhimatch_) {
1117       return true;
1118     }
1119   }
1120 
1121   return false;
1122 }
1123 
1124 // -----------------------------------------------------------------------------
1125 //
1126 bool JetPlusTrackCorrector::matchElectrons(TrackRefs::const_iterator& itrk,
1127                                            const edm::Handle<RecoElectrons>& elecs,
1128                                            const edm::Handle<RecoElectronIds>& elec_ids) const {
1129   if (elecs->empty()) {
1130     return false;
1131   }
1132 
1133   double deltaRMIN = 999.;
1134 
1135   uint32_t electron_index = 0;
1136   for (auto const& ielec : *elecs) {
1137     edm::Ref<RecoElectrons> electron_ref(elecs, electron_index);
1138     electron_index++;
1139 
1140     if ((*elec_ids)[electron_ref] < 1.e-6) {
1141       continue;
1142     }  //@@ Check for null value
1143 
1144     // DR matching b/w electron and track
1145     auto dR2 = deltaR2(ielec, **itrk);
1146     if (dR2 < deltaRMIN) {
1147       deltaRMIN = dR2;
1148     }
1149   }
1150   return deltaRMIN < electronDRmatch_ * electronDRmatch_;
1151 }
1152 
1153 bool JetPlusTrackCorrector::matchElectrons(TrackRefs::const_iterator& itrk,
1154                                            const edm::Handle<pat::ElectronCollection>& elecs) const {
1155   if (elecs->empty()) {
1156     return false;
1157   }
1158 
1159   double deltaRMIN = 999.;
1160   for (auto const& ielec : *elecs) {
1161     auto dR2 = deltaR2(ielec, **itrk);
1162     if (dR2 < deltaRMIN) {
1163       deltaRMIN = dR2;
1164     }
1165   }
1166   return deltaRMIN < electronDRmatch_ * electronDRmatch_;
1167 }
1168 
1169 // -----------------------------------------------------------------------------
1170 //
1171 void JetPlusTrackCorrector::rebuildJta(const reco::Jet& fJet,
1172                                        const JetTracksAssociations& jtV0,
1173                                        TrackRefs& tracksthis,
1174                                        TrackRefs& Excl) const {
1175   tracksthis = reco::JetTracksAssociation::getValue(jtV0, fJet);
1176 
1177   if (jetSplitMerge_ < 0)
1178     return;
1179 
1180   typedef std::vector<reco::JetBaseRef>::iterator JetBaseRefIterator;
1181   std::vector<reco::JetBaseRef> theJets = reco::JetTracksAssociation::allJets(jtV0);
1182 
1183   TrackRefs tracks = tracksthis;
1184   tracksthis.clear();
1185 
1186   int tr = 0;
1187 
1188   double jetEta = fJet.eta();
1189   double jetPhi = fJet.phi();
1190   double jetEtIn = 1.0 / fJet.et();
1191 
1192   for (TrackRefs::iterator it = tracks.begin(); it != tracks.end(); it++) {
1193     double trkEta = (**it).eta();
1194     double trkPhi = (**it).phi();
1195     double dR2this = deltaR2(jetEta, jetPhi, trkEta, trkPhi);
1196     //       double dfi = fabs(fJet.phi()-(**it).phi());
1197     //       if(dfi>4.*atan(1.))dfi = 8.*atan(1.)-dfi;
1198     //       double deta = fJet.eta() - (**it).eta();
1199     //       double dR2check = sqrt(dfi*dfi+deta*deta);
1200 
1201     double scalethis = dR2this;
1202     if (jetSplitMerge_ == 0)
1203       scalethis = 1. * jetEtIn;
1204     if (jetSplitMerge_ == 2)
1205       scalethis = dR2this * jetEtIn;
1206     tr++;
1207     int flag = 1;
1208     for (JetBaseRefIterator ii = theJets.begin(); ii != theJets.end(); ii++) {
1209       if (&(**ii) == &fJet) {
1210         continue;
1211       }
1212       double dR2 = deltaR2((*ii)->eta(), (*ii)->phi(), trkEta, trkPhi);
1213       double scale = dR2;
1214       if (jetSplitMerge_ == 0)
1215         scale = 1. / (**ii).et();
1216       if (jetSplitMerge_ == 2)
1217         scale = dR2 / (**ii).et();
1218       if (scale < scalethis)
1219         flag = 0;
1220 
1221       if (flag == 0) {
1222         break;
1223       }
1224     }
1225 
1226     if (flag == 1) {
1227       tracksthis.push_back(*it);
1228     } else {
1229       Excl.push_back(*it);
1230     }
1231   }
1232 
1233   return;
1234 }
1235 
1236 // -----------------------------------------------------------------------------
1237 //
1238 void JetPlusTrackCorrector::excludeJta(const reco::Jet& fJet,
1239                                        const JetTracksAssociations& jtV0,
1240                                        TrackRefs& tracksthis,
1241                                        const TrackRefs& Excl) const {
1242   tracksthis = reco::JetTracksAssociation::getValue(jtV0, fJet);
1243   if (Excl.empty())
1244     return;
1245   if (jetSplitMerge_ < 0)
1246     return;
1247 
1248   TrackRefs tracks = tracksthis;
1249   tracksthis.clear();
1250 
1251   for (TrackRefs::iterator it = tracks.begin(); it != tracks.end(); it++) {
1252     TrackRefs::iterator itold = find(Excl.begin(), Excl.end(), (*it));
1253     if (itold == Excl.end()) {
1254       tracksthis.push_back(*it);
1255     }
1256   }
1257 
1258   return;
1259 }
1260 
1261 //================================================================================================
1262 
1263 double JetPlusTrackCorrector::correctAA(const reco::Jet& fJet,
1264                                         const reco::TrackRefVector& trBgOutOfVertex,
1265                                         double& mConeSize,
1266                                         const reco::TrackRefVector& pioninin,
1267                                         const reco::TrackRefVector& pioninout,
1268                                         double ja,
1269                                         const reco::TrackRefVector& trBgOutOfCalo) const {
1270   double mScale = 1.;
1271   double NewResponse = fJet.energy();
1272 
1273   if (trBgOutOfVertex.empty())
1274     return mScale;
1275   double EnergyOfBackgroundCharged = 0.;
1276   double ResponseOfBackgroundCharged = 0.;
1277 
1278   //
1279   // calculate the mean response
1280   //
1281 
1282   //================= EnergyOfBackgroundCharged ==================>
1283   for (reco::TrackRefVector::iterator iBgtV = trBgOutOfVertex.begin(); iBgtV != trBgOutOfVertex.end(); iBgtV++) {
1284     double eta = fabs((**iBgtV).eta());
1285     double pt = fabs((**iBgtV).pt());
1286     uint32_t ieta = response_.etaBin(eta);
1287     uint32_t ipt = response_.ptBin(pt);
1288 
1289     if (fabs(fJet.eta() - (**iBgtV).eta()) > mConeSize)
1290       continue;
1291 
1292     double echarBg = sqrt((**iBgtV).px() * (**iBgtV).px() + (**iBgtV).py() * (**iBgtV).py() +
1293                           (**iBgtV).pz() * (**iBgtV).pz() + 0.14 * 0.14);
1294 
1295     EnergyOfBackgroundCharged += echarBg / efficiency_.value(ieta, ipt);
1296 
1297   }  // Energy BG tracks
1298 
1299   //============= ResponseOfBackgroundCharged =======================>
1300 
1301   for (reco::TrackRefVector::iterator iBgtC = trBgOutOfCalo.begin(); iBgtC != trBgOutOfCalo.end(); iBgtC++) {
1302     double eta = fabs((**iBgtC).eta());
1303     double pt = fabs((**iBgtC).pt());
1304     uint32_t ieta = response_.etaBin(eta);
1305     uint32_t ipt = response_.ptBin(pt);
1306 
1307     if (fabs(fJet.eta() - (**iBgtC).eta()) > mConeSize)
1308       continue;
1309 
1310     // Check bins (not for mips)
1311     if (ieta >= response_.nEtaBins()) {
1312       continue;
1313     }
1314     if (ipt >= response_.nPtBins()) {
1315       ipt = response_.nPtBins() - 1;
1316     }
1317 
1318     double echarBg = sqrt((**iBgtC).px() * (**iBgtC).px() + (**iBgtC).py() * (**iBgtC).py() +
1319                           (**iBgtC).pz() * (**iBgtC).pz() + 0.14 * 0.14);
1320 
1321     ResponseOfBackgroundCharged += echarBg * response_.value(ieta, ipt) / efficiency_.value(ieta, ipt);
1322 
1323   }  // Response of BG tracks
1324 
1325   //=================================================================>
1326 
1327   //=================================================================>
1328   // Look for in-out tracks
1329 
1330   double en = 0.;
1331   double px = 0.;
1332   double py = 0.;
1333   double pz = 0.;
1334 
1335   for (reco::TrackRefVector::const_iterator it = pioninout.begin(); it != pioninout.end(); it++) {
1336     px += (*it)->px();
1337     py += (*it)->py();
1338     pz += (*it)->pz();
1339     en += sqrt((*it)->p() * (*it)->p() + 0.14 * 0.14);
1340   }
1341 
1342   // Look for in-in tracks
1343 
1344   double en_in = 0.;
1345   double px_in = 0.;
1346   double py_in = 0.;
1347   double pz_in = 0.;
1348 
1349   for (reco::TrackRefVector::const_iterator it = pioninin.begin(); it != pioninin.end(); it++) {
1350     px_in += (*it)->px();
1351     py_in += (*it)->py();
1352     pz_in += (*it)->pz();
1353     en_in += sqrt((*it)->p() * (*it)->p() + 0.14 * 0.14);
1354   }
1355 
1356   //===================================================================>
1357 
1358   //=>
1359   double SquareEtaRingWithoutJets = ja;
1360 
1361   EnergyOfBackgroundCharged = EnergyOfBackgroundCharged / SquareEtaRingWithoutJets;
1362   ResponseOfBackgroundCharged = ResponseOfBackgroundCharged / SquareEtaRingWithoutJets;
1363 
1364   EnergyOfBackgroundCharged = M_PI * mConeSize * mConeSize * EnergyOfBackgroundCharged;
1365   ResponseOfBackgroundCharged = M_PI * mConeSize * mConeSize * ResponseOfBackgroundCharged;
1366 
1367   NewResponse = NewResponse - EnergyOfBackgroundCharged + ResponseOfBackgroundCharged;
1368 
1369   mScale = NewResponse / fJet.energy();
1370   if (mScale < 0.)
1371     mScale = 0.;
1372   return mScale;
1373 }
1374 // -----------------------------------------------------------------------------
1375 //
1376 Map::Map(std::string input, bool verbose) : eta_(), pt_(), data_() {
1377   // Some init
1378   clear();
1379   std::vector<Element> data;
1380 
1381   // Parse file
1382   std::string file = edm::FileInPath(input).fullPath();
1383   std::ifstream in(file.c_str());
1384   string line;
1385   uint32_t ieta_old = 0;
1386   while (std::getline(in, line)) {
1387     if (line.empty() || line[0] == '#') {
1388       continue;
1389     }
1390     std::istringstream ss(line);
1391     Element temp;
1392     ss >> temp.ieta_ >> temp.ipt_ >> temp.eta_ >> temp.pt_ >> temp.val_;
1393     data.push_back(temp);
1394     if (!ieta_old || temp.ieta_ != ieta_old) {
1395       if (eta_.size() < temp.ieta_ + 1) {
1396         eta_.resize(temp.ieta_ + 1, 0.);
1397       }
1398       eta_[temp.ieta_] = temp.eta_;
1399       ieta_old = temp.ieta_;
1400     }
1401     if (pt_.size() < temp.ipt_ + 1) {
1402       pt_.resize(temp.ipt_ + 1, 0.);
1403     }
1404     pt_[temp.ipt_] = temp.pt_;
1405   }
1406 
1407   // Populate container
1408   data_.resize(eta_.size(), VDouble(pt_.size(), 0.));
1409   std::vector<Element>::const_iterator idata = data.begin();
1410   std::vector<Element>::const_iterator jdata = data.end();
1411   for (; idata != jdata; ++idata) {
1412     data_[idata->ieta_][idata->ipt_] = idata->val_;
1413   }
1414 
1415   // Check
1416   if (data_.empty() || data_[0].empty()) {
1417     std::stringstream ss;
1418     ss << "[jpt::Map::" << __func__ << "]"
1419        << " Problem parsing map in location \"" << file << "\"! ";
1420     edm::LogError("JetPlusTrackCorrector") << ss.str();
1421   }
1422 
1423   // Check
1424   if (eta_.size() != data_.size() || pt_.size() != (data_.empty() ? 0 : data_[0].size())) {
1425     std::stringstream ss;
1426     ss << "[jpt::Map::" << __func__ << "]"
1427        << " Discrepancy b/w number of bins!";
1428     edm::LogError("JetPlusTrackCorrector") << ss.str();
1429   }
1430 
1431   // Debug
1432   if (verbose && edm::isDebugEnabled()) {
1433     std::stringstream ss;
1434     ss << "[jpt::Map::" << __func__ << "]"
1435        << " Parsed contents of map at location:" << std::endl
1436        << "\"" << file << "\"" << std::endl;
1437     print(ss);
1438     LogTrace("JetPlusTrackCorrector") << ss.str();
1439   }
1440 }
1441 
1442 // -----------------------------------------------------------------------------
1443 //
1444 Map::Map() : eta_(), pt_(), data_() { clear(); }
1445 
1446 // -----------------------------------------------------------------------------
1447 //
1448 Map::~Map() { clear(); }
1449 
1450 // -----------------------------------------------------------------------------
1451 //
1452 void Map::clear() {
1453   eta_.clear();
1454   pt_.clear();
1455   data_.clear();
1456 }
1457 // -----------------------------------------------------------------------------
1458 //
1459 double Map::eta(uint32_t eta_bin) const {
1460   if (!eta_.empty() && eta_bin < eta_.size()) {
1461     return eta_[eta_bin];
1462   } else {
1463     //    edm::LogWarning("JetPlusTrackCorrector")
1464     //      << "[jpt::Map::" << __func__ << "]"
1465     //      << " Trying to access element " << eta_bin
1466     //      << " of a vector with size " << eta_.size()
1467     //      << "!";
1468     return eta_[eta_.size() - 1];
1469   }
1470 }
1471 
1472 // -----------------------------------------------------------------------------
1473 //
1474 double Map::pt(uint32_t pt_bin) const {
1475   if (!pt_.empty() && pt_bin < pt_.size()) {
1476     return pt_[pt_bin];
1477   } else {
1478     //    edm::LogWarning("JetPlusTrackCorrector")
1479     //      << "[jpt::Map::" << __func__ << "]"
1480     //      << " Trying to access element " << pt_bin
1481     //      << " of a vector with size " << pt_.size()
1482     //      << "!";
1483     return pt_[pt_.size() - 1];
1484   }
1485 }
1486 
1487 // -----------------------------------------------------------------------------
1488 //
1489 double Map::binCenterEta(uint32_t eta_bin) const {
1490   if (!eta_.empty() && eta_bin + 1 < eta_.size()) {
1491     return eta_[eta_bin] + (eta_[eta_bin + 1] - eta_[eta_bin]) / 2.;
1492   } else {
1493     //    edm::LogWarning("JetPlusTrackCorrector")
1494     //      << "[jpt::Map::" << __func__ << "]"
1495     //      << " Trying to access element " << eta_bin+1
1496     //      << " of a vector with size " << eta_.size()
1497     //      << "!";
1498     return eta_[eta_.size() - 1];
1499   }
1500 }
1501 
1502 // -----------------------------------------------------------------------------
1503 //
1504 double Map::binCenterPt(uint32_t pt_bin) const {
1505   if (!pt_.empty() && pt_bin + 1 < pt_.size()) {
1506     return pt_[pt_bin] + (pt_[pt_bin + 1] - pt_[pt_bin]) / 2.;
1507   } else {
1508     //    edm::LogWarning("JetPlusTrackCorrector")
1509     //      << "[jpt::Map::" << __func__ << "]"
1510     //      << " Trying to access element " << pt_bin+1
1511     //      << " of a vector with size " << pt_.size()
1512     //      << "!";
1513     return pt_[pt_.size() - 1];
1514   }
1515 }
1516 
1517 // -----------------------------------------------------------------------------
1518 //
1519 uint32_t Map::etaBin(double val) const {
1520   val = fabs(val);
1521   for (uint32_t ieta = 0; ieta < nEtaBins() - 1; ++ieta) {  //@@ "-1" is bug?
1522     if (val > eta(ieta) && (ieta + 1 == nEtaBins() || val < eta(ieta + 1))) {
1523       return ieta;
1524     }
1525   }
1526   return nEtaBins();
1527 }
1528 
1529 // -----------------------------------------------------------------------------
1530 //
1531 uint32_t Map::ptBin(double val) const {
1532   val = fabs(val);
1533   for (uint32_t ipt = 0; ipt < nPtBins() - 1; ++ipt) {  //@@ "-1" is bug?
1534     if (val > pt(ipt) && ((ipt + 1) == nPtBins() || val < pt(ipt + 1))) {
1535       return ipt;
1536     }
1537   }
1538   return nPtBins();
1539 }
1540 
1541 // -----------------------------------------------------------------------------
1542 //
1543 double Map::value(uint32_t eta_bin, uint32_t pt_bin) const {
1544   if (eta_bin < data_.size() && pt_bin < (data_.empty() ? 0 : data_[0].size())) {
1545     return data_[eta_bin][pt_bin];
1546   } else {
1547     //    edm::LogWarning("JetPlusTrackCorrector")
1548     //      << "[jpt::Map::" << __func__ << "]"
1549     //      << " Trying to access element (" << eta_bin << "," << pt_bin << ")"
1550     //      << " of a vector with size (" << data_.size() << "," << ( data_.empty() ? 0 : data_[0].size() ) << ")"
1551     //      << "!";
1552     return 1.;
1553   }
1554 }
1555 
1556 // -----------------------------------------------------------------------------
1557 //
1558 void Map::print(std::stringstream& ss) const {
1559   ss << " Number of bins in eta : " << data_.size() << std::endl
1560      << " Number of bins in pt  : " << (data_.empty() ? 0 : data_[0].size()) << std::endl;
1561   VVDouble::const_iterator ieta = data_.begin();
1562   VVDouble::const_iterator jeta = data_.end();
1563   for (; ieta != jeta; ++ieta) {
1564     VDouble::const_iterator ipt = ieta->begin();
1565     VDouble::const_iterator jpt = ieta->end();
1566     for (; ipt != jpt; ++ipt) {
1567       uint32_t eta_bin = static_cast<uint32_t>(ieta - data_.begin());
1568       uint32_t pt_bin = static_cast<uint32_t>(ipt - ieta->begin());
1569       ss << " EtaBinNumber: " << eta_bin << " PtBinNumber: " << pt_bin << " EtaValue: " << eta_[eta_bin]
1570          << " PtValue: " << pt_[pt_bin] << " Value: " << data_[eta_bin][pt_bin] << std::endl;
1571     }
1572   }
1573 }
1574 
1575 // -----------------------------------------------------------------------------
1576 //
1577 MatchedTracks::MatchedTracks() : inVertexInCalo_(), outOfVertexInCalo_(), inVertexOutOfCalo_() { clear(); }
1578 
1579 // -----------------------------------------------------------------------------
1580 //
1581 MatchedTracks::~MatchedTracks() { clear(); }
1582 
1583 // -----------------------------------------------------------------------------
1584 //
1585 void MatchedTracks::clear() {
1586   inVertexInCalo_.clear();
1587   outOfVertexInCalo_.clear();
1588   inVertexOutOfCalo_.clear();
1589 }
1590 
1591 // -----------------------------------------------------------------------------
1592 //
1593 JetTracks::JetTracks() : vertex_(), caloFace_() { clear(); }
1594 
1595 // -----------------------------------------------------------------------------
1596 //
1597 JetTracks::~JetTracks() { clear(); }
1598 
1599 // -----------------------------------------------------------------------------
1600 //
1601 void JetTracks::clear() {
1602   vertex_.clear();
1603   caloFace_.clear();
1604 }
1605 
1606 // -----------------------------------------------------------------------------
1607 //
1608 Efficiency::Efficiency(const jpt::Map& response, const jpt::Map& efficiency, const jpt::Map& leakage)
1609     : response_(response), efficiency_(efficiency), leakage_(leakage) {
1610   reset();
1611 }
1612 
1613 // -----------------------------------------------------------------------------
1614 //
1615 double Efficiency::inConeCorr(uint32_t eta_bin, uint32_t pt_bin) const {
1616   if (check(eta_bin, pt_bin, __func__)) {
1617     return (outOfConeCorr(eta_bin, pt_bin) * leakage_.value(eta_bin, pt_bin) * response_.value(eta_bin, pt_bin));
1618   } else {
1619     return 0.;
1620   }
1621 }
1622 
1623 // -----------------------------------------------------------------------------
1624 //
1625 double Efficiency::outOfConeCorr(uint32_t eta_bin, uint32_t pt_bin) const {
1626   if (check(eta_bin, pt_bin, __func__)) {
1627     uint16_t ntrks = nTrks(eta_bin, pt_bin);
1628     double mean = meanE(eta_bin, pt_bin);
1629     double eff = (1. - efficiency_.value(eta_bin, pt_bin)) / efficiency_.value(eta_bin, pt_bin);
1630     if (!ntrks) {
1631       return 0.;
1632     }
1633     return (ntrks * eff * mean);
1634   } else {
1635     return 0.;
1636   }
1637 }
1638 
1639 // -----------------------------------------------------------------------------
1640 //
1641 uint16_t Efficiency::nTrks(uint32_t eta_bin, uint32_t pt_bin) const {
1642   if (check(eta_bin, pt_bin, __func__)) {
1643     return data_[eta_bin][pt_bin].first;
1644   } else {
1645     return 0;
1646   }
1647 }
1648 
1649 // -----------------------------------------------------------------------------
1650 //
1651 double Efficiency::sumE(uint32_t eta_bin, uint32_t pt_bin) const {
1652   if (check(eta_bin, pt_bin, __func__)) {
1653     return data_[eta_bin][pt_bin].second;
1654   } else {
1655     return 0.;
1656   }
1657 }
1658 
1659 // -----------------------------------------------------------------------------
1660 //
1661 double Efficiency::meanE(uint32_t eta_bin, uint32_t pt_bin) const {
1662   if (check(eta_bin, pt_bin, __func__)) {
1663     Pair tmp = data_[eta_bin][pt_bin];
1664     if (tmp.first) {
1665       return tmp.second / tmp.first;
1666     } else {
1667       return 0.;
1668     }
1669   } else {
1670     return 0.;
1671   }
1672 }
1673 
1674 // -----------------------------------------------------------------------------
1675 //
1676 void Efficiency::addE(uint32_t eta_bin, uint32_t pt_bin, double energy) {
1677   if (check(eta_bin, pt_bin, __func__)) {
1678     data_[eta_bin][pt_bin].first++;
1679     data_[eta_bin][pt_bin].second += energy;
1680   }
1681 }
1682 
1683 // -----------------------------------------------------------------------------
1684 //
1685 bool Efficiency::check(uint32_t eta_bin, uint32_t pt_bin, std::string method) const {
1686   if (eta_bin < data_.size() && pt_bin < (data_.empty() ? 0 : data_[0].size())) {
1687     return true;
1688   } else {
1689     //    edm::LogWarning("JetPlusTrackCorrector")
1690     //      << "[jpt::Efficiency::" << method << "]"
1691     //      << " Trying to access element (" << eta_bin << "," << pt_bin << ")"
1692     //      << " of a vector with size (" << data_.size() << "," << ( data_.empty() ? 0 : data_[0].size() ) << ")"
1693     //      << "!";
1694     return false;
1695   }
1696 }
1697 
1698 // -----------------------------------------------------------------------------
1699 //
1700 void Efficiency::reset() {
1701   data_.clear();
1702   data_.resize(response_.nEtaBins(), VPair(response_.nPtBins(), Pair(0, 0.)));
1703 }
1704 
1705 // -----------------------------------------------------------------------------
1706 //
1707 void Efficiency::print() const {
1708   if (edm::isDebugEnabled()) {
1709     std::stringstream ss;
1710     ss << "[jpt::Efficiency::" << __func__ << "]"
1711        << " Contents of maps:" << std::endl;
1712     ss << "Response map: " << std::endl;
1713     response_.print(ss);
1714     ss << "Efficiency map: " << std::endl;
1715     efficiency_.print(ss);
1716     ss << "Leakage map: " << std::endl;
1717     leakage_.print(ss);
1718     LogTrace("JetPlusTrackCorrector") << ss.str();
1719   }
1720 }