File indexing completed on 2024-12-20 03:14:05
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>& gsfEles,
1128 const edm::Handle<RecoElectronIds>& idDecisionMap) const {
1129 if (!gsfEles.isValid()) {
1130 return false;
1131 }
1132
1133 double deltaRMIN = 999.;
1134
1135 for (auto ele = gsfEles->begin(); ele != gsfEles->end(); ++ele) {
1136 const edm::Ptr<reco::GsfElectron> elePtr(gsfEles, ele - gsfEles->begin());
1137 bool passID = false;
1138 if (idDecisionMap.isValid())
1139 passID = (*idDecisionMap)[elePtr];
1140
1141 if (!passID) {
1142 continue;
1143 }
1144
1145
1146 auto dR2 = deltaR2(*ele, **itrk);
1147 if (dR2 < deltaRMIN) {
1148 deltaRMIN = dR2;
1149 }
1150 }
1151 return deltaRMIN < electronDRmatch_ * electronDRmatch_;
1152 }
1153
1154 bool JetPlusTrackCorrector::matchElectrons(TrackRefs::const_iterator& itrk,
1155 const edm::Handle<pat::ElectronCollection>& elecs) const {
1156 if (elecs->empty()) {
1157 return false;
1158 }
1159
1160 double deltaRMIN = 999.;
1161 for (auto const& ielec : *elecs) {
1162 auto dR2 = deltaR2(ielec, **itrk);
1163 if (dR2 < deltaRMIN) {
1164 deltaRMIN = dR2;
1165 }
1166 }
1167 return deltaRMIN < electronDRmatch_ * electronDRmatch_;
1168 }
1169
1170
1171
1172 void JetPlusTrackCorrector::rebuildJta(const reco::Jet& fJet,
1173 const JetTracksAssociations& jtV0,
1174 TrackRefs& tracksthis,
1175 TrackRefs& Excl) const {
1176 tracksthis = reco::JetTracksAssociation::getValue(jtV0, fJet);
1177
1178 if (jetSplitMerge_ < 0)
1179 return;
1180
1181 typedef std::vector<reco::JetBaseRef>::iterator JetBaseRefIterator;
1182 std::vector<reco::JetBaseRef> theJets = reco::JetTracksAssociation::allJets(jtV0);
1183
1184 TrackRefs tracks = tracksthis;
1185 tracksthis.clear();
1186
1187 double jetEta = fJet.eta();
1188 double jetPhi = fJet.phi();
1189 double jetEtIn = 1.0 / fJet.et();
1190
1191 for (TrackRefs::iterator it = tracks.begin(); it != tracks.end(); it++) {
1192 double trkEta = (**it).eta();
1193 double trkPhi = (**it).phi();
1194 double dR2this = deltaR2(jetEta, jetPhi, trkEta, trkPhi);
1195
1196
1197
1198
1199
1200 double scalethis = dR2this;
1201 if (jetSplitMerge_ == 0)
1202 scalethis = 1. * jetEtIn;
1203 if (jetSplitMerge_ == 2)
1204 scalethis = dR2this * jetEtIn;
1205 int flag = 1;
1206 for (JetBaseRefIterator ii = theJets.begin(); ii != theJets.end(); ii++) {
1207 if (&(**ii) == &fJet) {
1208 continue;
1209 }
1210 double dR2 = deltaR2((*ii)->eta(), (*ii)->phi(), trkEta, trkPhi);
1211 double scale = dR2;
1212 if (jetSplitMerge_ == 0)
1213 scale = 1. / (**ii).et();
1214 if (jetSplitMerge_ == 2)
1215 scale = dR2 / (**ii).et();
1216 if (scale < scalethis)
1217 flag = 0;
1218
1219 if (flag == 0) {
1220 break;
1221 }
1222 }
1223
1224 if (flag == 1) {
1225 tracksthis.push_back(*it);
1226 } else {
1227 Excl.push_back(*it);
1228 }
1229 }
1230
1231 return;
1232 }
1233
1234
1235
1236 void JetPlusTrackCorrector::excludeJta(const reco::Jet& fJet,
1237 const JetTracksAssociations& jtV0,
1238 TrackRefs& tracksthis,
1239 const TrackRefs& Excl) const {
1240 tracksthis = reco::JetTracksAssociation::getValue(jtV0, fJet);
1241 if (Excl.empty())
1242 return;
1243 if (jetSplitMerge_ < 0)
1244 return;
1245
1246 TrackRefs tracks = tracksthis;
1247 tracksthis.clear();
1248
1249 for (TrackRefs::iterator it = tracks.begin(); it != tracks.end(); it++) {
1250 TrackRefs::iterator itold = find(Excl.begin(), Excl.end(), (*it));
1251 if (itold == Excl.end()) {
1252 tracksthis.push_back(*it);
1253 }
1254 }
1255
1256 return;
1257 }
1258
1259
1260
1261 double JetPlusTrackCorrector::correctAA(const reco::Jet& fJet,
1262 const reco::TrackRefVector& trBgOutOfVertex,
1263 double& mConeSize,
1264 const reco::TrackRefVector& pioninin,
1265 const reco::TrackRefVector& pioninout,
1266 double ja,
1267 const reco::TrackRefVector& trBgOutOfCalo) const {
1268 double mScale = 1.;
1269 double NewResponse = fJet.energy();
1270
1271 if (trBgOutOfVertex.empty())
1272 return mScale;
1273 double EnergyOfBackgroundCharged = 0.;
1274 double ResponseOfBackgroundCharged = 0.;
1275
1276
1277
1278
1279
1280
1281 for (reco::TrackRefVector::iterator iBgtV = trBgOutOfVertex.begin(); iBgtV != trBgOutOfVertex.end(); iBgtV++) {
1282 double eta = fabs((**iBgtV).eta());
1283 double pt = fabs((**iBgtV).pt());
1284 uint32_t ieta = response_.etaBin(eta);
1285 uint32_t ipt = response_.ptBin(pt);
1286
1287 if (fabs(fJet.eta() - (**iBgtV).eta()) > mConeSize)
1288 continue;
1289
1290 double echarBg = sqrt((**iBgtV).px() * (**iBgtV).px() + (**iBgtV).py() * (**iBgtV).py() +
1291 (**iBgtV).pz() * (**iBgtV).pz() + 0.14 * 0.14);
1292
1293 EnergyOfBackgroundCharged += echarBg / efficiency_.value(ieta, ipt);
1294
1295 }
1296
1297
1298
1299 for (reco::TrackRefVector::iterator iBgtC = trBgOutOfCalo.begin(); iBgtC != trBgOutOfCalo.end(); iBgtC++) {
1300 double eta = fabs((**iBgtC).eta());
1301 double pt = fabs((**iBgtC).pt());
1302 uint32_t ieta = response_.etaBin(eta);
1303 uint32_t ipt = response_.ptBin(pt);
1304
1305 if (fabs(fJet.eta() - (**iBgtC).eta()) > mConeSize)
1306 continue;
1307
1308
1309 if (ieta >= response_.nEtaBins()) {
1310 continue;
1311 }
1312 if (ipt >= response_.nPtBins()) {
1313 ipt = response_.nPtBins() - 1;
1314 }
1315
1316 double echarBg = sqrt((**iBgtC).px() * (**iBgtC).px() + (**iBgtC).py() * (**iBgtC).py() +
1317 (**iBgtC).pz() * (**iBgtC).pz() + 0.14 * 0.14);
1318
1319 ResponseOfBackgroundCharged += echarBg * response_.value(ieta, ipt) / efficiency_.value(ieta, ipt);
1320
1321 }
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
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 }
