Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 09:59:48

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   double jetEta = fJet.eta();
1187   double jetPhi = fJet.phi();
1188   double jetEtIn = 1.0 / fJet.et();
1189 
1190   for (TrackRefs::iterator it = tracks.begin(); it != tracks.end(); it++) {
1191     double trkEta = (**it).eta();
1192     double trkPhi = (**it).phi();
1193     double dR2this = deltaR2(jetEta, jetPhi, trkEta, trkPhi);
1194     //       double dfi = fabs(fJet.phi()-(**it).phi());
1195     //       if(dfi>4.*atan(1.))dfi = 8.*atan(1.)-dfi;
1196     //       double deta = fJet.eta() - (**it).eta();
1197     //       double dR2check = sqrt(dfi*dfi+deta*deta);
1198 
1199     double scalethis = dR2this;
1200     if (jetSplitMerge_ == 0)
1201       scalethis = 1. * jetEtIn;
1202     if (jetSplitMerge_ == 2)
1203       scalethis = dR2this * jetEtIn;
1204     int flag = 1;
1205     for (JetBaseRefIterator ii = theJets.begin(); ii != theJets.end(); ii++) {
1206       if (&(**ii) == &fJet) {
1207         continue;
1208       }
1209       double dR2 = deltaR2((*ii)->eta(), (*ii)->phi(), trkEta, trkPhi);
1210       double scale = dR2;
1211       if (jetSplitMerge_ == 0)
1212         scale = 1. / (**ii).et();
1213       if (jetSplitMerge_ == 2)
1214         scale = dR2 / (**ii).et();
1215       if (scale < scalethis)
1216         flag = 0;
1217 
1218       if (flag == 0) {
1219         break;
1220       }
1221     }
1222 
1223     if (flag == 1) {
1224       tracksthis.push_back(*it);
1225     } else {
1226       Excl.push_back(*it);
1227     }
1228   }
1229 
1230   return;
1231 }
1232 
1233 // -----------------------------------------------------------------------------
1234 //
1235 void JetPlusTrackCorrector::excludeJta(const reco::Jet& fJet,
1236                                        const JetTracksAssociations& jtV0,
1237                                        TrackRefs& tracksthis,
1238                                        const TrackRefs& Excl) const {
1239   tracksthis = reco::JetTracksAssociation::getValue(jtV0, fJet);
1240   if (Excl.empty())
1241     return;
1242   if (jetSplitMerge_ < 0)
1243     return;
1244 
1245   TrackRefs tracks = tracksthis;
1246   tracksthis.clear();
1247 
1248   for (TrackRefs::iterator it = tracks.begin(); it != tracks.end(); it++) {
1249     TrackRefs::iterator itold = find(Excl.begin(), Excl.end(), (*it));
1250     if (itold == Excl.end()) {
1251       tracksthis.push_back(*it);
1252     }
1253   }
1254 
1255   return;
1256 }
1257 
1258 //================================================================================================
1259 
1260 double JetPlusTrackCorrector::correctAA(const reco::Jet& fJet,
1261                                         const reco::TrackRefVector& trBgOutOfVertex,
1262                                         double& mConeSize,
1263                                         const reco::TrackRefVector& pioninin,
1264                                         const reco::TrackRefVector& pioninout,
1265                                         double ja,
1266                                         const reco::TrackRefVector& trBgOutOfCalo) const {
1267   double mScale = 1.;
1268   double NewResponse = fJet.energy();
1269 
1270   if (trBgOutOfVertex.empty())
1271     return mScale;
1272   double EnergyOfBackgroundCharged = 0.;
1273   double ResponseOfBackgroundCharged = 0.;
1274 
1275   //
1276   // calculate the mean response
1277   //
1278 
1279   //================= EnergyOfBackgroundCharged ==================>
1280   for (reco::TrackRefVector::iterator iBgtV = trBgOutOfVertex.begin(); iBgtV != trBgOutOfVertex.end(); iBgtV++) {
1281     double eta = fabs((**iBgtV).eta());
1282     double pt = fabs((**iBgtV).pt());
1283     uint32_t ieta = response_.etaBin(eta);
1284     uint32_t ipt = response_.ptBin(pt);
1285 
1286     if (fabs(fJet.eta() - (**iBgtV).eta()) > mConeSize)
1287       continue;
1288 
1289     double echarBg = sqrt((**iBgtV).px() * (**iBgtV).px() + (**iBgtV).py() * (**iBgtV).py() +
1290                           (**iBgtV).pz() * (**iBgtV).pz() + 0.14 * 0.14);
1291 
1292     EnergyOfBackgroundCharged += echarBg / efficiency_.value(ieta, ipt);
1293 
1294   }  // Energy BG tracks
1295 
1296   //============= ResponseOfBackgroundCharged =======================>
1297 
1298   for (reco::TrackRefVector::iterator iBgtC = trBgOutOfCalo.begin(); iBgtC != trBgOutOfCalo.end(); iBgtC++) {
1299     double eta = fabs((**iBgtC).eta());
1300     double pt = fabs((**iBgtC).pt());
1301     uint32_t ieta = response_.etaBin(eta);
1302     uint32_t ipt = response_.ptBin(pt);
1303 
1304     if (fabs(fJet.eta() - (**iBgtC).eta()) > mConeSize)
1305       continue;
1306 
1307     // Check bins (not for mips)
1308     if (ieta >= response_.nEtaBins()) {
1309       continue;
1310     }
1311     if (ipt >= response_.nPtBins()) {
1312       ipt = response_.nPtBins() - 1;
1313     }
1314 
1315     double echarBg = sqrt((**iBgtC).px() * (**iBgtC).px() + (**iBgtC).py() * (**iBgtC).py() +
1316                           (**iBgtC).pz() * (**iBgtC).pz() + 0.14 * 0.14);
1317 
1318     ResponseOfBackgroundCharged += echarBg * response_.value(ieta, ipt) / efficiency_.value(ieta, ipt);
1319 
1320   }  // Response of BG tracks
1321 
1322   //=================================================================>
1323 
1324   /*
1325   //=================================================================>
1326   // Look for in-out tracks
1327 
1328   double en = 0.;
1329   double px = 0.;
1330   double py = 0.;
1331   double pz = 0.;
1332 
1333   for (reco::TrackRefVector::const_iterator it = pioninout.begin(); it != pioninout.end(); it++) {
1334     px += (*it)->px();
1335     py += (*it)->py();
1336     pz += (*it)->pz();
1337     en += sqrt((*it)->p() * (*it)->p() + 0.14 * 0.14);
1338   }
1339 
1340   // Look for in-in tracks
1341 
1342   double en_in = 0.;
1343   double px_in = 0.;
1344   double py_in = 0.;
1345   double pz_in = 0.;
1346 
1347   for (reco::TrackRefVector::const_iterator it = pioninin.begin(); it != pioninin.end(); it++) {
1348     px_in += (*it)->px();
1349     py_in += (*it)->py();
1350     pz_in += (*it)->pz();
1351     en_in += sqrt((*it)->p() * (*it)->p() + 0.14 * 0.14);
1352   }
1353 
1354   //===================================================================>
1355 
1356   //=>
1357 */
1358   double SquareEtaRingWithoutJets = ja;
1359 
1360   EnergyOfBackgroundCharged = EnergyOfBackgroundCharged / SquareEtaRingWithoutJets;
1361   ResponseOfBackgroundCharged = ResponseOfBackgroundCharged / SquareEtaRingWithoutJets;
1362 
1363   EnergyOfBackgroundCharged = M_PI * mConeSize * mConeSize * EnergyOfBackgroundCharged;
1364   ResponseOfBackgroundCharged = M_PI * mConeSize * mConeSize * ResponseOfBackgroundCharged;
1365 
1366   NewResponse = NewResponse - EnergyOfBackgroundCharged + ResponseOfBackgroundCharged;
1367 
1368   mScale = NewResponse / fJet.energy();
1369   if (mScale < 0.)
1370     mScale = 0.;
1371   return mScale;
1372 }
1373 // -----------------------------------------------------------------------------
1374 //
1375 Map::Map(std::string input, bool verbose) : eta_(), pt_(), data_() {
1376   // Some init
1377   clear();
1378   std::vector<Element> data;
1379 
1380   // Parse file
1381   std::string file = edm::FileInPath(input).fullPath();
1382   std::ifstream in(file.c_str());
1383   string line;
1384   uint32_t ieta_old = 0;
1385   while (std::getline(in, line)) {
1386     if (line.empty() || line[0] == '#') {
1387       continue;
1388     }
1389     std::istringstream ss(line);
1390     Element temp;
1391     ss >> temp.ieta_ >> temp.ipt_ >> temp.eta_ >> temp.pt_ >> temp.val_;
1392     data.push_back(temp);
1393     if (!ieta_old || temp.ieta_ != ieta_old) {
1394       if (eta_.size() < temp.ieta_ + 1) {
1395         eta_.resize(temp.ieta_ + 1, 0.);
1396       }
1397       eta_[temp.ieta_] = temp.eta_;
1398       ieta_old = temp.ieta_;
1399     }
1400     if (pt_.size() < temp.ipt_ + 1) {
1401       pt_.resize(temp.ipt_ + 1, 0.);
1402     }
1403     pt_[temp.ipt_] = temp.pt_;
1404   }
1405 
1406   // Populate container
1407   data_.resize(eta_.size(), VDouble(pt_.size(), 0.));
1408   std::vector<Element>::const_iterator idata = data.begin();
1409   std::vector<Element>::const_iterator jdata = data.end();
1410   for (; idata != jdata; ++idata) {
1411     data_[idata->ieta_][idata->ipt_] = idata->val_;
1412   }
1413 
1414   // Check
1415   if (data_.empty() || data_[0].empty()) {
1416     std::stringstream ss;
1417     ss << "[jpt::Map::" << __func__ << "]"
1418        << " Problem parsing map in location \"" << file << "\"! ";
1419     edm::LogError("JetPlusTrackCorrector") << ss.str();
1420   }
1421 
1422   // Check
1423   if (eta_.size() != data_.size() || pt_.size() != (data_.empty() ? 0 : data_[0].size())) {
1424     std::stringstream ss;
1425     ss << "[jpt::Map::" << __func__ << "]"
1426        << " Discrepancy b/w number of bins!";
1427     edm::LogError("JetPlusTrackCorrector") << ss.str();
1428   }
1429 
1430   // Debug
1431   if (verbose && edm::isDebugEnabled()) {
1432     std::stringstream ss;
1433     ss << "[jpt::Map::" << __func__ << "]"
1434        << " Parsed contents of map at location:" << std::endl
1435        << "\"" << file << "\"" << std::endl;
1436     print(ss);
1437     LogTrace("JetPlusTrackCorrector") << ss.str();
1438   }
1439 }
1440 
1441 // -----------------------------------------------------------------------------
1442 //
1443 Map::Map() : eta_(), pt_(), data_() { clear(); }
1444 
1445 // -----------------------------------------------------------------------------
1446 //
1447 Map::~Map() { clear(); }
1448 
1449 // -----------------------------------------------------------------------------
1450 //
1451 void Map::clear() {
1452   eta_.clear();
1453   pt_.clear();
1454   data_.clear();
1455 }
1456 // -----------------------------------------------------------------------------
1457 //
1458 double Map::eta(uint32_t eta_bin) const {
1459   if (!eta_.empty() && eta_bin < eta_.size()) {
1460     return eta_[eta_bin];
1461   } else {
1462     //    edm::LogWarning("JetPlusTrackCorrector")
1463     //      << "[jpt::Map::" << __func__ << "]"
1464     //      << " Trying to access element " << eta_bin
1465     //      << " of a vector with size " << eta_.size()
1466     //      << "!";
1467     return eta_[eta_.size() - 1];
1468   }
1469 }
1470 
1471 // -----------------------------------------------------------------------------
1472 //
1473 double Map::pt(uint32_t pt_bin) const {
1474   if (!pt_.empty() && pt_bin < pt_.size()) {
1475     return pt_[pt_bin];
1476   } else {
1477     //    edm::LogWarning("JetPlusTrackCorrector")
1478     //      << "[jpt::Map::" << __func__ << "]"
1479     //      << " Trying to access element " << pt_bin
1480     //      << " of a vector with size " << pt_.size()
1481     //      << "!";
1482     return pt_[pt_.size() - 1];
1483   }
1484 }
1485 
1486 // -----------------------------------------------------------------------------
1487 //
1488 double Map::binCenterEta(uint32_t eta_bin) const {
1489   if (!eta_.empty() && eta_bin + 1 < eta_.size()) {
1490     return eta_[eta_bin] + (eta_[eta_bin + 1] - eta_[eta_bin]) / 2.;
1491   } else {
1492     //    edm::LogWarning("JetPlusTrackCorrector")
1493     //      << "[jpt::Map::" << __func__ << "]"
1494     //      << " Trying to access element " << eta_bin+1
1495     //      << " of a vector with size " << eta_.size()
1496     //      << "!";
1497     return eta_[eta_.size() - 1];
1498   }
1499 }
1500 
1501 // -----------------------------------------------------------------------------
1502 //
1503 double Map::binCenterPt(uint32_t pt_bin) const {
1504   if (!pt_.empty() && pt_bin + 1 < pt_.size()) {
1505     return pt_[pt_bin] + (pt_[pt_bin + 1] - pt_[pt_bin]) / 2.;
1506   } else {
1507     //    edm::LogWarning("JetPlusTrackCorrector")
1508     //      << "[jpt::Map::" << __func__ << "]"
1509     //      << " Trying to access element " << pt_bin+1
1510     //      << " of a vector with size " << pt_.size()
1511     //      << "!";
1512     return pt_[pt_.size() - 1];
1513   }
1514 }
1515 
1516 // -----------------------------------------------------------------------------
1517 //
1518 uint32_t Map::etaBin(double val) const {
1519   val = fabs(val);
1520   for (uint32_t ieta = 0; ieta < nEtaBins() - 1; ++ieta) {  //@@ "-1" is bug?
1521     if (val > eta(ieta) && (ieta + 1 == nEtaBins() || val < eta(ieta + 1))) {
1522       return ieta;
1523     }
1524   }
1525   return nEtaBins();
1526 }
1527 
1528 // -----------------------------------------------------------------------------
1529 //
1530 uint32_t Map::ptBin(double val) const {
1531   val = fabs(val);
1532   for (uint32_t ipt = 0; ipt < nPtBins() - 1; ++ipt) {  //@@ "-1" is bug?
1533     if (val > pt(ipt) && ((ipt + 1) == nPtBins() || val < pt(ipt + 1))) {
1534       return ipt;
1535     }
1536   }
1537   return nPtBins();
1538 }
1539 
1540 // -----------------------------------------------------------------------------
1541 //
1542 double Map::value(uint32_t eta_bin, uint32_t pt_bin) const {
1543   if (eta_bin < data_.size() && pt_bin < (data_.empty() ? 0 : data_[0].size())) {
1544     return data_[eta_bin][pt_bin];
1545   } else {
1546     //    edm::LogWarning("JetPlusTrackCorrector")
1547     //      << "[jpt::Map::" << __func__ << "]"
1548     //      << " Trying to access element (" << eta_bin << "," << pt_bin << ")"
1549     //      << " of a vector with size (" << data_.size() << "," << ( data_.empty() ? 0 : data_[0].size() ) << ")"
1550     //      << "!";
1551     return 1.;
1552   }
1553 }
1554 
1555 // -----------------------------------------------------------------------------
1556 //
1557 void Map::print(std::stringstream& ss) const {
1558   ss << " Number of bins in eta : " << data_.size() << std::endl
1559      << " Number of bins in pt  : " << (data_.empty() ? 0 : data_[0].size()) << std::endl;
1560   VVDouble::const_iterator ieta = data_.begin();
1561   VVDouble::const_iterator jeta = data_.end();
1562   for (; ieta != jeta; ++ieta) {
1563     VDouble::const_iterator ipt = ieta->begin();
1564     VDouble::const_iterator jpt = ieta->end();
1565     for (; ipt != jpt; ++ipt) {
1566       uint32_t eta_bin = static_cast<uint32_t>(ieta - data_.begin());
1567       uint32_t pt_bin = static_cast<uint32_t>(ipt - ieta->begin());
1568       ss << " EtaBinNumber: " << eta_bin << " PtBinNumber: " << pt_bin << " EtaValue: " << eta_[eta_bin]
1569          << " PtValue: " << pt_[pt_bin] << " Value: " << data_[eta_bin][pt_bin] << std::endl;
1570     }
1571   }
1572 }
1573 
1574 // -----------------------------------------------------------------------------
1575 //
1576 MatchedTracks::MatchedTracks() : inVertexInCalo_(), outOfVertexInCalo_(), inVertexOutOfCalo_() { clear(); }
1577 
1578 // -----------------------------------------------------------------------------
1579 //
1580 MatchedTracks::~MatchedTracks() { clear(); }
1581 
1582 // -----------------------------------------------------------------------------
1583 //
1584 void MatchedTracks::clear() {
1585   inVertexInCalo_.clear();
1586   outOfVertexInCalo_.clear();
1587   inVertexOutOfCalo_.clear();
1588 }
1589 
1590 // -----------------------------------------------------------------------------
1591 //
1592 JetTracks::JetTracks() : vertex_(), caloFace_() { clear(); }
1593 
1594 // -----------------------------------------------------------------------------
1595 //
1596 JetTracks::~JetTracks() { clear(); }
1597 
1598 // -----------------------------------------------------------------------------
1599 //
1600 void JetTracks::clear() {
1601   vertex_.clear();
1602   caloFace_.clear();
1603 }
1604 
1605 // -----------------------------------------------------------------------------
1606 //
1607 Efficiency::Efficiency(const jpt::Map& response, const jpt::Map& efficiency, const jpt::Map& leakage)
1608     : response_(response), efficiency_(efficiency), leakage_(leakage) {
1609   reset();
1610 }
1611 
1612 // -----------------------------------------------------------------------------
1613 //
1614 double Efficiency::inConeCorr(uint32_t eta_bin, uint32_t pt_bin) const {
1615   if (check(eta_bin, pt_bin, __func__)) {
1616     return (outOfConeCorr(eta_bin, pt_bin) * leakage_.value(eta_bin, pt_bin) * response_.value(eta_bin, pt_bin));
1617   } else {
1618     return 0.;
1619   }
1620 }
1621 
1622 // -----------------------------------------------------------------------------
1623 //
1624 double Efficiency::outOfConeCorr(uint32_t eta_bin, uint32_t pt_bin) const {
1625   if (check(eta_bin, pt_bin, __func__)) {
1626     uint16_t ntrks = nTrks(eta_bin, pt_bin);
1627     double mean = meanE(eta_bin, pt_bin);
1628     double eff = (1. - efficiency_.value(eta_bin, pt_bin)) / efficiency_.value(eta_bin, pt_bin);
1629     if (!ntrks) {
1630       return 0.;
1631     }
1632     return (ntrks * eff * mean);
1633   } else {
1634     return 0.;
1635   }
1636 }
1637 
1638 // -----------------------------------------------------------------------------
1639 //
1640 uint16_t Efficiency::nTrks(uint32_t eta_bin, uint32_t pt_bin) const {
1641   if (check(eta_bin, pt_bin, __func__)) {
1642     return data_[eta_bin][pt_bin].first;
1643   } else {
1644     return 0;
1645   }
1646 }
1647 
1648 // -----------------------------------------------------------------------------
1649 //
1650 double Efficiency::sumE(uint32_t eta_bin, uint32_t pt_bin) const {
1651   if (check(eta_bin, pt_bin, __func__)) {
1652     return data_[eta_bin][pt_bin].second;
1653   } else {
1654     return 0.;
1655   }
1656 }
1657 
1658 // -----------------------------------------------------------------------------
1659 //
1660 double Efficiency::meanE(uint32_t eta_bin, uint32_t pt_bin) const {
1661   if (check(eta_bin, pt_bin, __func__)) {
1662     Pair tmp = data_[eta_bin][pt_bin];
1663     if (tmp.first) {
1664       return tmp.second / tmp.first;
1665     } else {
1666       return 0.;
1667     }
1668   } else {
1669     return 0.;
1670   }
1671 }
1672 
1673 // -----------------------------------------------------------------------------
1674 //
1675 void Efficiency::addE(uint32_t eta_bin, uint32_t pt_bin, double energy) {
1676   if (check(eta_bin, pt_bin, __func__)) {
1677     data_[eta_bin][pt_bin].first++;
1678     data_[eta_bin][pt_bin].second += energy;
1679   }
1680 }
1681 
1682 // -----------------------------------------------------------------------------
1683 //
1684 bool Efficiency::check(uint32_t eta_bin, uint32_t pt_bin, std::string method) const {
1685   if (eta_bin < data_.size() && pt_bin < (data_.empty() ? 0 : data_[0].size())) {
1686     return true;
1687   } else {
1688     //    edm::LogWarning("JetPlusTrackCorrector")
1689     //      << "[jpt::Efficiency::" << method << "]"
1690     //      << " Trying to access element (" << eta_bin << "," << pt_bin << ")"
1691     //      << " of a vector with size (" << data_.size() << "," << ( data_.empty() ? 0 : data_[0].size() ) << ")"
1692     //      << "!";
1693     return false;
1694   }
1695 }
1696 
1697 // -----------------------------------------------------------------------------
1698 //
1699 void Efficiency::reset() {
1700   data_.clear();
1701   data_.resize(response_.nEtaBins(), VPair(response_.nPtBins(), Pair(0, 0.)));
1702 }
1703 
1704 // -----------------------------------------------------------------------------
1705 //
1706 void Efficiency::print() const {
1707   if (edm::isDebugEnabled()) {
1708     std::stringstream ss;
1709     ss << "[jpt::Efficiency::" << __func__ << "]"
1710        << " Contents of maps:" << std::endl;
1711     ss << "Response map: " << std::endl;
1712     response_.print(ss);
1713     ss << "Efficiency map: " << std::endl;
1714     efficiency_.print(ss);
1715     ss << "Leakage map: " << std::endl;
1716     leakage_.print(ss);
1717     LogTrace("JetPlusTrackCorrector") << ss.str();
1718   }
1719 }