Back to home page

Project CMSSW displayed by LXR
 
 

    

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     // compute iso/miniiso
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_;          // only save cands with pT>pT_cut_
0096     const float pT_cut_noIso_;    // above this pT, don't apply any iso cut
0097     const float pfIsolation_DR_;  // isolation radius
0098     const float pfIsolation_DZ_;  // used in determining if pfcand is from PV or PU
0099     const float absIso_cut_;      // save if ANY of absIso, relIso, or miniRelIso pass the cuts
0100     const float relIso_cut_;
0101     const float miniRelIso_cut_;
0102     const float caloJet_DR_;          // save energy of nearest calojet within caloJet_DR_
0103     const float pflepoverlap_DR_;     // pf lepton overlap radius
0104     const float pflepoverlap_pTmin_;  // pf lepton overlap min pT (only look at PF candidates with pT>pflepoverlap_pTmin_)
0105     const float pcRefNearest_DR_;     // radius for nearest charged packed candidate
0106     const float pcRefNearest_pTmin_;  // min pT for nearest charged packed candidate
0107     const float pfneutralsum_DR_;     // pf lepton overlap radius
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 }  // namespace pat
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   // TrackAssociator parameters
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   // packedPFCandidate collection
0182   edm::Handle<pat::PackedCandidateCollection> pc_h;
0183   iEvent.getByToken(pc_, pc_h);
0184   const pat::PackedCandidateCollection* pc = pc_h.product();
0185 
0186   // lostTracks collection
0187   edm::Handle<pat::PackedCandidateCollection> lt_h;
0188   iEvent.getByToken(lt_, lt_h);
0189   const pat::PackedCandidateCollection* lt = lt_h.product();
0190 
0191   // generalTracks collection
0192   edm::Handle<reco::TrackCollection> gt_h;
0193   iEvent.getByToken(gt_, gt_h);
0194   const reco::TrackCollection* generalTracks = gt_h.product();
0195 
0196   // get the primary vertex
0197   edm::Handle<reco::VertexCollection> pvs;
0198   iEvent.getByToken(pv_, pvs);
0199   const reco::Vertex& pv = (*pvs)[0];
0200 
0201   // generalTracks-->packedPFCandidate association
0202   edm::Handle<edm::Association<pat::PackedCandidateCollection>> gt2pc;
0203   iEvent.getByToken(gt2pc_, gt2pc);
0204 
0205   // generalTracks-->lostTracks association
0206   edm::Handle<edm::Association<pat::PackedCandidateCollection>> gt2lt;
0207   iEvent.getByToken(gt2lt_, gt2lt);
0208 
0209   // packedPFCandidates-->particleFlow(reco::PFCandidate) association
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   // associate generalTracks with their DeDx data (estimator for strip dE/dx)
0217   edm::Handle<edm::ValueMap<reco::DeDxData>> gt2dedxStrip;
0218   iEvent.getByToken(gt2dedxStrip_, gt2dedxStrip);
0219 
0220   // associate generalTracks with their DeDx data (estimator for pixel dE/dx)
0221   edm::Handle<edm::ValueMap<reco::DeDxData>> gt2dedxPixel;
0222   iEvent.getByToken(gt2dedxPixel_, gt2dedxPixel);
0223 
0224   // associate generalTracks with their DeDx hit info (used to estimate pixel dE/dx)
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   //add general tracks
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     // Determine if this general track is associated with anything in packedPFCandidates or lostTracks
0254     // Sometimes, a track gets associated w/ a neutral pfCand.
0255     // In this case, ignore the pfCand and take from lostTracks
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;  //to avoid counting packedPFCands in their own isolation
0265     int ltCandInd;  //to avoid pointing lost track to itself when looking for closest
0266 
0267     // get the four-momentum and charge
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;  //assume pion mass
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     // get the isolation of the track
0304     pat::PFIsolation isolationDR03;
0305     pat::PFIsolation miniIso;
0306     getIsolation(polarP4, pc, pfCandInd, isolationDR03, miniIso);
0307 
0308     // isolation cut
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     // get the rest after the pt/iso cuts. Saves some runtime
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();  //NULL reference
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     // if no dEdx info exists, just store -1
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     // get the associated ecal/hcal detectors
0378     TrackDetMatchInfo trackDetInfo = getTrackDetMatchInfo(iEvent, iSetup, gentk);
0379 
0380     // fill ecal/hcal status vectors
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   // there are some number of pfcandidates with no associated track
0442   // (mostly electrons, with a handful of muons)
0443   // here we find these and store. Track-specific variables get some default values
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     // already counted if it has a track reference in the generalTracks collection
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     // get the isolation of the track
0467     pat::PFIsolation isolationDR03;
0468     pat::PFIsolation miniIso;
0469     getIsolation(polarP4, pc, ipc, isolationDR03, miniIso);
0470 
0471     // isolation cut
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     // fill with default values
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);  // these never have dE/dx hit info, as there's no track
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;      // standard isolation
0563   float chmiso = 0, nhmiso = 0, phmiso = 0, pumiso = 0;  // mini isolation
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)  //don't count itself
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       // charged cands from PV get added to trackIso
0575       if (id == 211 && fromPV)
0576         chiso += pt;
0577       // charged cands not from PV get added to pileup iso
0578       else if (id == 211)
0579         puiso += pt;
0580       // neutral hadron iso
0581       if (id == 130)
0582         nhiso += pt;
0583       // photon iso
0584       if (id == 22)
0585         phiso += pt;
0586     }
0587     // same for mini isolation
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 //get overlap of isolated track with a PF lepton
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)  // exclude neutral candidates
0616       continue;
0617     if (!(fromPV))  // exclude candidates not from PV
0618       continue;
0619     if (pt < pflepoverlap_pTmin_)  // exclude pf candidates w/ pT below threshold
0620       continue;
0621 
0622     float dr = deltaR(p4, pf);
0623     if (dr > pflepoverlap_DR_)  // exclude pf candidates far from isolated track
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 //get ref to nearest pf packed candidate
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)  //don't count itself
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)  // exclude neutral candidates
0654       continue;
0655     if (pt < pcRefNearest_pTmin_)  // exclude candidates w/ pT below threshold
0656       continue;
0657     if (dr > dr_min && dr > dr_min_pu)  // exclude too far candidates
0658       continue;
0659 
0660     if (fromPV) {  // Priority to candidates from PV
0661       if (dr < dr_min) {
0662         dr_min = dr;
0663         pc_ref_idx = int(pf_it - pc->begin());
0664       }
0665     } else {  // Otherwise, store candidate from non-PV if no candidate from PV found
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)  // If no candidate from PV was found, store candidate from non-PV (if found)
0675     pc_ref_idx = pc_ref_idx_pu;
0676 }
0677 
0678 // get PF neutral pT sum around isolated track
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)  //don't count itself
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       // neutral hadron sum
0693       if (id == 130)
0694         nhsum += pt;
0695       // photon iso
0696       if (id == 22)
0697         phsum += pt;
0698     }
0699   }
0700 
0701   nsum = nhsum + phsum;
0702 
0703   return nsum;
0704 }
0705 
0706 // get the estimated DeDx in either the pixels or strips (or both)
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     // probably shouldn't happen
0723     if (!isPixel && !isStrip)
0724       continue;
0725 
0726     // shape selection for strips
0727     if (isStrip && !deDxTools::shapeSelection(*(hitInfo->stripCluster(ih))))
0728       continue;
0729 
0730     float Norm = 0;
0731     if (isPixel)
0732       Norm = 3.61e-06;  //compute the normalization factor to get the energy in MeV/mm
0733     if (isStrip)
0734       Norm = 3.61e-06 * 265;
0735 
0736     //save the dE/dx in MeV/mm to a vector.
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   //build the harmonic 2 dE/dx estimator
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   // can't use the associate() using reco::Track directly, since
0760   // track->extra() is non-null but segfaults when trying to use it
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);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!