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