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
0169 corrected = fJet.p4();
0170
0171
0172 validMatches = matchTracks(fJetcalo, event, setup, pions, muons, elecs);
0173 if (!validMatches) {
0174 return 1.;
0175 }
0176
0177
0178 if (!canCorrect(fJet)) {
0179 return 1.;
0180 }
0181
0182
0183 if (verbose_) {
0184 edm::LogInfo("JetPlusTrackCorrector") << "[JetPlusTrackCorrector::" << __func__ << "]"
0185 << " Applying JPT corrections...";
0186 }
0187
0188
0189 if (usePions_) {
0190 corrected += pionCorrection(fJet.p4(), pions);
0191 }
0192
0193
0194 if (useMuons_) {
0195 corrected += muonCorrection(fJet.p4(), muons);
0196 }
0197
0198
0199 if (useElecs_) {
0200 corrected += elecCorrection(fJet.p4(), elecs);
0201 }
0202
0203
0204 if (vectorial_ && !vecResponse_) {
0205 if (fabs(corrected.eta()) < 2.) {
0206 corrected = jetDirFromTracks(corrected, pions, muons, elecs);
0207 }
0208 }
0209
0210
0211 double scale = checkScale(fJet.p4(), corrected);
0212
0213
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.);
0229 edm::LogVerbatim("JetPlusTrackCorrector") << ss.str();
0230 }
0231
0232
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
0275 JetTracks jet_tracks;
0276 bool ok = jetTrackAssociation(fJet, event, setup, jet_tracks);
0277
0278 if (!ok) {
0279 return false;
0280 }
0281
0282
0283 matchTracks(jet_tracks, event, pions, muons, elecs);
0284
0285
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
0304 trks.clear();
0305
0306
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
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
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
0353 if (trks.vertex_.empty()) {
0354 return false;
0355 }
0356
0357
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
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
0403 pions.clear();
0404 muons.clear();
0405 elecs.clear();
0406
0407
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
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
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
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
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
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 }
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
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
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
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
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;
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;
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
0801 P4 corr;
0802
0803
0804 bool tracks_present_inin = false;
0805
0806
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
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
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
0842 }
0843 }
0844
0845
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
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
0862 }
0863 }
0864
0865
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
0888 P4 correction;
0889
0890
0891 eff.reset();
0892
0893 double theSumResp = 0;
0894 double theSumPt = 0;
0895 double theSumEnergy = 0;
0896
0897
0898 if (!tracks.empty()) {
0899 TrackRefs::iterator itrk = tracks.begin();
0900 TrackRefs::iterator jtrk = tracks.end();
0901
0902 for (; itrk != jtrk; ++itrk) {
0903
0904 if (in_cone_at_calo_face && is_pion && (*itrk)->pt() >= 50.) {
0905 continue;
0906 }
0907
0908
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
0919 if (in_cone_at_vertex) {
0920 correction += inner;
0921 }
0922
0923
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
0930 if (is_pion && (ieta == response_.nEtaBins() || ipt == response_.nPtBins())) {
0931 continue;
0932 }
0933
0934
0935 P4 outer;
0936 if (in_cone_at_calo_face) {
0937 if (vectorial_ && vecResponse_) {
0938
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();
0945 outer_phi = (*itrk)->outerPosition().phi();
0946 }
0947 outer = PtEtaPhiM(outer_pt, outer_eta, outer_phi, mass);
0948
0949 if (!is_pion) {
0950 outer *= (outer.energy() > 0. ? mip / outer.energy() : 1.);
0951 }
0952 else {
0953 outer *= (outer.energy() > 0. ? inner.energy() / outer.energy() : 1.);
0954 }
0955 } else {
0956
0957 if (!is_pion) {
0958 outer = (jet.energy() > 0. ? mip / jet.energy() : 1.) * jet;
0959 }
0960 else {
0961 outer = inner;
0962 }
0963 }
0964 if (is_pion) {
0965 outer *= response_.value(ieta, ipt);
0966 }
0967 correction -= outer;
0968
0969
0970 theSumResp += response_.value(ieta, ipt);
0971 }
0972
0973
0974 theSumPt += inner.pt();
0975 theSumEnergy += inner.energy();
0976
0977
0978 if (is_pion) {
0979 eff.addE(ieta, ipt, inner.energy());
0980 }
0981
0982
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 }
0997 }
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
1016 P4 correction;
1017
1018 double theSumResp = 0;
1019 double theSumPt = 0;
1020 double theSumEnergy = 0;
1021
1022
1023 for (uint32_t ieta = 0; ieta < response_.nEtaBins() - 1; ++ieta) {
1024 for (uint32_t ipt = 0; ipt < response_.nPtBins() - 1; ++ipt) {
1025
1026 if (!eff.nTrks(ieta, ipt)) {
1027 continue;
1028 }
1029
1030 for (uint16_t ii = 0; ii < 2; ++ii) {
1031
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
1042 P4 corr_p4;
1043 if (vectorial_ && vecResponse_) {
1044 double corr_eta = response_.binCenterEta(ieta);
1045 double corr_phi = jet.phi();
1046 double corr_pt = response_.binCenterPt(ipt);
1047 corr_p4 =
1048 PtEtaPhiM(corr_pt, corr_eta, corr_phi, pionMass_);
1049 corr_p4 *= (corr_p4.energy() > 0. ? corr / corr_p4.energy() : 1.);
1050 } else {
1051 corr_p4 = (jet.energy() > 0. ? corr / jet.energy() : 1.) * jet;
1052 }
1053
1054
1055 if (ii == 0) {
1056 correction += corr_p4;
1057 theSumPt += corr_p4.pt();
1058 theSumEnergy += corr_p4.energy();
1059 }
1060 else if (ii == 1) {
1061 correction -= corr_p4;
1062 theSumResp += corr_p4.energy();
1063 }
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 }
1143
1144
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
1197
1198
1199
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
1280
1281
1282
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 }
1298
1299
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
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 }
1324
1325
1326
1327
1328
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
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
1378 clear();
1379 std::vector<Element> data;
1380
1381
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
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
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
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
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
1464
1465
1466
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
1479
1480
1481
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
1494
1495
1496
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
1509
1510
1511
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) {
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) {
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
1548
1549
1550
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
1690
1691
1692
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 }
