File indexing completed on 2024-04-06 12:23:49
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 crossedHcalStatus.reserve(trackDetInfo.crossedHcalIds.size());
0383 for (auto const& did : trackDetInfo.crossedHcalIds) {
0384 crossedHcalStatus.push_back(hcalQ->getValues(did.rawId())->getValue());
0385 }
0386 std::vector<uint16_t> crossedEcalStatus;
0387 crossedEcalStatus.reserve(trackDetInfo.crossedEcalIds.size());
0388 for (auto const& did : trackDetInfo.crossedEcalIds) {
0389 crossedEcalStatus.push_back(ecalS->find(did.rawId())->getStatusCode());
0390 }
0391
0392 int deltaEta = int((trackDetInfo.trkGlobPosAtEcal.eta() - gentk.eta()) / 0.5 * 250);
0393 int deltaPhi = int((trackDetInfo.trkGlobPosAtEcal.phi() - gentk.phi()) / 0.5 * 250);
0394 if (deltaEta < -250)
0395 deltaEta = -250;
0396 if (deltaEta > 250)
0397 deltaEta = 250;
0398 if (deltaPhi < -250)
0399 deltaPhi = -250;
0400 if (deltaPhi > 250)
0401 deltaPhi = 250;
0402
0403 outPtrP->push_back(pat::IsolatedTrack(isolationDR03,
0404 miniIso,
0405 caloJetEm,
0406 caloJetHad,
0407 pfLepOverlap,
0408 pfNeutralSum,
0409 p4,
0410 charge,
0411 pdgId,
0412 dz,
0413 dxy,
0414 dzError,
0415 dxyError,
0416 gentk.hitPattern(),
0417 dEdxStrip,
0418 dEdxPixel,
0419 fromPV,
0420 trackQuality,
0421 crossedEcalStatus,
0422 crossedHcalStatus,
0423 deltaEta,
0424 deltaPhi,
0425 refToCand,
0426 refToNearestPF,
0427 refToNearestLostTrack));
0428 outPtrP->back().setStatus(prescaled);
0429
0430 if (saveDeDxHitInfo_) {
0431 const auto& dedxRef = (*gt2dedxHitInfo)[tkref];
0432 if (saveDeDxHitInfoCut_(outPtrP->back()) && dedxRef.isNonnull()) {
0433 outDeDxC->push_back(*dedxRef);
0434 dEdXass.push_back(outDeDxC->size() - 1);
0435 } else {
0436 dEdXass.push_back(-1);
0437 }
0438 }
0439 }
0440
0441
0442
0443
0444 for (unsigned int ipc = 0; ipc < pc->size(); ipc++) {
0445 const pat::PackedCandidate& pfCand = pc->at(ipc);
0446 pat::PackedCandidateRef pcref = pat::PackedCandidateRef(pc_h, ipc);
0447 reco::PFCandidateRef pfref = (*pc2pf)[pcref];
0448
0449
0450 if (pfref.get()->trackRef().isNonnull() && pfref.get()->trackRef().id() == gt_h.id())
0451 continue;
0452
0453 PolarLorentzVector polarP4;
0454 pat::PackedCandidateRef refToCand;
0455 int pdgId, charge, fromPV;
0456 float dz, dxy, dzError, dxyError;
0457
0458 polarP4 = pfCand.polarP4();
0459 charge = pfCand.charge();
0460
0461 if (polarP4.pt() < pT_cut_)
0462 continue;
0463 if (charge == 0)
0464 continue;
0465
0466
0467 pat::PFIsolation isolationDR03;
0468 pat::PFIsolation miniIso;
0469 getIsolation(polarP4, pc, ipc, isolationDR03, miniIso);
0470
0471
0472 if (polarP4.pt() < pT_cut_noIso_ && !(isolationDR03.chargedHadronIso() < absIso_cut_ ||
0473 isolationDR03.chargedHadronIso() / polarP4.pt() < relIso_cut_ ||
0474 miniIso.chargedHadronIso() / polarP4.pt() < miniRelIso_cut_))
0475 continue;
0476
0477 pdgId = pfCand.pdgId();
0478 dz = pfCand.dz();
0479 dxy = pfCand.dxy();
0480 if (pfCand.hasTrackDetails()) {
0481 dzError = pfCand.dzError();
0482 dxyError = pfCand.dxyError();
0483 } else {
0484 dzError = 0;
0485 dxyError = 0;
0486 }
0487 fromPV = pfCand.fromPV();
0488 refToCand = pcref;
0489
0490 float caloJetEm, caloJetHad;
0491 getCaloJetEnergy(polarP4, caloJets.product(), caloJetEm, caloJetHad);
0492
0493 bool pfLepOverlap = getPFLeptonOverlap(polarP4, pc);
0494 float pfNeutralSum = getPFNeutralSum(polarP4, pc, ipc);
0495
0496 pat::PackedCandidateRef refToNearestPF = pat::PackedCandidateRef();
0497 int refToNearestPF_idx = -1;
0498 getNearestPCRef(polarP4, pc, ipc, refToNearestPF_idx);
0499 if (refToNearestPF_idx != -1)
0500 refToNearestPF = pat::PackedCandidateRef(pc_h, refToNearestPF_idx);
0501
0502 pat::PackedCandidateRef refToNearestLostTrack = pat::PackedCandidateRef();
0503 int refToNearestLostTrack_idx = -1;
0504 getNearestPCRef(polarP4, lt, -1, refToNearestLostTrack_idx);
0505 if (refToNearestLostTrack_idx != -1)
0506 refToNearestLostTrack = pat::PackedCandidateRef(lt_h, refToNearestLostTrack_idx);
0507
0508
0509 reco::HitPattern hp;
0510 float dEdxPixel = -1, dEdxStrip = -1;
0511 int trackQuality = 0;
0512 std::vector<uint16_t> ecalStatus;
0513 std::vector<uint32_t> hcalStatus;
0514 int deltaEta = 0;
0515 int deltaPhi = 0;
0516
0517 outPtrP->push_back(pat::IsolatedTrack(isolationDR03,
0518 miniIso,
0519 caloJetEm,
0520 caloJetHad,
0521 pfLepOverlap,
0522 pfNeutralSum,
0523 pfCand.p4(),
0524 charge,
0525 pdgId,
0526 dz,
0527 dxy,
0528 dzError,
0529 dxyError,
0530 hp,
0531 dEdxStrip,
0532 dEdxPixel,
0533 fromPV,
0534 trackQuality,
0535 ecalStatus,
0536 hcalStatus,
0537 deltaEta,
0538 deltaPhi,
0539 refToCand,
0540 refToNearestPF,
0541 refToNearestLostTrack));
0542
0543 dEdXass.push_back(-1);
0544 }
0545
0546 auto orphHandle = iEvent.put(std::move(outPtrP));
0547 if (saveDeDxHitInfo_) {
0548 auto dedxOH = iEvent.put(std::move(outDeDxC));
0549 auto dedxMatch = std::make_unique<reco::DeDxHitInfoAss>(dedxOH);
0550 reco::DeDxHitInfoAss::Filler filler(*dedxMatch);
0551 filler.insert(orphHandle, dEdXass.begin(), dEdXass.end());
0552 filler.fill();
0553 iEvent.put(std::move(dedxMatch));
0554 }
0555 }
0556
0557 void pat::PATIsolatedTrackProducer::getIsolation(const PolarLorentzVector& p4,
0558 const pat::PackedCandidateCollection* pc,
0559 int pc_idx,
0560 pat::PFIsolation& iso,
0561 pat::PFIsolation& miniiso) const {
0562 float chiso = 0, nhiso = 0, phiso = 0, puiso = 0;
0563 float chmiso = 0, nhmiso = 0, phmiso = 0, pumiso = 0;
0564 float miniDR = std::max(miniIsoParams_[0], std::min(miniIsoParams_[1], miniIsoParams_[2] / p4.pt()));
0565 for (pat::PackedCandidateCollection::const_iterator pf_it = pc->begin(); pf_it != pc->end(); pf_it++) {
0566 if (int(pf_it - pc->begin()) == pc_idx)
0567 continue;
0568 int id = std::abs(pf_it->pdgId());
0569 bool fromPV = (pf_it->fromPV() > 1 || fabs(pf_it->dz()) < pfIsolation_DZ_);
0570 float pt = pf_it->p4().pt();
0571 float dr = deltaR(p4, *pf_it);
0572
0573 if (dr < pfIsolation_DR_) {
0574
0575 if (id == 211 && fromPV)
0576 chiso += pt;
0577
0578 else if (id == 211)
0579 puiso += pt;
0580
0581 if (id == 130)
0582 nhiso += pt;
0583
0584 if (id == 22)
0585 phiso += pt;
0586 }
0587
0588 if (dr < miniDR) {
0589 if (id == 211 && fromPV)
0590 chmiso += pt;
0591 else if (id == 211)
0592 pumiso += pt;
0593 if (id == 130)
0594 nhmiso += pt;
0595 if (id == 22)
0596 phmiso += pt;
0597 }
0598 }
0599
0600 iso = pat::PFIsolation(chiso, nhiso, phiso, puiso);
0601 miniiso = pat::PFIsolation(chmiso, nhmiso, phmiso, pumiso);
0602 }
0603
0604
0605 bool pat::PATIsolatedTrackProducer::getPFLeptonOverlap(const PolarLorentzVector& p4,
0606 const pat::PackedCandidateCollection* pc) const {
0607 bool isOverlap = false;
0608 float dr_min = pflepoverlap_DR_;
0609 int id_drmin = 0;
0610 for (const auto& pf : *pc) {
0611 int id = std::abs(pf.pdgId());
0612 int charge = std::abs(pf.charge());
0613 bool fromPV = (pf.fromPV() > 1 || std::abs(pf.dz()) < pfIsolation_DZ_);
0614 float pt = pf.pt();
0615 if (charge == 0)
0616 continue;
0617 if (!(fromPV))
0618 continue;
0619 if (pt < pflepoverlap_pTmin_)
0620 continue;
0621
0622 float dr = deltaR(p4, pf);
0623 if (dr > pflepoverlap_DR_)
0624 continue;
0625
0626 if (dr < dr_min) {
0627 dr_min = dr;
0628 id_drmin = id;
0629 }
0630 }
0631
0632 if (dr_min < pflepoverlap_DR_ && (id_drmin == 11 || id_drmin == 13))
0633 isOverlap = true;
0634
0635 return isOverlap;
0636 }
0637
0638
0639 void pat::PATIsolatedTrackProducer::getNearestPCRef(const PolarLorentzVector& p4,
0640 const pat::PackedCandidateCollection* pc,
0641 int pc_idx,
0642 int& pc_ref_idx) const {
0643 float dr_min = pcRefNearest_DR_;
0644 float dr_min_pu = pcRefNearest_DR_;
0645 int pc_ref_idx_pu = -1;
0646 for (pat::PackedCandidateCollection::const_iterator pf_it = pc->begin(); pf_it != pc->end(); pf_it++) {
0647 if (int(pf_it - pc->begin()) == pc_idx)
0648 continue;
0649 int charge = std::abs(pf_it->charge());
0650 bool fromPV = (pf_it->fromPV() > 1 || fabs(pf_it->dz()) < pfIsolation_DZ_);
0651 float pt = pf_it->p4().pt();
0652 float dr = deltaR(p4, *pf_it);
0653 if (charge == 0)
0654 continue;
0655 if (pt < pcRefNearest_pTmin_)
0656 continue;
0657 if (dr > dr_min && dr > dr_min_pu)
0658 continue;
0659
0660 if (fromPV) {
0661 if (dr < dr_min) {
0662 dr_min = dr;
0663 pc_ref_idx = int(pf_it - pc->begin());
0664 }
0665 } else {
0666 if (dr < dr_min_pu) {
0667 dr_min_pu = dr;
0668 pc_ref_idx_pu = int(pf_it - pc->begin());
0669 }
0670 }
0671 }
0672
0673 if (pc_ref_idx == -1 &&
0674 pc_ref_idx_pu != -1)
0675 pc_ref_idx = pc_ref_idx_pu;
0676 }
0677
0678
0679 float pat::PATIsolatedTrackProducer::getPFNeutralSum(const PolarLorentzVector& p4,
0680 const pat::PackedCandidateCollection* pc,
0681 int pc_idx) const {
0682 float nsum = 0;
0683 float nhsum = 0, phsum = 0;
0684 for (pat::PackedCandidateCollection::const_iterator pf_it = pc->begin(); pf_it != pc->end(); pf_it++) {
0685 if (int(pf_it - pc->begin()) == pc_idx)
0686 continue;
0687 int id = std::abs(pf_it->pdgId());
0688 float pt = pf_it->p4().pt();
0689 float dr = deltaR(p4, *pf_it);
0690
0691 if (dr < pfneutralsum_DR_) {
0692
0693 if (id == 130)
0694 nhsum += pt;
0695
0696 if (id == 22)
0697 phsum += pt;
0698 }
0699 }
0700
0701 nsum = nhsum + phsum;
0702
0703 return nsum;
0704 }
0705
0706
0707 float pat::PATIsolatedTrackProducer::getDeDx(const reco::DeDxHitInfo* hitInfo, bool doPixel, bool doStrip) const {
0708 if (hitInfo == nullptr) {
0709 return -1;
0710 }
0711
0712 std::vector<float> charge_vec;
0713 for (unsigned int ih = 0; ih < hitInfo->size(); ih++) {
0714 bool isPixel = (hitInfo->pixelCluster(ih) != nullptr);
0715 bool isStrip = (hitInfo->stripCluster(ih) != nullptr);
0716
0717 if (isPixel && !doPixel)
0718 continue;
0719 if (isStrip && !doStrip)
0720 continue;
0721
0722
0723 if (!isPixel && !isStrip)
0724 continue;
0725
0726
0727 if (isStrip && !deDxTools::shapeSelection(*(hitInfo->stripCluster(ih))))
0728 continue;
0729
0730 float Norm = 0;
0731 if (isPixel)
0732 Norm = 3.61e-06;
0733 if (isStrip)
0734 Norm = 3.61e-06 * 265;
0735
0736
0737 charge_vec.push_back(Norm * hitInfo->charge(ih) / hitInfo->pathlength(ih));
0738 }
0739
0740 int size = charge_vec.size();
0741 float result = 0.0;
0742
0743
0744 float expo = -2;
0745 for (int i = 0; i < size; i++) {
0746 result += pow(charge_vec[i], expo);
0747 }
0748 result = (size > 0) ? pow(result / size, 1. / expo) : 0.0;
0749
0750 return result;
0751 }
0752
0753 TrackDetMatchInfo pat::PATIsolatedTrackProducer::getTrackDetMatchInfo(const edm::Event& iEvent,
0754 const edm::EventSetup& iSetup,
0755 const reco::Track& track) {
0756 auto const& bField = iSetup.getData(bFieldToken_);
0757 FreeTrajectoryState initialState = trajectoryStateTransform::initialFreeState(track, &bField);
0758
0759
0760
0761 return trackAssociator_.associate(iEvent, iSetup, trackAssocParameters_, &initialState);
0762 }
0763
0764 void pat::PATIsolatedTrackProducer::getCaloJetEnergy(const PolarLorentzVector& p4,
0765 const reco::CaloJetCollection* cJets,
0766 float& caloJetEm,
0767 float& caloJetHad) const {
0768 float nearestDR = 999;
0769 int ind = -1;
0770 for (unsigned int i = 0; i < cJets->size(); i++) {
0771 float dR = deltaR(cJets->at(i), p4);
0772 if (dR < caloJet_DR_ && dR < nearestDR) {
0773 nearestDR = dR;
0774 ind = i;
0775 }
0776 }
0777
0778 if (ind == -1) {
0779 caloJetEm = 0;
0780 caloJetHad = 0;
0781 } else {
0782 const reco::CaloJet& cJet = cJets->at(ind);
0783 caloJetEm = cJet.emEnergyInEB() + cJet.emEnergyInEE() + cJet.emEnergyInHF();
0784 caloJetHad = cJet.hadEnergyInHB() + cJet.hadEnergyInHE() + cJet.hadEnergyInHF();
0785 }
0786 }
0787
0788 using pat::PATIsolatedTrackProducer;
0789 #include "FWCore/Framework/interface/MakerMacros.h"
0790 DEFINE_FWK_MODULE(PATIsolatedTrackProducer);