File indexing completed on 2023-01-17 23:57:55
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
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
1360
1361 double SquareEtaRingWithoutJets = ja;
1362
1363 EnergyOfBackgroundCharged = EnergyOfBackgroundCharged / SquareEtaRingWithoutJets;
1364 ResponseOfBackgroundCharged = ResponseOfBackgroundCharged / SquareEtaRingWithoutJets;
1365
1366 EnergyOfBackgroundCharged = M_PI * mConeSize * mConeSize * EnergyOfBackgroundCharged;
1367 ResponseOfBackgroundCharged = M_PI * mConeSize * mConeSize * ResponseOfBackgroundCharged;
1368
1369 NewResponse = NewResponse - EnergyOfBackgroundCharged + ResponseOfBackgroundCharged;
1370
1371 mScale = NewResponse / fJet.energy();
1372 if (mScale < 0.)
1373 mScale = 0.;
1374 return mScale;
1375 }
1376
1377
1378 Map::Map(std::string input, bool verbose) : eta_(), pt_(), data_() {
1379
1380 clear();
1381 std::vector<Element> data;
1382
1383
1384 std::string file = edm::FileInPath(input).fullPath();
1385 std::ifstream in(file.c_str());
1386 string line;
1387 uint32_t ieta_old = 0;
1388 while (std::getline(in, line)) {
1389 if (line.empty() || line[0] == '#') {
1390 continue;
1391 }
1392 std::istringstream ss(line);
1393 Element temp;
1394 ss >> temp.ieta_ >> temp.ipt_ >> temp.eta_ >> temp.pt_ >> temp.val_;
1395 data.push_back(temp);
1396 if (!ieta_old || temp.ieta_ != ieta_old) {
1397 if (eta_.size() < temp.ieta_ + 1) {
1398 eta_.resize(temp.ieta_ + 1, 0.);
1399 }
1400 eta_[temp.ieta_] = temp.eta_;
1401 ieta_old = temp.ieta_;
1402 }
1403 if (pt_.size() < temp.ipt_ + 1) {
1404 pt_.resize(temp.ipt_ + 1, 0.);
1405 }
1406 pt_[temp.ipt_] = temp.pt_;
1407 }
1408
1409
1410 data_.resize(eta_.size(), VDouble(pt_.size(), 0.));
1411 std::vector<Element>::const_iterator idata = data.begin();
1412 std::vector<Element>::const_iterator jdata = data.end();
1413 for (; idata != jdata; ++idata) {
1414 data_[idata->ieta_][idata->ipt_] = idata->val_;
1415 }
1416
1417
1418 if (data_.empty() || data_[0].empty()) {
1419 std::stringstream ss;
1420 ss << "[jpt::Map::" << __func__ << "]"
1421 << " Problem parsing map in location \"" << file << "\"! ";
1422 edm::LogError("JetPlusTrackCorrector") << ss.str();
1423 }
1424
1425
1426 if (eta_.size() != data_.size() || pt_.size() != (data_.empty() ? 0 : data_[0].size())) {
1427 std::stringstream ss;
1428 ss << "[jpt::Map::" << __func__ << "]"
1429 << " Discrepancy b/w number of bins!";
1430 edm::LogError("JetPlusTrackCorrector") << ss.str();
1431 }
1432
1433
1434 if (verbose && edm::isDebugEnabled()) {
1435 std::stringstream ss;
1436 ss << "[jpt::Map::" << __func__ << "]"
1437 << " Parsed contents of map at location:" << std::endl
1438 << "\"" << file << "\"" << std::endl;
1439 print(ss);
1440 LogTrace("JetPlusTrackCorrector") << ss.str();
1441 }
1442 }
1443
1444
1445
1446 Map::Map() : eta_(), pt_(), data_() { clear(); }
1447
1448
1449
1450 Map::~Map() { clear(); }
1451
1452
1453
1454 void Map::clear() {
1455 eta_.clear();
1456 pt_.clear();
1457 data_.clear();
1458 }
1459
1460
1461 double Map::eta(uint32_t eta_bin) const {
1462 if (!eta_.empty() && eta_bin < eta_.size()) {
1463 return eta_[eta_bin];
1464 } else {
1465
1466
1467
1468
1469
1470 return eta_[eta_.size() - 1];
1471 }
1472 }
1473
1474
1475
1476 double Map::pt(uint32_t pt_bin) const {
1477 if (!pt_.empty() && pt_bin < pt_.size()) {
1478 return pt_[pt_bin];
1479 } else {
1480
1481
1482
1483
1484
1485 return pt_[pt_.size() - 1];
1486 }
1487 }
1488
1489
1490
1491 double Map::binCenterEta(uint32_t eta_bin) const {
1492 if (!eta_.empty() && eta_bin + 1 < eta_.size()) {
1493 return eta_[eta_bin] + (eta_[eta_bin + 1] - eta_[eta_bin]) / 2.;
1494 } else {
1495
1496
1497
1498
1499
1500 return eta_[eta_.size() - 1];
1501 }
1502 }
1503
1504
1505
1506 double Map::binCenterPt(uint32_t pt_bin) const {
1507 if (!pt_.empty() && pt_bin + 1 < pt_.size()) {
1508 return pt_[pt_bin] + (pt_[pt_bin + 1] - pt_[pt_bin]) / 2.;
1509 } else {
1510
1511
1512
1513
1514
1515 return pt_[pt_.size() - 1];
1516 }
1517 }
1518
1519
1520
1521 uint32_t Map::etaBin(double val) const {
1522 val = fabs(val);
1523 for (uint32_t ieta = 0; ieta < nEtaBins() - 1; ++ieta) {
1524 if (val > eta(ieta) && (ieta + 1 == nEtaBins() || val < eta(ieta + 1))) {
1525 return ieta;
1526 }
1527 }
1528 return nEtaBins();
1529 }
1530
1531
1532
1533 uint32_t Map::ptBin(double val) const {
1534 val = fabs(val);
1535 for (uint32_t ipt = 0; ipt < nPtBins() - 1; ++ipt) {
1536 if (val > pt(ipt) && ((ipt + 1) == nPtBins() || val < pt(ipt + 1))) {
1537 return ipt;
1538 }
1539 }
1540 return nPtBins();
1541 }
1542
1543
1544
1545 double Map::value(uint32_t eta_bin, uint32_t pt_bin) const {
1546 if (eta_bin < data_.size() && pt_bin < (data_.empty() ? 0 : data_[0].size())) {
1547 return data_[eta_bin][pt_bin];
1548 } else {
1549
1550
1551
1552
1553
1554 return 1.;
1555 }
1556 }
1557
1558
1559
1560 void Map::print(std::stringstream& ss) const {
1561 ss << " Number of bins in eta : " << data_.size() << std::endl
1562 << " Number of bins in pt : " << (data_.empty() ? 0 : data_[0].size()) << std::endl;
1563 VVDouble::const_iterator ieta = data_.begin();
1564 VVDouble::const_iterator jeta = data_.end();
1565 for (; ieta != jeta; ++ieta) {
1566 VDouble::const_iterator ipt = ieta->begin();
1567 VDouble::const_iterator jpt = ieta->end();
1568 for (; ipt != jpt; ++ipt) {
1569 uint32_t eta_bin = static_cast<uint32_t>(ieta - data_.begin());
1570 uint32_t pt_bin = static_cast<uint32_t>(ipt - ieta->begin());
1571 ss << " EtaBinNumber: " << eta_bin << " PtBinNumber: " << pt_bin << " EtaValue: " << eta_[eta_bin]
1572 << " PtValue: " << pt_[pt_bin] << " Value: " << data_[eta_bin][pt_bin] << std::endl;
1573 }
1574 }
1575 }
1576
1577
1578
1579 MatchedTracks::MatchedTracks() : inVertexInCalo_(), outOfVertexInCalo_(), inVertexOutOfCalo_() { clear(); }
1580
1581
1582
1583 MatchedTracks::~MatchedTracks() { clear(); }
1584
1585
1586
1587 void MatchedTracks::clear() {
1588 inVertexInCalo_.clear();
1589 outOfVertexInCalo_.clear();
1590 inVertexOutOfCalo_.clear();
1591 }
1592
1593
1594
1595 JetTracks::JetTracks() : vertex_(), caloFace_() { clear(); }
1596
1597
1598
1599 JetTracks::~JetTracks() { clear(); }
1600
1601
1602
1603 void JetTracks::clear() {
1604 vertex_.clear();
1605 caloFace_.clear();
1606 }
1607
1608
1609
1610 Efficiency::Efficiency(const jpt::Map& response, const jpt::Map& efficiency, const jpt::Map& leakage)
1611 : response_(response), efficiency_(efficiency), leakage_(leakage) {
1612 reset();
1613 }
1614
1615
1616
1617 double Efficiency::inConeCorr(uint32_t eta_bin, uint32_t pt_bin) const {
1618 if (check(eta_bin, pt_bin, __func__)) {
1619 return (outOfConeCorr(eta_bin, pt_bin) * leakage_.value(eta_bin, pt_bin) * response_.value(eta_bin, pt_bin));
1620 } else {
1621 return 0.;
1622 }
1623 }
1624
1625
1626
1627 double Efficiency::outOfConeCorr(uint32_t eta_bin, uint32_t pt_bin) const {
1628 if (check(eta_bin, pt_bin, __func__)) {
1629 uint16_t ntrks = nTrks(eta_bin, pt_bin);
1630 double mean = meanE(eta_bin, pt_bin);
1631 double eff = (1. - efficiency_.value(eta_bin, pt_bin)) / efficiency_.value(eta_bin, pt_bin);
1632 if (!ntrks) {
1633 return 0.;
1634 }
1635 return (ntrks * eff * mean);
1636 } else {
1637 return 0.;
1638 }
1639 }
1640
1641
1642
1643 uint16_t Efficiency::nTrks(uint32_t eta_bin, uint32_t pt_bin) const {
1644 if (check(eta_bin, pt_bin, __func__)) {
1645 return data_[eta_bin][pt_bin].first;
1646 } else {
1647 return 0;
1648 }
1649 }
1650
1651
1652
1653 double Efficiency::sumE(uint32_t eta_bin, uint32_t pt_bin) const {
1654 if (check(eta_bin, pt_bin, __func__)) {
1655 return data_[eta_bin][pt_bin].second;
1656 } else {
1657 return 0.;
1658 }
1659 }
1660
1661
1662
1663 double Efficiency::meanE(uint32_t eta_bin, uint32_t pt_bin) const {
1664 if (check(eta_bin, pt_bin, __func__)) {
1665 Pair tmp = data_[eta_bin][pt_bin];
1666 if (tmp.first) {
1667 return tmp.second / tmp.first;
1668 } else {
1669 return 0.;
1670 }
1671 } else {
1672 return 0.;
1673 }
1674 }
1675
1676
1677
1678 void Efficiency::addE(uint32_t eta_bin, uint32_t pt_bin, double energy) {
1679 if (check(eta_bin, pt_bin, __func__)) {
1680 data_[eta_bin][pt_bin].first++;
1681 data_[eta_bin][pt_bin].second += energy;
1682 }
1683 }
1684
1685
1686
1687 bool Efficiency::check(uint32_t eta_bin, uint32_t pt_bin, std::string method) const {
1688 if (eta_bin < data_.size() && pt_bin < (data_.empty() ? 0 : data_[0].size())) {
1689 return true;
1690 } else {
1691
1692
1693
1694
1695
1696 return false;
1697 }
1698 }
1699
1700
1701
1702 void Efficiency::reset() {
1703 data_.clear();
1704 data_.resize(response_.nEtaBins(), VPair(response_.nPtBins(), Pair(0, 0.)));
1705 }
1706
1707
1708
1709 void Efficiency::print() const {
1710 if (edm::isDebugEnabled()) {
1711 std::stringstream ss;
1712 ss << "[jpt::Efficiency::" << __func__ << "]"
1713 << " Contents of maps:" << std::endl;
1714 ss << "Response map: " << std::endl;
1715 response_.print(ss);
1716 ss << "Efficiency map: " << std::endl;
1717 efficiency_.print(ss);
1718 ss << "Leakage map: " << std::endl;
1719 leakage_.print(ss);
1720 LogTrace("JetPlusTrackCorrector") << ss.str();
1721 }
1722 }
