File indexing completed on 2022-03-18 02:30:11
0001 #include <string>
0002
0003 #include "FWCore/Framework/interface/stream/EDProducer.h"
0004 #include "FWCore/Framework/interface/ESHandle.h"
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/Framework/interface/EventSetup.h"
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/Utilities/interface/Exception.h"
0010
0011 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0012
0013 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0014 #include "DataFormats/PatCandidates/interface/IsolatedTrack.h"
0015 #include "DataFormats/PatCandidates/interface/PFIsolation.h"
0016 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0017 #include "DataFormats/Common/interface/View.h"
0018 #include "DataFormats/Common/interface/RefToPtr.h"
0019 #include "DataFormats/Math/interface/deltaR.h"
0020 #include "DataFormats/JetReco/interface/CaloJet.h"
0021 #include "DataFormats/TrackReco/interface/Track.h"
0022 #include "DataFormats/TrackReco/interface/TrackExtraFwd.h"
0023 #include "DataFormats/TrackReco/interface/DeDxData.h"
0024 #include "DataFormats/TrackReco/interface/DeDxHitInfo.h"
0025 #include "DataFormats/VertexReco/interface/Vertex.h"
0026 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0027 #include "RecoTracker/DeDx/interface/DeDxTools.h"
0028 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0029
0030 #include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
0031 #include "TrackingTools/TrackAssociator/interface/TrackAssociatorParameters.h"
0032 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0033
0034 #include "CondFormats/HcalObjects/interface/HcalChannelQuality.h"
0035 #include "CondFormats/HcalObjects/interface/HcalChannelStatus.h"
0036 #include "CondFormats/EcalObjects/interface/EcalChannelStatus.h"
0037 #include "CondFormats/DataRecord/interface/HcalChannelQualityRcd.h"
0038 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0039
0040 #include "MagneticField/Engine/interface/MagneticField.h"
0041
0042 namespace pat {
0043
0044 class PATIsolatedTrackProducer : public edm::stream::EDProducer<> {
0045 public:
0046 typedef pat::IsolatedTrack::PolarLorentzVector PolarLorentzVector;
0047 typedef pat::IsolatedTrack::LorentzVector LorentzVector;
0048
0049 explicit PATIsolatedTrackProducer(const edm::ParameterSet&);
0050 ~PATIsolatedTrackProducer() override;
0051
0052 void produce(edm::Event&, const edm::EventSetup&) override;
0053
0054
0055 void getIsolation(const PolarLorentzVector& p4,
0056 const pat::PackedCandidateCollection* pc,
0057 int pc_idx,
0058 pat::PFIsolation& iso,
0059 pat::PFIsolation& miniiso) const;
0060
0061 bool getPFLeptonOverlap(const PolarLorentzVector& p4, const pat::PackedCandidateCollection* pc) const;
0062
0063 float getPFNeutralSum(const PolarLorentzVector& p4, const pat::PackedCandidateCollection* pc, int pc_idx) const;
0064
0065 void getNearestPCRef(const PolarLorentzVector& p4,
0066 const pat::PackedCandidateCollection* pc,
0067 int pc_idx,
0068 int& pc_ref_idx) const;
0069
0070 float getDeDx(const reco::DeDxHitInfo* hitInfo, bool doPixel, bool doStrip) const;
0071
0072 TrackDetMatchInfo getTrackDetMatchInfo(const edm::Event&, const edm::EventSetup&, const reco::Track&);
0073
0074 void getCaloJetEnergy(const PolarLorentzVector&, const reco::CaloJetCollection*, float&, float&) const;
0075
0076 private:
0077 const edm::EDGetTokenT<pat::PackedCandidateCollection> pc_;
0078 const edm::EDGetTokenT<pat::PackedCandidateCollection> lt_;
0079 const edm::EDGetTokenT<reco::TrackCollection> gt_;
0080 const edm::EDGetTokenT<reco::VertexCollection> pv_;
0081 const edm::EDGetTokenT<edm::Association<pat::PackedCandidateCollection>> gt2pc_;
0082 const edm::EDGetTokenT<edm::Association<pat::PackedCandidateCollection>> gt2lt_;
0083 const edm::EDGetTokenT<edm::Association<reco::PFCandidateCollection>> pc2pf_;
0084 const edm::EDGetTokenT<reco::CaloJetCollection> caloJets_;
0085 const edm::EDGetTokenT<edm::ValueMap<reco::DeDxData>> gt2dedxStrip_;
0086 const edm::EDGetTokenT<edm::ValueMap<reco::DeDxData>> gt2dedxPixel_;
0087 const edm::EDGetTokenT<reco::DeDxHitInfoAss> gt2dedxHitInfo_;
0088 const edm::ESGetToken<HcalChannelQuality, HcalChannelQualityRcd> hcalQToken_;
0089 const edm::ESGetToken<EcalChannelStatus, EcalChannelStatusRcd> ecalSToken_;
0090 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> bFieldToken_;
0091 const bool addPrescaledDeDxTracks_;
0092 const edm::EDGetTokenT<edm::ValueMap<int>> gt2dedxHitInfoPrescale_;
0093 const bool usePrecomputedDeDxStrip_;
0094 const bool usePrecomputedDeDxPixel_;
0095 const float pT_cut_;
0096 const float pT_cut_noIso_;
0097 const float pfIsolation_DR_;
0098 const float pfIsolation_DZ_;
0099 const float absIso_cut_;
0100 const float relIso_cut_;
0101 const float miniRelIso_cut_;
0102 const float caloJet_DR_;
0103 const float pflepoverlap_DR_;
0104 const float pflepoverlap_pTmin_;
0105 const float pcRefNearest_DR_;
0106 const float pcRefNearest_pTmin_;
0107 const float pfneutralsum_DR_;
0108 const bool useHighPurity_;
0109 const bool saveDeDxHitInfo_;
0110 StringCutObjectSelector<pat::IsolatedTrack> saveDeDxHitInfoCut_;
0111
0112 std::vector<double> miniIsoParams_;
0113
0114 TrackDetectorAssociator trackAssociator_;
0115 TrackAssociatorParameters trackAssocParameters_;
0116 };
0117 }
0118
0119 pat::PATIsolatedTrackProducer::PATIsolatedTrackProducer(const edm::ParameterSet& iConfig)
0120 : pc_(consumes<pat::PackedCandidateCollection>(iConfig.getParameter<edm::InputTag>("packedPFCandidates"))),
0121 lt_(consumes<pat::PackedCandidateCollection>(iConfig.getParameter<edm::InputTag>("lostTracks"))),
0122 gt_(consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("generalTracks"))),
0123 pv_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("primaryVertices"))),
0124 gt2pc_(consumes<edm::Association<pat::PackedCandidateCollection>>(
0125 iConfig.getParameter<edm::InputTag>("packedPFCandidates"))),
0126 gt2lt_(consumes<edm::Association<pat::PackedCandidateCollection>>(
0127 iConfig.getParameter<edm::InputTag>("lostTracks"))),
0128 pc2pf_(consumes<edm::Association<reco::PFCandidateCollection>>(
0129 iConfig.getParameter<edm::InputTag>("packedPFCandidates"))),
0130 caloJets_(consumes<reco::CaloJetCollection>(iConfig.getParameter<edm::InputTag>("caloJets"))),
0131 gt2dedxStrip_(consumes<edm::ValueMap<reco::DeDxData>>(iConfig.getParameter<edm::InputTag>("dEdxDataStrip"))),
0132 gt2dedxPixel_(consumes<edm::ValueMap<reco::DeDxData>>(iConfig.getParameter<edm::InputTag>("dEdxDataPixel"))),
0133 gt2dedxHitInfo_(consumes<reco::DeDxHitInfoAss>(iConfig.getParameter<edm::InputTag>("dEdxHitInfo"))),
0134 hcalQToken_(esConsumes(edm::ESInputTag("", "withTopo"))),
0135 ecalSToken_(esConsumes()),
0136 bFieldToken_(esConsumes()),
0137 addPrescaledDeDxTracks_(iConfig.getParameter<bool>("addPrescaledDeDxTracks")),
0138 gt2dedxHitInfoPrescale_(addPrescaledDeDxTracks_ ? consumes<edm::ValueMap<int>>(
0139 iConfig.getParameter<edm::InputTag>("dEdxHitInfoPrescale"))
0140 : edm::EDGetTokenT<edm::ValueMap<int>>()),
0141 usePrecomputedDeDxStrip_(iConfig.getParameter<bool>("usePrecomputedDeDxStrip")),
0142 usePrecomputedDeDxPixel_(iConfig.getParameter<bool>("usePrecomputedDeDxPixel")),
0143 pT_cut_(iConfig.getParameter<double>("pT_cut")),
0144 pT_cut_noIso_(iConfig.getParameter<double>("pT_cut_noIso")),
0145 pfIsolation_DR_(iConfig.getParameter<double>("pfIsolation_DR")),
0146 pfIsolation_DZ_(iConfig.getParameter<double>("pfIsolation_DZ")),
0147 absIso_cut_(iConfig.getParameter<double>("absIso_cut")),
0148 relIso_cut_(iConfig.getParameter<double>("relIso_cut")),
0149 miniRelIso_cut_(iConfig.getParameter<double>("miniRelIso_cut")),
0150 caloJet_DR_(iConfig.getParameter<double>("caloJet_DR")),
0151 pflepoverlap_DR_(iConfig.getParameter<double>("pflepoverlap_DR")),
0152 pflepoverlap_pTmin_(iConfig.getParameter<double>("pflepoverlap_pTmin")),
0153 pcRefNearest_DR_(iConfig.getParameter<double>("pcRefNearest_DR")),
0154 pcRefNearest_pTmin_(iConfig.getParameter<double>("pcRefNearest_pTmin")),
0155 pfneutralsum_DR_(iConfig.getParameter<double>("pfneutralsum_DR")),
0156 useHighPurity_(iConfig.getParameter<bool>("useHighPurity")),
0157 saveDeDxHitInfo_(iConfig.getParameter<bool>("saveDeDxHitInfo")),
0158 saveDeDxHitInfoCut_(iConfig.getParameter<std::string>("saveDeDxHitInfoCut")) {
0159
0160 edm::ParameterSet parameters = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
0161 edm::ConsumesCollector iC = consumesCollector();
0162 trackAssocParameters_.loadParameters(parameters, iC);
0163
0164 trackAssociator_.useDefaultPropagator();
0165
0166 miniIsoParams_ = iConfig.getParameter<std::vector<double>>("miniIsoParams");
0167 if (miniIsoParams_.size() != 3)
0168 throw cms::Exception("ParameterError") << "miniIsoParams must have exactly 3 elements.\n";
0169
0170 produces<pat::IsolatedTrackCollection>();
0171
0172 if (saveDeDxHitInfo_) {
0173 produces<reco::DeDxHitInfoCollection>();
0174 produces<reco::DeDxHitInfoAss>();
0175 }
0176 }
0177
0178 pat::PATIsolatedTrackProducer::~PATIsolatedTrackProducer() {}
0179
0180 void pat::PATIsolatedTrackProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0181
0182 edm::Handle<pat::PackedCandidateCollection> pc_h;
0183 iEvent.getByToken(pc_, pc_h);
0184 const pat::PackedCandidateCollection* pc = pc_h.product();
0185
0186
0187 edm::Handle<pat::PackedCandidateCollection> lt_h;
0188 iEvent.getByToken(lt_, lt_h);
0189 const pat::PackedCandidateCollection* lt = lt_h.product();
0190
0191
0192 edm::Handle<reco::TrackCollection> gt_h;
0193 iEvent.getByToken(gt_, gt_h);
0194 const reco::TrackCollection* generalTracks = gt_h.product();
0195
0196
0197 edm::Handle<reco::VertexCollection> pvs;
0198 iEvent.getByToken(pv_, pvs);
0199 const reco::Vertex& pv = (*pvs)[0];
0200
0201
0202 edm::Handle<edm::Association<pat::PackedCandidateCollection>> gt2pc;
0203 iEvent.getByToken(gt2pc_, gt2pc);
0204
0205
0206 edm::Handle<edm::Association<pat::PackedCandidateCollection>> gt2lt;
0207 iEvent.getByToken(gt2lt_, gt2lt);
0208
0209
0210 edm::Handle<edm::Association<reco::PFCandidateCollection>> pc2pf;
0211 iEvent.getByToken(pc2pf_, pc2pf);
0212
0213 edm::Handle<reco::CaloJetCollection> caloJets;
0214 iEvent.getByToken(caloJets_, caloJets);
0215
0216
0217 edm::Handle<edm::ValueMap<reco::DeDxData>> gt2dedxStrip;
0218 iEvent.getByToken(gt2dedxStrip_, gt2dedxStrip);
0219
0220
0221 edm::Handle<edm::ValueMap<reco::DeDxData>> gt2dedxPixel;
0222 iEvent.getByToken(gt2dedxPixel_, gt2dedxPixel);
0223
0224
0225 edm::Handle<reco::DeDxHitInfoAss> gt2dedxHitInfo;
0226 iEvent.getByToken(gt2dedxHitInfo_, gt2dedxHitInfo);
0227 edm::Handle<edm::ValueMap<int>> gt2dedxHitInfoPrescale;
0228 if (addPrescaledDeDxTracks_) {
0229 iEvent.getByToken(gt2dedxHitInfoPrescale_, gt2dedxHitInfoPrescale);
0230 }
0231
0232 const HcalChannelQuality* hcalQ = &iSetup.getData(hcalQToken_);
0233
0234 const EcalChannelStatus* ecalS = &iSetup.getData(ecalSToken_);
0235
0236 auto outDeDxC = std::make_unique<reco::DeDxHitInfoCollection>();
0237 std::vector<int> dEdXass;
0238
0239 auto outPtrP = std::make_unique<std::vector<pat::IsolatedTrack>>();
0240
0241
0242 for (unsigned int igt = 0; igt < generalTracks->size(); igt++) {
0243 const reco::Track& gentk = (*gt_h)[igt];
0244 if (useHighPurity_)
0245 if (!(gentk.quality(reco::TrackBase::qualityByName("highPurity"))))
0246 continue;
0247 reco::TrackRef tkref = reco::TrackRef(gt_h, igt);
0248 pat::PackedCandidateRef pcref = (*gt2pc)[tkref];
0249 pat::PackedCandidateRef ltref = (*gt2lt)[tkref];
0250 const pat::PackedCandidate* pfCand = pcref.get();
0251 const pat::PackedCandidate* lostTrack = ltref.get();
0252
0253
0254
0255
0256 bool isInPackedCands = (pcref.isNonnull() && pcref.id() == pc_h.id() && pfCand->charge() != 0);
0257 bool isInLostTracks = (ltref.isNonnull() && ltref.id() == lt_h.id());
0258
0259 PolarLorentzVector polarP4;
0260 LorentzVector p4;
0261 pat::PackedCandidateRef refToCand;
0262 int pdgId, charge, fromPV;
0263 float dz, dxy, dzError, dxyError;
0264 int pfCandInd;
0265 int ltCandInd;
0266
0267
0268 if (isInPackedCands) {
0269 p4 = pfCand->p4();
0270 polarP4 = pfCand->p4();
0271 charge = pfCand->charge();
0272 pfCandInd = pcref.key();
0273 ltCandInd = -1;
0274 } else if (isInLostTracks) {
0275 p4 = lostTrack->p4();
0276 polarP4 = lostTrack->p4();
0277 charge = lostTrack->charge();
0278 pfCandInd = -1;
0279 ltCandInd = ltref.key();
0280 } else {
0281 double m = 0.13957018;
0282 double E = sqrt(m * m + gentk.p() * gentk.p());
0283 p4.SetPxPyPzE(gentk.px(), gentk.py(), gentk.pz(), E);
0284 polarP4.SetCoordinates(gentk.pt(), gentk.eta(), gentk.phi(), m);
0285 charge = gentk.charge();
0286 pfCandInd = -1;
0287 ltCandInd = -1;
0288 }
0289
0290 int prescaled = 0;
0291 if (addPrescaledDeDxTracks_) {
0292 const auto& dedxRef = (*gt2dedxHitInfo)[tkref];
0293 if (dedxRef.isNonnull()) {
0294 prescaled = (*gt2dedxHitInfoPrescale)[dedxRef];
0295 }
0296 }
0297
0298 if (polarP4.pt() < pT_cut_ && prescaled <= 1)
0299 continue;
0300 if (charge == 0)
0301 continue;
0302
0303
0304 pat::PFIsolation isolationDR03;
0305 pat::PFIsolation miniIso;
0306 getIsolation(polarP4, pc, pfCandInd, isolationDR03, miniIso);
0307
0308
0309 if (polarP4.pt() < pT_cut_noIso_ && prescaled <= 1 &&
0310 !(isolationDR03.chargedHadronIso() < absIso_cut_ ||
0311 isolationDR03.chargedHadronIso() / polarP4.pt() < relIso_cut_ ||
0312 miniIso.chargedHadronIso() / polarP4.pt() < miniRelIso_cut_))
0313 continue;
0314
0315
0316 if (isInPackedCands) {
0317 pdgId = pfCand->pdgId();
0318 dz = pfCand->dz();
0319 dxy = pfCand->dxy();
0320 dzError = pfCand->hasTrackDetails() ? pfCand->dzError() : gentk.dzError();
0321 dxyError = pfCand->hasTrackDetails() ? pfCand->dxyError() : gentk.dxyError();
0322 fromPV = pfCand->fromPV();
0323 refToCand = pcref;
0324 } else if (isInLostTracks) {
0325 pdgId = lostTrack->pdgId();
0326 dz = lostTrack->dz();
0327 dxy = lostTrack->dxy();
0328 dzError = lostTrack->hasTrackDetails() ? lostTrack->dzError() : gentk.dzError();
0329 dxyError = lostTrack->hasTrackDetails() ? lostTrack->dxyError() : gentk.dxyError();
0330 fromPV = lostTrack->fromPV();
0331 refToCand = ltref;
0332 } else {
0333 pdgId = 0;
0334 dz = gentk.dz(pv.position());
0335 dxy = gentk.dxy(pv.position());
0336 dzError = gentk.dzError();
0337 dxyError = gentk.dxyError();
0338 fromPV = -1;
0339 refToCand = pat::PackedCandidateRef();
0340 }
0341
0342 float caloJetEm, caloJetHad;
0343 getCaloJetEnergy(polarP4, caloJets.product(), caloJetEm, caloJetHad);
0344
0345 bool pfLepOverlap = getPFLeptonOverlap(polarP4, pc);
0346 float pfNeutralSum = getPFNeutralSum(polarP4, pc, pfCandInd);
0347
0348 pat::PackedCandidateRef refToNearestPF = pat::PackedCandidateRef();
0349 int refToNearestPF_idx = -1;
0350 getNearestPCRef(polarP4, pc, pfCandInd, refToNearestPF_idx);
0351 if (refToNearestPF_idx != -1)
0352 refToNearestPF = pat::PackedCandidateRef(pc_h, refToNearestPF_idx);
0353
0354 pat::PackedCandidateRef refToNearestLostTrack = pat::PackedCandidateRef();
0355 int refToNearestLostTrack_idx = -1;
0356 getNearestPCRef(polarP4, lt, ltCandInd, refToNearestLostTrack_idx);
0357 if (refToNearestLostTrack_idx != -1)
0358 refToNearestLostTrack = pat::PackedCandidateRef(lt_h, refToNearestLostTrack_idx);
0359
0360
0361 float dEdxPixel = -1, dEdxStrip = -1;
0362 if (usePrecomputedDeDxStrip_ && gt2dedxStrip.isValid() && gt2dedxStrip->contains(tkref.id())) {
0363 dEdxStrip = (*gt2dedxStrip)[tkref].dEdx();
0364 } else if (gt2dedxHitInfo.isValid() && gt2dedxHitInfo->contains(tkref.id())) {
0365 const reco::DeDxHitInfo* hitInfo = (*gt2dedxHitInfo)[tkref].get();
0366 dEdxStrip = getDeDx(hitInfo, false, true);
0367 }
0368 if (usePrecomputedDeDxPixel_ && gt2dedxPixel.isValid() && gt2dedxPixel->contains(tkref.id())) {
0369 dEdxPixel = (*gt2dedxPixel)[tkref].dEdx();
0370 } else if (gt2dedxHitInfo.isValid() && gt2dedxHitInfo->contains(tkref.id())) {
0371 const reco::DeDxHitInfo* hitInfo = (*gt2dedxHitInfo)[tkref].get();
0372 dEdxPixel = getDeDx(hitInfo, true, false);
0373 }
0374
0375 int trackQuality = gentk.qualityMask();
0376
0377
0378 TrackDetMatchInfo trackDetInfo = getTrackDetMatchInfo(iEvent, iSetup, gentk);
0379
0380
0381 std::vector<uint32_t> crossedHcalStatus;
0382 for (auto const& did : trackDetInfo.crossedHcalIds) {
0383 crossedHcalStatus.push_back(hcalQ->getValues(did.rawId())->getValue());
0384 }
0385 std::vector<uint16_t> crossedEcalStatus;
0386 for (auto const& did : trackDetInfo.crossedEcalIds) {
0387 crossedEcalStatus.push_back(ecalS->find(did.rawId())->getStatusCode());
0388 }
0389
0390 int deltaEta = int((trackDetInfo.trkGlobPosAtEcal.eta() - gentk.eta()) / 0.5 * 250);
0391 int deltaPhi = int((trackDetInfo.trkGlobPosAtEcal.phi() - gentk.phi()) / 0.5 * 250);
0392 if (deltaEta < -250)
0393 deltaEta = -250;
0394 if (deltaEta > 250)
0395 deltaEta = 250;
0396 if (deltaPhi < -250)
0397 deltaPhi = -250;
0398 if (deltaPhi > 250)
0399 deltaPhi = 250;
0400
0401 outPtrP->push_back(pat::IsolatedTrack(isolationDR03,
0402 miniIso,
0403 caloJetEm,
0404 caloJetHad,
0405 pfLepOverlap,
0406 pfNeutralSum,
0407 p4,
0408 charge,
0409 pdgId,
0410 dz,
0411 dxy,
0412 dzError,
0413 dxyError,
0414 gentk.hitPattern(),
0415 dEdxStrip,
0416 dEdxPixel,
0417 fromPV,
0418 trackQuality,
0419 crossedEcalStatus,
0420 crossedHcalStatus,
0421 deltaEta,
0422 deltaPhi,
0423 refToCand,
0424 refToNearestPF,
0425 refToNearestLostTrack));
0426 outPtrP->back().setStatus(prescaled);
0427
0428 if (saveDeDxHitInfo_) {
0429 const auto& dedxRef = (*gt2dedxHitInfo)[tkref];
0430 if (saveDeDxHitInfoCut_(outPtrP->back()) && dedxRef.isNonnull()) {
0431 outDeDxC->push_back(*dedxRef);
0432 dEdXass.push_back(outDeDxC->size() - 1);
0433 } else {
0434 dEdXass.push_back(-1);
0435 }
0436 }
0437 }
0438
0439
0440
0441
0442 for (unsigned int ipc = 0; ipc < pc->size(); ipc++) {
0443 const pat::PackedCandidate& pfCand = pc->at(ipc);
0444 pat::PackedCandidateRef pcref = pat::PackedCandidateRef(pc_h, ipc);
0445 reco::PFCandidateRef pfref = (*pc2pf)[pcref];
0446
0447
0448 if (pfref.get()->trackRef().isNonnull() && pfref.get()->trackRef().id() == gt_h.id())
0449 continue;
0450
0451 PolarLorentzVector polarP4;
0452 pat::PackedCandidateRef refToCand;
0453 int pdgId, charge, fromPV;
0454 float dz, dxy, dzError, dxyError;
0455
0456 polarP4 = pfCand.polarP4();
0457 charge = pfCand.charge();
0458
0459 if (polarP4.pt() < pT_cut_)
0460 continue;
0461 if (charge == 0)
0462 continue;
0463
0464
0465 pat::PFIsolation isolationDR03;
0466 pat::PFIsolation miniIso;
0467 getIsolation(polarP4, pc, ipc, isolationDR03, miniIso);
0468
0469
0470 if (polarP4.pt() < pT_cut_noIso_ && !(isolationDR03.chargedHadronIso() < absIso_cut_ ||
0471 isolationDR03.chargedHadronIso() / polarP4.pt() < relIso_cut_ ||
0472 miniIso.chargedHadronIso() / polarP4.pt() < miniRelIso_cut_))
0473 continue;
0474
0475 pdgId = pfCand.pdgId();
0476 dz = pfCand.dz();
0477 dxy = pfCand.dxy();
0478 if (pfCand.hasTrackDetails()) {
0479 dzError = pfCand.dzError();
0480 dxyError = pfCand.dxyError();
0481 } else {
0482 dzError = 0;
0483 dxyError = 0;
0484 }
0485 fromPV = pfCand.fromPV();
0486 refToCand = pcref;
0487
0488 float caloJetEm, caloJetHad;
0489 getCaloJetEnergy(polarP4, caloJets.product(), caloJetEm, caloJetHad);
0490
0491 bool pfLepOverlap = getPFLeptonOverlap(polarP4, pc);
0492 float pfNeutralSum = getPFNeutralSum(polarP4, pc, ipc);
0493
0494 pat::PackedCandidateRef refToNearestPF = pat::PackedCandidateRef();
0495 int refToNearestPF_idx = -1;
0496 getNearestPCRef(polarP4, pc, ipc, refToNearestPF_idx);
0497 if (refToNearestPF_idx != -1)
0498 refToNearestPF = pat::PackedCandidateRef(pc_h, refToNearestPF_idx);
0499
0500 pat::PackedCandidateRef refToNearestLostTrack = pat::PackedCandidateRef();
0501 int refToNearestLostTrack_idx = -1;
0502 getNearestPCRef(polarP4, lt, -1, refToNearestLostTrack_idx);
0503 if (refToNearestLostTrack_idx != -1)
0504 refToNearestLostTrack = pat::PackedCandidateRef(lt_h, refToNearestLostTrack_idx);
0505
0506
0507 reco::HitPattern hp;
0508 float dEdxPixel = -1, dEdxStrip = -1;
0509 int trackQuality = 0;
0510 std::vector<uint16_t> ecalStatus;
0511 std::vector<uint32_t> hcalStatus;
0512 int deltaEta = 0;
0513 int deltaPhi = 0;
0514
0515 outPtrP->push_back(pat::IsolatedTrack(isolationDR03,
0516 miniIso,
0517 caloJetEm,
0518 caloJetHad,
0519 pfLepOverlap,
0520 pfNeutralSum,
0521 pfCand.p4(),
0522 charge,
0523 pdgId,
0524 dz,
0525 dxy,
0526 dzError,
0527 dxyError,
0528 hp,
0529 dEdxStrip,
0530 dEdxPixel,
0531 fromPV,
0532 trackQuality,
0533 ecalStatus,
0534 hcalStatus,
0535 deltaEta,
0536 deltaPhi,
0537 refToCand,
0538 refToNearestPF,
0539 refToNearestLostTrack));
0540
0541 dEdXass.push_back(-1);
0542 }
0543
0544 auto orphHandle = iEvent.put(std::move(outPtrP));
0545 if (saveDeDxHitInfo_) {
0546 auto dedxOH = iEvent.put(std::move(outDeDxC));
0547 auto dedxMatch = std::make_unique<reco::DeDxHitInfoAss>(dedxOH);
0548 reco::DeDxHitInfoAss::Filler filler(*dedxMatch);
0549 filler.insert(orphHandle, dEdXass.begin(), dEdXass.end());
0550 filler.fill();
0551 iEvent.put(std::move(dedxMatch));
0552 }
0553 }
0554
0555 void pat::PATIsolatedTrackProducer::getIsolation(const PolarLorentzVector& p4,
0556 const pat::PackedCandidateCollection* pc,
0557 int pc_idx,
0558 pat::PFIsolation& iso,
0559 pat::PFIsolation& miniiso) const {
0560 float chiso = 0, nhiso = 0, phiso = 0, puiso = 0;
0561 float chmiso = 0, nhmiso = 0, phmiso = 0, pumiso = 0;
0562 float miniDR = std::max(miniIsoParams_[0], std::min(miniIsoParams_[1], miniIsoParams_[2] / p4.pt()));
0563 for (pat::PackedCandidateCollection::const_iterator pf_it = pc->begin(); pf_it != pc->end(); pf_it++) {
0564 if (int(pf_it - pc->begin()) == pc_idx)
0565 continue;
0566 int id = std::abs(pf_it->pdgId());
0567 bool fromPV = (pf_it->fromPV() > 1 || fabs(pf_it->dz()) < pfIsolation_DZ_);
0568 float pt = pf_it->p4().pt();
0569 float dr = deltaR(p4, *pf_it);
0570
0571 if (dr < pfIsolation_DR_) {
0572
0573 if (id == 211 && fromPV)
0574 chiso += pt;
0575
0576 else if (id == 211)
0577 puiso += pt;
0578
0579 if (id == 130)
0580 nhiso += pt;
0581
0582 if (id == 22)
0583 phiso += pt;
0584 }
0585
0586 if (dr < miniDR) {
0587 if (id == 211 && fromPV)
0588 chmiso += pt;
0589 else if (id == 211)
0590 pumiso += pt;
0591 if (id == 130)
0592 nhmiso += pt;
0593 if (id == 22)
0594 phmiso += pt;
0595 }
0596 }
0597
0598 iso = pat::PFIsolation(chiso, nhiso, phiso, puiso);
0599 miniiso = pat::PFIsolation(chmiso, nhmiso, phmiso, pumiso);
0600 }
0601
0602
0603 bool pat::PATIsolatedTrackProducer::getPFLeptonOverlap(const PolarLorentzVector& p4,
0604 const pat::PackedCandidateCollection* pc) const {
0605 bool isOverlap = false;
0606 float dr_min = pflepoverlap_DR_;
0607 int id_drmin = 0;
0608 for (const auto& pf : *pc) {
0609 int id = std::abs(pf.pdgId());
0610 int charge = std::abs(pf.charge());
0611 bool fromPV = (pf.fromPV() > 1 || std::abs(pf.dz()) < pfIsolation_DZ_);
0612 float pt = pf.pt();
0613 if (charge == 0)
0614 continue;
0615 if (!(fromPV))
0616 continue;
0617 if (pt < pflepoverlap_pTmin_)
0618 continue;
0619
0620 float dr = deltaR(p4, pf);
0621 if (dr > pflepoverlap_DR_)
0622 continue;
0623
0624 if (dr < dr_min) {
0625 dr_min = dr;
0626 id_drmin = id;
0627 }
0628 }
0629
0630 if (dr_min < pflepoverlap_DR_ && (id_drmin == 11 || id_drmin == 13))
0631 isOverlap = true;
0632
0633 return isOverlap;
0634 }
0635
0636
0637 void pat::PATIsolatedTrackProducer::getNearestPCRef(const PolarLorentzVector& p4,
0638 const pat::PackedCandidateCollection* pc,
0639 int pc_idx,
0640 int& pc_ref_idx) const {
0641 float dr_min = pcRefNearest_DR_;
0642 float dr_min_pu = pcRefNearest_DR_;
0643 int pc_ref_idx_pu = -1;
0644 for (pat::PackedCandidateCollection::const_iterator pf_it = pc->begin(); pf_it != pc->end(); pf_it++) {
0645 if (int(pf_it - pc->begin()) == pc_idx)
0646 continue;
0647 int charge = std::abs(pf_it->charge());
0648 bool fromPV = (pf_it->fromPV() > 1 || fabs(pf_it->dz()) < pfIsolation_DZ_);
0649 float pt = pf_it->p4().pt();
0650 float dr = deltaR(p4, *pf_it);
0651 if (charge == 0)
0652 continue;
0653 if (pt < pcRefNearest_pTmin_)
0654 continue;
0655 if (dr > dr_min && dr > dr_min_pu)
0656 continue;
0657
0658 if (fromPV) {
0659 if (dr < dr_min) {
0660 dr_min = dr;
0661 pc_ref_idx = int(pf_it - pc->begin());
0662 }
0663 } else {
0664 if (dr < dr_min_pu) {
0665 dr_min_pu = dr;
0666 pc_ref_idx_pu = int(pf_it - pc->begin());
0667 }
0668 }
0669 }
0670
0671 if (pc_ref_idx == -1 &&
0672 pc_ref_idx_pu != -1)
0673 pc_ref_idx = pc_ref_idx_pu;
0674 }
0675
0676
0677 float pat::PATIsolatedTrackProducer::getPFNeutralSum(const PolarLorentzVector& p4,
0678 const pat::PackedCandidateCollection* pc,
0679 int pc_idx) const {
0680 float nsum = 0;
0681 float nhsum = 0, phsum = 0;
0682 for (pat::PackedCandidateCollection::const_iterator pf_it = pc->begin(); pf_it != pc->end(); pf_it++) {
0683 if (int(pf_it - pc->begin()) == pc_idx)
0684 continue;
0685 int id = std::abs(pf_it->pdgId());
0686 float pt = pf_it->p4().pt();
0687 float dr = deltaR(p4, *pf_it);
0688
0689 if (dr < pfneutralsum_DR_) {
0690
0691 if (id == 130)
0692 nhsum += pt;
0693
0694 if (id == 22)
0695 phsum += pt;
0696 }
0697 }
0698
0699 nsum = nhsum + phsum;
0700
0701 return nsum;
0702 }
0703
0704
0705 float pat::PATIsolatedTrackProducer::getDeDx(const reco::DeDxHitInfo* hitInfo, bool doPixel, bool doStrip) const {
0706 if (hitInfo == nullptr) {
0707 return -1;
0708 }
0709
0710 std::vector<float> charge_vec;
0711 for (unsigned int ih = 0; ih < hitInfo->size(); ih++) {
0712 bool isPixel = (hitInfo->pixelCluster(ih) != nullptr);
0713 bool isStrip = (hitInfo->stripCluster(ih) != nullptr);
0714
0715 if (isPixel && !doPixel)
0716 continue;
0717 if (isStrip && !doStrip)
0718 continue;
0719
0720
0721 if (!isPixel && !isStrip)
0722 continue;
0723
0724
0725 if (isStrip && !deDxTools::shapeSelection(*(hitInfo->stripCluster(ih))))
0726 continue;
0727
0728 float Norm = 0;
0729 if (isPixel)
0730 Norm = 3.61e-06;
0731 if (isStrip)
0732 Norm = 3.61e-06 * 265;
0733
0734
0735 charge_vec.push_back(Norm * hitInfo->charge(ih) / hitInfo->pathlength(ih));
0736 }
0737
0738 int size = charge_vec.size();
0739 float result = 0.0;
0740
0741
0742 float expo = -2;
0743 for (int i = 0; i < size; i++) {
0744 result += pow(charge_vec[i], expo);
0745 }
0746 result = (size > 0) ? pow(result / size, 1. / expo) : 0.0;
0747
0748 return result;
0749 }
0750
0751 TrackDetMatchInfo pat::PATIsolatedTrackProducer::getTrackDetMatchInfo(const edm::Event& iEvent,
0752 const edm::EventSetup& iSetup,
0753 const reco::Track& track) {
0754 auto const& bField = iSetup.getData(bFieldToken_);
0755 FreeTrajectoryState initialState = trajectoryStateTransform::initialFreeState(track, &bField);
0756
0757
0758
0759 return trackAssociator_.associate(iEvent, iSetup, trackAssocParameters_, &initialState);
0760 }
0761
0762 void pat::PATIsolatedTrackProducer::getCaloJetEnergy(const PolarLorentzVector& p4,
0763 const reco::CaloJetCollection* cJets,
0764 float& caloJetEm,
0765 float& caloJetHad) const {
0766 float nearestDR = 999;
0767 int ind = -1;
0768 for (unsigned int i = 0; i < cJets->size(); i++) {
0769 float dR = deltaR(cJets->at(i), p4);
0770 if (dR < caloJet_DR_ && dR < nearestDR) {
0771 nearestDR = dR;
0772 ind = i;
0773 }
0774 }
0775
0776 if (ind == -1) {
0777 caloJetEm = 0;
0778 caloJetHad = 0;
0779 } else {
0780 const reco::CaloJet& cJet = cJets->at(ind);
0781 caloJetEm = cJet.emEnergyInEB() + cJet.emEnergyInEE() + cJet.emEnergyInHF();
0782 caloJetHad = cJet.hadEnergyInHB() + cJet.hadEnergyInHE() + cJet.hadEnergyInHF();
0783 }
0784 }
0785
0786 using pat::PATIsolatedTrackProducer;
0787 #include "FWCore/Framework/interface/MakerMacros.h"
0788 DEFINE_FWK_MODULE(PATIsolatedTrackProducer);