File indexing completed on 2024-04-06 12:23:50
0001 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0002 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0003 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0004 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0005 #include "DataFormats/Common/interface/Association.h"
0006 #include "DataFormats/MuonReco/interface/Muon.h"
0007 #include "DataFormats/VertexReco/interface/Vertex.h"
0008 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0009 #include "FWCore/Framework/interface/global/EDProducer.h"
0010 #include "FWCore/Framework/interface/ESHandle.h"
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0014
0015 namespace {
0016 bool passesQuality(const reco::Track& trk, const std::vector<reco::TrackBase::TrackQuality>& allowedQuals) {
0017 for (const auto& qual : allowedQuals) {
0018 if (trk.quality(qual))
0019 return true;
0020 }
0021 return false;
0022 }
0023 }
0024
0025 namespace pat {
0026 class PATLostTracks : public edm::global::EDProducer<> {
0027 public:
0028 explicit PATLostTracks(const edm::ParameterSet&);
0029 ~PATLostTracks() override;
0030
0031 void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0032
0033 private:
0034 enum class TrkStatus { NOTUSED = 0, PFCAND, PFCANDNOTRKPROPS, PFELECTRON, PFPOSITRON, VTX };
0035 bool passTrkCuts(const reco::Track& tr) const;
0036 void addPackedCandidate(std::vector<pat::PackedCandidate>& cands,
0037 const reco::TrackRef& trk,
0038 const reco::VertexRef& pvSlimmed,
0039 const reco::VertexRefProd& pvSlimmedColl,
0040 const TrkStatus& trkStatus,
0041 const pat::PackedCandidate::PVAssociationQuality& pvAssocQuality,
0042 edm::Handle<reco::MuonCollection> muons) const;
0043 std::pair<int, pat::PackedCandidate::PVAssociationQuality> associateTrkToVtx(const reco::VertexCollection& vertices,
0044 const reco::TrackRef& trk) const;
0045
0046 private:
0047 const edm::EDGetTokenT<reco::PFCandidateCollection> cands_;
0048 const edm::EDGetTokenT<edm::Association<pat::PackedCandidateCollection>> map_;
0049 const edm::EDGetTokenT<reco::TrackCollection> tracks_;
0050 const edm::EDGetTokenT<reco::VertexCollection> vertices_;
0051 const edm::EDGetTokenT<reco::VertexCompositeCandidateCollection> kshorts_;
0052 const edm::EDGetTokenT<reco::VertexCompositeCandidateCollection> lambdas_;
0053 const edm::EDGetTokenT<reco::VertexCollection> pv_;
0054 const edm::EDGetTokenT<reco::VertexCollection> pvOrigs_;
0055 const double minPt_;
0056 const double minHits_;
0057 const double minPixelHits_;
0058 const double minPtToStoreProps_;
0059 const double minPtToStoreLowQualityProps_;
0060 const int covarianceVersion_;
0061 const std::vector<int> covariancePackingSchemas_;
0062 std::vector<reco::TrackBase::TrackQuality> qualsToAutoAccept_;
0063 const edm::EDGetTokenT<reco::MuonCollection> muons_;
0064 StringCutObjectSelector<reco::Track, false> passThroughCut_;
0065 const double maxDzForPrimaryAssignment_;
0066 const double maxDzSigForPrimaryAssignment_;
0067 const double maxDzErrorForPrimaryAssignment_;
0068 const double maxDxyForNotReconstructedPrimary_;
0069 const double maxDxySigForNotReconstructedPrimary_;
0070 const bool useLegacySetup_;
0071 const bool xiSelection_;
0072 const double xiMassCut_;
0073 };
0074 }
0075
0076 pat::PATLostTracks::PATLostTracks(const edm::ParameterSet& iConfig)
0077 : cands_(consumes<reco::PFCandidateCollection>(iConfig.getParameter<edm::InputTag>("inputCandidates"))),
0078 map_(consumes<edm::Association<pat::PackedCandidateCollection>>(
0079 iConfig.getParameter<edm::InputTag>("packedPFCandidates"))),
0080 tracks_(consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("inputTracks"))),
0081 vertices_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("secondaryVertices"))),
0082 kshorts_(consumes<reco::VertexCompositeCandidateCollection>(iConfig.getParameter<edm::InputTag>("kshorts"))),
0083 lambdas_(consumes<reco::VertexCompositeCandidateCollection>(iConfig.getParameter<edm::InputTag>("lambdas"))),
0084 pv_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("primaryVertices"))),
0085 pvOrigs_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("originalVertices"))),
0086 minPt_(iConfig.getParameter<double>("minPt")),
0087 minHits_(iConfig.getParameter<uint32_t>("minHits")),
0088 minPixelHits_(iConfig.getParameter<uint32_t>("minPixelHits")),
0089 minPtToStoreProps_(iConfig.getParameter<double>("minPtToStoreProps")),
0090 minPtToStoreLowQualityProps_(iConfig.getParameter<double>("minPtToStoreLowQualityProps")),
0091 covarianceVersion_(iConfig.getParameter<int>("covarianceVersion")),
0092 covariancePackingSchemas_(iConfig.getParameter<std::vector<int>>("covariancePackingSchemas")),
0093 muons_(consumes<reco::MuonCollection>(iConfig.getParameter<edm::InputTag>("muons"))),
0094 passThroughCut_(iConfig.getParameter<std::string>("passThroughCut")),
0095 maxDzForPrimaryAssignment_(
0096 iConfig.getParameter<edm::ParameterSet>("pvAssignment").getParameter<double>("maxDzForPrimaryAssignment")),
0097 maxDzSigForPrimaryAssignment_(
0098 iConfig.getParameter<edm::ParameterSet>("pvAssignment").getParameter<double>("maxDzSigForPrimaryAssignment")),
0099 maxDzErrorForPrimaryAssignment_(iConfig.getParameter<edm::ParameterSet>("pvAssignment")
0100 .getParameter<double>("maxDzErrorForPrimaryAssignment")),
0101 maxDxyForNotReconstructedPrimary_(iConfig.getParameter<edm::ParameterSet>("pvAssignment")
0102 .getParameter<double>("maxDxyForNotReconstructedPrimary")),
0103 maxDxySigForNotReconstructedPrimary_(iConfig.getParameter<edm::ParameterSet>("pvAssignment")
0104 .getParameter<double>("maxDxySigForNotReconstructedPrimary")),
0105 useLegacySetup_(iConfig.getParameter<bool>("useLegacySetup")),
0106 xiSelection_(iConfig.getParameter<bool>("xiSelection")),
0107 xiMassCut_(iConfig.getParameter<double>("xiMassCut")) {
0108 std::vector<std::string> trkQuals(iConfig.getParameter<std::vector<std::string>>("qualsToAutoAccept"));
0109 std::transform(
0110 trkQuals.begin(), trkQuals.end(), std::back_inserter(qualsToAutoAccept_), reco::TrackBase::qualityByName);
0111
0112 if (std::find(qualsToAutoAccept_.begin(), qualsToAutoAccept_.end(), reco::TrackBase::undefQuality) !=
0113 qualsToAutoAccept_.end()) {
0114 std::ostringstream msg;
0115 msg << " PATLostTracks has a quality requirement which resolves to undefQuality. This usually means a typo and is "
0116 "therefore treated a config error\nquality requirements:\n ";
0117 for (const auto& trkQual : trkQuals)
0118 msg << trkQual << " ";
0119 throw cms::Exception("Configuration") << msg.str();
0120 }
0121
0122 produces<std::vector<reco::Track>>();
0123 produces<std::vector<pat::PackedCandidate>>();
0124 produces<std::vector<pat::PackedCandidate>>("eleTracks");
0125 produces<edm::Association<pat::PackedCandidateCollection>>();
0126 }
0127
0128 pat::PATLostTracks::~PATLostTracks() {}
0129
0130 void pat::PATLostTracks::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0131 edm::Handle<reco::PFCandidateCollection> cands;
0132 iEvent.getByToken(cands_, cands);
0133
0134 edm::Handle<edm::Association<pat::PackedCandidateCollection>> pf2pc;
0135 iEvent.getByToken(map_, pf2pc);
0136
0137 edm::Handle<reco::TrackCollection> tracks;
0138 iEvent.getByToken(tracks_, tracks);
0139
0140 edm::Handle<reco::VertexCollection> vertices;
0141 iEvent.getByToken(vertices_, vertices);
0142
0143 edm::Handle<reco::MuonCollection> muons;
0144 iEvent.getByToken(muons_, muons);
0145
0146 edm::Handle<reco::VertexCompositeCandidateCollection> kshorts;
0147 iEvent.getByToken(kshorts_, kshorts);
0148 edm::Handle<reco::VertexCompositeCandidateCollection> lambdas;
0149 iEvent.getByToken(lambdas_, lambdas);
0150
0151 edm::Handle<reco::VertexCollection> pvs;
0152 iEvent.getByToken(pv_, pvs);
0153 reco::VertexRef pv(pvs.id());
0154 reco::VertexRefProd pvRefProd(pvs);
0155 edm::Handle<reco::VertexCollection> pvOrigs;
0156 iEvent.getByToken(pvOrigs_, pvOrigs);
0157
0158 auto outPtrTrks = std::make_unique<std::vector<reco::Track>>();
0159 auto outPtrTrksAsCands = std::make_unique<std::vector<pat::PackedCandidate>>();
0160 auto outPtrEleTrksAsCands = std::make_unique<std::vector<pat::PackedCandidate>>();
0161
0162 std::vector<TrkStatus> trkStatus(tracks->size(), TrkStatus::NOTUSED);
0163
0164
0165
0166
0167 for (unsigned int ic = 0, nc = cands->size(); ic < nc; ++ic) {
0168 edm::Ref<reco::PFCandidateCollection> r(cands, ic);
0169 const reco::PFCandidate& cand = (*cands)[ic];
0170 if (cand.charge() && cand.trackRef().isNonnull() && cand.trackRef().id() == tracks.id()) {
0171 if (cand.pdgId() == 11)
0172 trkStatus[cand.trackRef().key()] = TrkStatus::PFELECTRON;
0173 else if (cand.pdgId() == -11)
0174 trkStatus[cand.trackRef().key()] = TrkStatus::PFPOSITRON;
0175 else if ((*pf2pc)[r]->numberOfHits() > 0)
0176 trkStatus[cand.trackRef().key()] = TrkStatus::PFCAND;
0177 else
0178 trkStatus[cand.trackRef().key()] = TrkStatus::PFCANDNOTRKPROPS;
0179 }
0180 }
0181
0182
0183 for (const auto& secVert : *vertices) {
0184 for (auto trkIt = secVert.tracks_begin(); trkIt != secVert.tracks_end(); trkIt++) {
0185 if (trkStatus[trkIt->key()] == TrkStatus::NOTUSED)
0186 trkStatus[trkIt->key()] = TrkStatus::VTX;
0187 }
0188 }
0189 for (const auto& v0 : *kshorts) {
0190 for (size_t dIdx = 0; dIdx < v0.numberOfDaughters(); dIdx++) {
0191 size_t key = (dynamic_cast<const reco::RecoChargedCandidate*>(v0.daughter(dIdx)))->track().key();
0192 if (trkStatus[key] == TrkStatus::NOTUSED)
0193 trkStatus[key] = TrkStatus::VTX;
0194 }
0195 }
0196 for (const auto& v0 : *lambdas) {
0197 double protonCharge = 0;
0198 for (size_t dIdx = 0; dIdx < v0.numberOfDaughters(); dIdx++) {
0199 size_t key = (dynamic_cast<const reco::RecoChargedCandidate*>(v0.daughter(dIdx)))->track().key();
0200 if (trkStatus[key] == TrkStatus::NOTUSED)
0201 trkStatus[key] = TrkStatus::VTX;
0202 protonCharge += v0.daughter(dIdx)->charge() * v0.daughter(dIdx)->momentum().mag2();
0203 }
0204 if (xiSelection_) {
0205
0206 math::XYZTLorentzVector p4Lambda(v0.px(), v0.py(), v0.pz(), sqrt(v0.momentum().mag2() + v0.mass() * v0.mass()));
0207
0208 for (unsigned int trkIndx = 0; trkIndx < tracks->size(); trkIndx++) {
0209 reco::TrackRef trk(tracks, trkIndx);
0210 if ((*trk).charge() * protonCharge < 0 || trkStatus[trkIndx] != TrkStatus::NOTUSED)
0211 continue;
0212
0213 math::XYZTLorentzVector p4pi(
0214 trk->px(), trk->py(), trk->pz(), sqrt(trk->momentum().mag2() + 0.01947995518));
0215 if ((p4Lambda + p4pi).M() < xiMassCut_)
0216 trkStatus[trkIndx] = TrkStatus::VTX;
0217 }
0218 }
0219 }
0220 std::vector<int> mapping(tracks->size(), -1);
0221 int lostTrkIndx = 0;
0222 for (unsigned int trkIndx = 0; trkIndx < tracks->size(); trkIndx++) {
0223 reco::TrackRef trk(tracks, trkIndx);
0224 if (trkStatus[trkIndx] == TrkStatus::VTX || (trkStatus[trkIndx] == TrkStatus::NOTUSED && passTrkCuts(*trk))) {
0225 outPtrTrks->emplace_back(*trk);
0226
0227 std::pair<int, pat::PackedCandidate::PVAssociationQuality> pvAsso = associateTrkToVtx(*pvOrigs, trk);
0228 const reco::VertexRef& pvOrigRef = reco::VertexRef(pvOrigs, pvAsso.first);
0229 if (pvOrigRef.isNonnull()) {
0230 pv = reco::VertexRef(pvs, pvOrigRef.key());
0231 } else if (!pvs->empty()) {
0232 pv = reco::VertexRef(pvs, 0);
0233 }
0234 addPackedCandidate(*outPtrTrksAsCands, trk, pv, pvRefProd, trkStatus[trkIndx], pvAsso.second, muons);
0235
0236
0237
0238 mapping[trkIndx] = lostTrkIndx;
0239 lostTrkIndx++;
0240 } else if ((trkStatus[trkIndx] == TrkStatus::PFELECTRON || trkStatus[trkIndx] == TrkStatus::PFPOSITRON) &&
0241 passTrkCuts(*trk)) {
0242
0243 std::pair<int, pat::PackedCandidate::PVAssociationQuality> pvAsso = associateTrkToVtx(*pvOrigs, trk);
0244 const reco::VertexRef& pvOrigRef = reco::VertexRef(pvOrigs, pvAsso.first);
0245 if (pvOrigRef.isNonnull()) {
0246 pv = reco::VertexRef(pvs, pvOrigRef.key());
0247 } else if (!pvs->empty()) {
0248 pv = reco::VertexRef(pvs, 0);
0249 }
0250 addPackedCandidate(*outPtrEleTrksAsCands, trk, pv, pvRefProd, trkStatus[trkIndx], pvAsso.second, muons);
0251 }
0252 }
0253
0254 iEvent.put(std::move(outPtrTrks));
0255 iEvent.put(std::move(outPtrEleTrksAsCands), "eleTracks");
0256 edm::OrphanHandle<pat::PackedCandidateCollection> oh = iEvent.put(std::move(outPtrTrksAsCands));
0257 auto tk2pc = std::make_unique<edm::Association<pat::PackedCandidateCollection>>(oh);
0258 edm::Association<pat::PackedCandidateCollection>::Filler tk2pcFiller(*tk2pc);
0259 tk2pcFiller.insert(tracks, mapping.begin(), mapping.end());
0260 tk2pcFiller.fill();
0261 iEvent.put(std::move(tk2pc));
0262 }
0263
0264 bool pat::PATLostTracks::passTrkCuts(const reco::Track& tr) const {
0265 const bool passTrkHits = tr.pt() > minPt_ && tr.numberOfValidHits() >= minHits_ &&
0266 tr.hitPattern().numberOfValidPixelHits() >= minPixelHits_;
0267 const bool passTrkQual = passesQuality(tr, qualsToAutoAccept_);
0268
0269 return passTrkHits || passTrkQual || passThroughCut_(tr);
0270 }
0271
0272 void pat::PATLostTracks::addPackedCandidate(std::vector<pat::PackedCandidate>& cands,
0273 const reco::TrackRef& trk,
0274 const reco::VertexRef& pvSlimmed,
0275 const reco::VertexRefProd& pvSlimmedColl,
0276 const pat::PATLostTracks::TrkStatus& trkStatus,
0277 const pat::PackedCandidate::PVAssociationQuality& pvAssocQuality,
0278 edm::Handle<reco::MuonCollection> muons) const {
0279 const float mass = 0.13957018;
0280
0281 int id = 211 * trk->charge();
0282 if (trkStatus == TrkStatus::PFELECTRON)
0283 id = 11;
0284 else if (trkStatus == TrkStatus::PFPOSITRON)
0285 id = -11;
0286
0287
0288 const reco::Muon* muon(nullptr);
0289 for (auto& mu : *muons) {
0290 if (reco::TrackRef(mu.innerTrack()) == trk) {
0291 id = -13 * trk->charge();
0292 muon = μ
0293 break;
0294 }
0295 }
0296
0297 pat::PackedCandidate::LostInnerHits lostHits = pat::PackedCandidate::noLostInnerHits;
0298 int nlost = trk->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS);
0299 if (nlost == 0) {
0300 if (trk->hitPattern().hasValidHitInPixelLayer(PixelSubdetector::SubDetector::PixelBarrel, 1)) {
0301 lostHits = pat::PackedCandidate::validHitInFirstPixelBarrelLayer;
0302 }
0303 } else {
0304 lostHits = (nlost == 1 ? pat::PackedCandidate::oneLostInnerHit : pat::PackedCandidate::moreLostInnerHits);
0305 }
0306
0307 reco::Candidate::PolarLorentzVector p4(trk->pt(), trk->eta(), trk->phi(), mass);
0308 cands.emplace_back(
0309 pat::PackedCandidate(p4, trk->vertex(), trk->pt(), trk->eta(), trk->phi(), id, pvSlimmedColl, pvSlimmed.key()));
0310
0311 cands.back().setTrackHighPurity(trk->quality(reco::TrackBase::highPurity));
0312 cands.back().setCovarianceVersion(covarianceVersion_);
0313 cands.back().setLostInnerHits(lostHits);
0314 if (trk->pt() > minPtToStoreProps_ || trkStatus == TrkStatus::VTX) {
0315 cands.back().setTrkAlgo(static_cast<uint8_t>(trk->algo()), static_cast<uint8_t>(trk->originalAlgo()));
0316 cands.back().setFirstHit(trk->hitPattern().getHitPattern(reco::HitPattern::TRACK_HITS, 0));
0317 if (useLegacySetup_ || std::abs(id) == 11 || trkStatus == TrkStatus::VTX) {
0318 cands.back().setTrackProperties(*trk, covariancePackingSchemas_[4], covarianceVersion_);
0319 } else {
0320 if (trk->hitPattern().numberOfValidPixelHits() > 0) {
0321 cands.back().setTrackProperties(
0322 *trk, covariancePackingSchemas_[0], covarianceVersion_);
0323 } else {
0324 cands.back().setTrackProperties(
0325 *trk, covariancePackingSchemas_[1], covarianceVersion_);
0326 }
0327 }
0328 } else if (!useLegacySetup_ && trk->pt() > minPtToStoreLowQualityProps_) {
0329 if (trk->hitPattern().numberOfValidPixelHits() > 0) {
0330 cands.back().setTrackProperties(
0331 *trk, covariancePackingSchemas_[2], covarianceVersion_);
0332 } else {
0333 cands.back().setTrackProperties(
0334 *trk, covariancePackingSchemas_[3], covarianceVersion_);
0335 }
0336 }
0337 cands.back().setAssociationQuality(pvAssocQuality);
0338
0339 if (muon)
0340 cands.back().setMuonID(muon->isStandAloneMuon(), muon->isGlobalMuon());
0341 }
0342
0343 std::pair<int, pat::PackedCandidate::PVAssociationQuality> pat::PATLostTracks::associateTrkToVtx(
0344 const reco::VertexCollection& vertices, const reco::TrackRef& trk) const {
0345
0346
0347
0348 if (useLegacySetup_) {
0349 float w = vertices[0].trackWeight(trk);
0350 if (w > 0.5) {
0351 return std::pair<int, pat::PackedCandidate::PVAssociationQuality>(0, pat::PackedCandidate::UsedInFitTight);
0352 } else if (w > 0.) {
0353 return std::pair<int, pat::PackedCandidate::PVAssociationQuality>(0,
0354 pat::PackedCandidate::NotReconstructedPrimary);
0355 } else {
0356 return std::pair<int, pat::PackedCandidate::PVAssociationQuality>(-1,
0357 pat::PackedCandidate::NotReconstructedPrimary);
0358 }
0359 }
0360
0361
0362
0363
0364 int iVtxMaxWeight = -1;
0365 int iVtxMinDzDist = -1;
0366 size_t idx = 0;
0367 float maxWeight = 0;
0368 double minDz = std::numeric_limits<double>::max();
0369 double minDzSig = std::numeric_limits<double>::max();
0370 for (auto const& vtx : vertices) {
0371 float w = vtx.trackWeight(trk);
0372 double dz = std::abs(trk->dz(vtx.position()));
0373 double dzSig = dz / trk->dzError();
0374 if (w > maxWeight) {
0375 maxWeight = w;
0376 iVtxMaxWeight = idx;
0377 }
0378 if (dzSig < minDzSig) {
0379 minDzSig = dzSig;
0380 minDz = dz;
0381 iVtxMinDzDist = idx;
0382 }
0383 idx++;
0384 }
0385
0386 if (iVtxMaxWeight >= 0) {
0387 if (maxWeight > 0.5) {
0388 return std::pair<int, pat::PackedCandidate::PVAssociationQuality>(iVtxMaxWeight,
0389 pat::PackedCandidate::UsedInFitTight);
0390 } else {
0391 return std::pair<int, pat::PackedCandidate::PVAssociationQuality>(iVtxMaxWeight,
0392 pat::PackedCandidate::UsedInFitLoose);
0393 }
0394 }
0395
0396 if (minDz < maxDzForPrimaryAssignment_) {
0397 const double add_cov = vertices[iVtxMinDzDist].covariance(2, 2);
0398 const double dzErr = sqrt(trk->dzError() * trk->dzError() + add_cov);
0399 if (minDz / dzErr < maxDzSigForPrimaryAssignment_ && trk->dzError() < maxDzErrorForPrimaryAssignment_) {
0400 return std::pair<int, pat::PackedCandidate::PVAssociationQuality>(iVtxMinDzDist,
0401 pat::PackedCandidate::CompatibilityDz);
0402 }
0403 }
0404
0405
0406 if (!vertices.empty() && std::abs(trk->dxy(vertices[0].position())) < maxDxyForNotReconstructedPrimary_ &&
0407 std::abs(trk->dxy(vertices[0].position()) / trk->dxyError()) < maxDxySigForNotReconstructedPrimary_)
0408 return std::pair<int, pat::PackedCandidate::PVAssociationQuality>(iVtxMinDzDist,
0409 pat::PackedCandidate::NotReconstructedPrimary);
0410
0411 return std::pair<int, pat::PackedCandidate::PVAssociationQuality>(iVtxMinDzDist, pat::PackedCandidate::OtherDeltaZ);
0412 }
0413
0414 using pat::PATLostTracks;
0415 #include "FWCore/Framework/interface/MakerMacros.h"
0416 DEFINE_FWK_MODULE(PATLostTracks);