Back to home page

Project CMSSW displayed by LXR
 
 

    

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     // 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     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   // there are some number of pfcandidates with no associated track
0440   // (mostly electrons, with a handful of muons)
0441   // here we find these and store. Track-specific variables get some default values
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     // already counted if it has a track reference in the generalTracks collection
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     // get the isolation of the track
0465     pat::PFIsolation isolationDR03;
0466     pat::PFIsolation miniIso;
0467     getIsolation(polarP4, pc, ipc, isolationDR03, miniIso);
0468 
0469     // isolation cut
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     // fill with default values
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);  // these never have dE/dx hit info, as there's no track
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;      // standard isolation
0561   float chmiso = 0, nhmiso = 0, phmiso = 0, pumiso = 0;  // mini isolation
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)  //don't count itself
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       // charged cands from PV get added to trackIso
0573       if (id == 211 && fromPV)
0574         chiso += pt;
0575       // charged cands not from PV get added to pileup iso
0576       else if (id == 211)
0577         puiso += pt;
0578       // neutral hadron iso
0579       if (id == 130)
0580         nhiso += pt;
0581       // photon iso
0582       if (id == 22)
0583         phiso += pt;
0584     }
0585     // same for mini isolation
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 //get overlap of isolated track with a PF lepton
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)  // exclude neutral candidates
0614       continue;
0615     if (!(fromPV))  // exclude candidates not from PV
0616       continue;
0617     if (pt < pflepoverlap_pTmin_)  // exclude pf candidates w/ pT below threshold
0618       continue;
0619 
0620     float dr = deltaR(p4, pf);
0621     if (dr > pflepoverlap_DR_)  // exclude pf candidates far from isolated track
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 //get ref to nearest pf packed candidate
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)  //don't count itself
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)  // exclude neutral candidates
0652       continue;
0653     if (pt < pcRefNearest_pTmin_)  // exclude candidates w/ pT below threshold
0654       continue;
0655     if (dr > dr_min && dr > dr_min_pu)  // exclude too far candidates
0656       continue;
0657 
0658     if (fromPV) {  // Priority to candidates from PV
0659       if (dr < dr_min) {
0660         dr_min = dr;
0661         pc_ref_idx = int(pf_it - pc->begin());
0662       }
0663     } else {  // Otherwise, store candidate from non-PV if no candidate from PV found
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)  // If no candidate from PV was found, store candidate from non-PV (if found)
0673     pc_ref_idx = pc_ref_idx_pu;
0674 }
0675 
0676 // get PF neutral pT sum around isolated track
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)  //don't count itself
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       // neutral hadron sum
0691       if (id == 130)
0692         nhsum += pt;
0693       // photon iso
0694       if (id == 22)
0695         phsum += pt;
0696     }
0697   }
0698 
0699   nsum = nhsum + phsum;
0700 
0701   return nsum;
0702 }
0703 
0704 // get the estimated DeDx in either the pixels or strips (or both)
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     // probably shouldn't happen
0721     if (!isPixel && !isStrip)
0722       continue;
0723 
0724     // shape selection for strips
0725     if (isStrip && !deDxTools::shapeSelection(*(hitInfo->stripCluster(ih))))
0726       continue;
0727 
0728     float Norm = 0;
0729     if (isPixel)
0730       Norm = 3.61e-06;  //compute the normalization factor to get the energy in MeV/mm
0731     if (isStrip)
0732       Norm = 3.61e-06 * 265;
0733 
0734     //save the dE/dx in MeV/mm to a vector.
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   //build the harmonic 2 dE/dx estimator
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   // can't use the associate() using reco::Track directly, since
0758   // track->extra() is non-null but segfaults when trying to use it
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);
[ 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!