FinalTrackSelectors/plugins/SingleLongTrackProducer.cc

0001 // user includes
0002 #include "DataFormats/Math/interface/deltaR.h"
0003 #include "DataFormats/MuonReco/interface/Muon.h"
0004 #include "DataFormats/TrackReco/interface/Track.h"
0005 #include "DataFormats/VertexReco/interface/Vertex.h"
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "FWCore/Framework/interface/stream/EDProducer.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0012
0013 // ROOT includes
0014 #include "TLorentzVector.h"
0015
0016 class SingleLongTrackProducer : public edm::stream::EDProducer<> {
0017 public:
0018   explicit SingleLongTrackProducer(const edm::ParameterSet &);
0019   ~SingleLongTrackProducer() override = default;
0020
0021   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0022
0023 private:
0024   void produce(edm::Event &, const edm::EventSetup &) override;
0025
0026   const edm::EDGetTokenT<std::vector<reco::Track>> tracksToken;
0027   const edm::EDGetTokenT<std::vector<reco::Muon>> muonsToken;
0028   const edm::EDGetTokenT<reco::VertexCollection> PrimVtxToken;
0029
0030   const int minNumberOfLayers;
0031   const double matchInDr;
0032   const bool onlyValidHits;
0033   const bool debug;
0034   const double minPt;
0035   const double maxEta;
0036   const double maxDxy;
0037   const double maxDz;
0038 };
0039
0040 SingleLongTrackProducer::SingleLongTrackProducer(const edm::ParameterSet &iConfig)
0041     : tracksToken{consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("allTracks"))},
0042       muonsToken{consumes<std::vector<reco::Muon>>(iConfig.getParameter<edm::InputTag>("matchMuons"))},
0043       PrimVtxToken{consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("PrimaryVertex"))},
0044       minNumberOfLayers{iConfig.getParameter<int>("minNumberOfLayers")},
0045       matchInDr{iConfig.getParameter<double>("requiredDr")},
0046       onlyValidHits{iConfig.getParameter<bool>("onlyValidHits")},
0047       debug{iConfig.getParameter<bool>("debug")},
0048       minPt{iConfig.getParameter<double>("minPt")},
0049       maxEta{iConfig.getParameter<double>("maxEta")},
0050       maxDxy{iConfig.getParameter<double>("maxDxy")},
0051       maxDz{iConfig.getParameter<double>("maxDz")} {
0052   produces<reco::TrackCollection>("").setBranchAlias("");
0053 }
0054
0055 void SingleLongTrackProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0056   using namespace edm;
0057
0058   // register output collection:
0059   std::unique_ptr<reco::TrackCollection> goodTracks(new reco::TrackCollection);
0060
0061   // register input collections:
0062   const auto &tracks = iEvent.getHandle(tracksToken);
0063   if (!tracks.isValid()) {
0064     edm::LogError("SingleLongTrackProducer")
0065         << "Input track collection is not valid.\n Returning empty output track collection.";
0066     iEvent.put(std::move(goodTracks), "");
0067     return;
0068   }
0069
0070   const auto &muons = iEvent.getHandle(muonsToken);
0071   if (!muons.isValid() && matchInDr > 0.) {
0072     edm::LogError("SingleLongTrackProducer")
0073         << "Input muon collection is not valid.\n Returning empty output track collection.";
0074     iEvent.put(std::move(goodTracks), "");
0075     return;
0076   }
0077
0078   const auto &vertices = iEvent.getHandle(PrimVtxToken);
0079   if (!vertices.isValid()) {
0080     edm::LogError("SingleLongTrackProducer")
0081         << "Input vertex collection is not valid.\n Returning empty output track collection.";
0082     iEvent.put(std::move(goodTracks), "");
0083     return;
0084   }
0085
0086   // Preselection of long quality tracks
0087   const reco::Vertex &vtx = vertices->at(0);
0088   std::vector<reco::Track> selTracks;
0089   reco::Track bestTrack;
0090   unsigned int tMuon = 0;
0091   double fitProb = std::numeric_limits<double>::max();
0092
0093   for (const auto &track : *tracks) {
0094     const reco::HitPattern &hitpattern = track.hitPattern();
0095     double dR2min = std::numeric_limits<double>::max();
0096     double chiNdof = track.normalizedChi2();
0097     double dxy = std::abs(track.dxy(vtx.position()));
0098     double dz = std::abs(track.dz(vtx.position()));
0099
0100     if (track.pt() < minPt || std::abs(track.eta()) > maxEta ||
0101         hitpattern.trackerLayersWithMeasurement() < minNumberOfLayers) {
0102       continue;
0103     }
0104
0105     if (matchInDr > 0.0) {
0106       for (const auto &muon : *muons) {
0107         if (muon.isTrackerMuon()) {
0108           tMuon++;
0109           double dr2 = reco::deltaR2(track, *(muon.innerTrack()));
0110           dR2min = std::min(dR2min, dr2);
0111         }
0112       }
0113       if (dR2min >= matchInDr * matchInDr)
0114         continue;
0115     }
0116
0117     if (dxy < maxDxy && dz < maxDz && track.validFraction() >= 1.0 && chiNdof < fitProb) {
0118       fitProb = chiNdof;
0119       bestTrack = track;
0120     }
0121
0122     if (debug) {
0123       edm::LogPrint("SingleLongTrackProducer") << "deltaR2 (general) track to matched Track: " << dR2min
0124                                                << " chi2Ndof: " << chiNdof << " best Track chi2Ndof: " << fitProb;
0125     }
0126   }
0127
0128   // Only push if bestTrack was set
0129   if (bestTrack.recHitsOk()) {
0130     selTracks.push_back(bestTrack);
0131   }
0132
0133   if (debug)
0134     edm::LogPrint("SingleLongTrackProducer")
0135         << " number of Tracker Muons: " << tMuon << ", thereof " << selTracks.size() << " tracks passed preselection.";
0136
0137   // check hits validity in preselected tracks
0138   bool hitIsNotValid{false};
0139
0140   for (const auto &track : selTracks) {
0141     // Check validity of track extra and rechits
0142     if (track.recHitsOk()) {
0143       reco::HitPattern hitpattern = track.hitPattern();
0144       auto hb = track.recHitsBegin();
0145
0146       // Checking if rechits are valid
0147       for (unsigned int h = 0; h < track.recHitsSize(); h++) {
0148         auto recHit = *(hb + h);
0149         auto const &hit = *recHit;
0150         if (onlyValidHits && !hit.isValid()) {
0151           hitIsNotValid = true;
0152           continue;
0153         }
0154       }
0155
0156       if (hitIsNotValid == true)
0157         break;  // (Un)Comment this line to (not) allow for events with not valid hits
0158
0159       // Checking if hitpattern hits are valid
0160       for (unsigned int h = 0; h < track.recHitsSize(); h++) {
0161         uint32_t pHit = hitpattern.getHitPattern(reco::HitPattern::TRACK_HITS, h);
0162         if (onlyValidHits && !(hitpattern.validHitFilter(pHit))) {
0163           if (debug)
0164             edm::LogPrint("SingleLongTrackProducer") << "hit not valid: " << h;
0165           continue;
0166         }
0167       }
0168       goodTracks->push_back(track);
0169     } else
0170       break;
0171   }
0172
0173   if (debug) {
0174     auto const &moduleType = moduleDescription().moduleName();
0175     auto const &moduleLabel = moduleDescription().moduleLabel();
0176     edm::LogPrint("SingleLongTrackProducer") << "[" << moduleType << "] (" << moduleLabel << ") "
0177                                              << " output track size: " << goodTracks.get()->size();
0178   }
0179
0180   // save track collection in event:
0181   iEvent.put(std::move(goodTracks), "");
0182 }
0183
0184 void SingleLongTrackProducer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0185   edm::ParameterSetDescription desc;
0186   desc.add<edm::InputTag>("allTracks", edm::InputTag("generalTracks"))->setComment("input track collection");
0187   desc.add<edm::InputTag>("matchMuons", edm::InputTag("earlyMuons"))->setComment("input muon collection for matching");
0188   desc.add<edm::InputTag>("PrimaryVertex", edm::InputTag("offlinePrimaryVertices"))
0189       ->setComment("input primary vertex collection");
0190   desc.add<int>("minNumberOfLayers", 10)->setComment("minimum number of layers");
0191   desc.add<double>("requiredDr", 0.01)->setComment("matching muons deltaR. If negative do not match");
0192   desc.add<bool>("onlyValidHits", true)->setComment("use only valid hits");
0193   desc.add<bool>("debug", false)->setComment("verbose?");
0194   desc.add<double>("minPt", 15.0)->setComment("minimum pT");
0195   desc.add<double>("maxEta", 2.2)->setComment("maximum pseudorapidity (absolute value)");
0196   desc.add<double>("maxDxy", 0.02)->setComment("maximum transverse impact parameter");
0197   desc.add<double>("maxDz", 0.5)->setComment("maximum longitudinal impact parameter");
0198   descriptions.addWithDefaultLabel(desc);
0199 }
0200
0201 #include "FWCore/Framework/interface/MakerMacros.h"
0202 DEFINE_FWK_MODULE(SingleLongTrackProducer);