RecoAlgos/src/TrackWithVertexSelector.cc

0001 #include "CommonTools/RecoAlgos/interface/TrackWithVertexSelector.h"
0002 #include "FWCore/Utilities/interface/isFinite.h"
0003 //
0004 // constructors and destructor
0005 //
0006
0007 namespace {
0008   constexpr float fakeBeamSpotTimeWidth = 0.175f;
0009 }
0010
0011 TrackWithVertexSelector::TrackWithVertexSelector(const edm::ParameterSet &iConfig, edm::ConsumesCollector &iC)
0012     : numberOfValidHits_(iConfig.getParameter<uint32_t>("numberOfValidHits")),
0013       numberOfValidPixelHits_(iConfig.getParameter<uint32_t>("numberOfValidPixelHits")),
0014       numberOfLostHits_(iConfig.getParameter<uint32_t>("numberOfLostHits")),
0015       normalizedChi2_(iConfig.getParameter<double>("normalizedChi2")),
0016       ptMin_(iConfig.getParameter<double>("ptMin")),
0017       ptMax_(iConfig.getParameter<double>("ptMax")),
0018       etaMin_(iConfig.getParameter<double>("etaMin")),
0019       etaMax_(iConfig.getParameter<double>("etaMax")),
0020       dzMax_(iConfig.getParameter<double>("dzMax")),
0021       d0Max_(iConfig.getParameter<double>("d0Max")),
0022       ptErrorCut_(iConfig.getParameter<double>("ptErrorCut")),
0023       quality_(iConfig.getParameter<std::string>("quality")),
0024       nVertices_(iConfig.getParameter<bool>("useVtx") ? iConfig.getParameter<uint32_t>("nVertices") : 0),
0025       vertexToken_(iC.consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertexTag"))),
0026       timesToken_(iC.consumes<edm::ValueMap<float> >(iConfig.getParameter<edm::InputTag>("timesTag"))),
0027       timeResosToken_(iC.consumes<edm::ValueMap<float> >(iConfig.getParameter<edm::InputTag>("timeResosTag"))),
0028       vtxFallback_(iConfig.getParameter<bool>("vtxFallback")),
0029       zetaVtx_(iConfig.getParameter<double>("zetaVtx")),
0030       rhoVtx_(iConfig.getParameter<double>("rhoVtx")),
0031       nSigmaDtVertex_(iConfig.getParameter<double>("nSigmaDtVertex")) {}
0032
0033 TrackWithVertexSelector::~TrackWithVertexSelector() {}
0034
0035 void TrackWithVertexSelector::init(const edm::Event &event) {
0036   edm::Handle<reco::VertexCollection> hVtx;
0037   event.getByToken(vertexToken_, hVtx);
0038   vcoll_ = nVertices_ > 0 ? hVtx.product() : nullptr;
0039
0040   edm::Handle<edm::ValueMap<float> > hTimes;
0041   event.getByToken(timesToken_, hTimes);
0042   timescoll_ = hTimes.isValid() ? hTimes.product() : nullptr;
0043
0044   edm::Handle<edm::ValueMap<float> > hTimeResos;
0045   event.getByToken(timeResosToken_, hTimeResos);
0046   timeresoscoll_ = hTimeResos.isValid() ? hTimeResos.product() : nullptr;
0047 }
0048
0049 bool TrackWithVertexSelector::testTrack(const reco::Track &t) const {
0050   using std::abs;
0051   if ((t.numberOfValidHits() >= numberOfValidHits_) &&
0052       (static_cast<unsigned int>(t.hitPattern().numberOfValidPixelHits()) >= numberOfValidPixelHits_) &&
0053       (t.numberOfLostHits() <= numberOfLostHits_) && (t.normalizedChi2() <= normalizedChi2_) &&
0054       (t.ptError() / t.pt() * std::max(1., t.normalizedChi2()) <= ptErrorCut_) &&
0055       (t.quality(t.qualityByName(quality_))) && (t.pt() >= ptMin_) && (t.pt() <= ptMax_) && (abs(t.eta()) <= etaMax_) &&
0056       (abs(t.eta()) >= etaMin_) && (abs(t.dz()) <= dzMax_) && (abs(t.d0()) <= d0Max_)) {
0057     return true;
0058   }
0059   return false;
0060 }
0061
0062 bool TrackWithVertexSelector::testTrack(const reco::TrackRef &tref) const { return testTrack(*tref); }
0063
0064 bool TrackWithVertexSelector::testVertices(const reco::Track &t, const reco::VertexCollection &vtxs) const {
0065   bool ok = false;
0066   if (!vtxs.empty()) {
0067     unsigned int tested = 1;
0068     for (reco::VertexCollection::const_iterator it = vtxs.begin(), ed = vtxs.end(); it != ed; ++it) {
0069       if ((std::abs(t.dxy(it->position())) < rhoVtx_) && (std::abs(t.dz(it->position())) < zetaVtx_)) {
0070         ok = true;
0071         break;
0072       }
0073       if (tested++ >= nVertices_)
0074         break;
0075     }
0076   } else if (vtxFallback_) {
0077     return ((std::abs(t.vertex().z()) < 15.9) && (t.vertex().Rho() < 0.2));
0078   }
0079   return ok;
0080 }
0081
0082 bool TrackWithVertexSelector::testVertices(const reco::TrackRef &tref, const reco::VertexCollection &vtxs) const {
0083   const auto &t = *tref;
0084   const bool timeAvailable = timescoll_ != nullptr && timeresoscoll_ != nullptr;
0085   bool ok = false;
0086   if (!vtxs.empty()) {
0087     unsigned int tested = 1;
0088     for (reco::VertexCollection::const_iterator it = vtxs.begin(), ed = vtxs.end(); it != ed; ++it) {
0089       const bool useTime = timeAvailable && it->t() != 0.;
0090       float time = useTime ? (*timescoll_)[tref] : -1.f;
0091       float timeReso = useTime ? (*timeresoscoll_)[tref] : -1.f;
0092       timeReso = (timeReso > 1e-6 ? timeReso : fakeBeamSpotTimeWidth);
0093
0094       if (edm::isNotFinite(time)) {
0095         time = 0.0;
0096         timeReso = 2.0 * fakeBeamSpotTimeWidth;
0097       }
0098
0099       const double vtxSigmaT2 = it->tError() * it->tError();
0100       const double vtxTrackErr = std::sqrt(vtxSigmaT2 + timeReso * timeReso);
0101
0102       if ((std::abs(t.dxy(it->position())) < rhoVtx_) && (std::abs(t.dz(it->position())) < zetaVtx_) &&
0103           (!useTime || (std::abs(time - it->t()) / vtxTrackErr < nSigmaDtVertex_))) {
0104         ok = true;
0105         break;
0106       }
0107       if (tested++ >= nVertices_)
0108         break;
0109     }
0110   } else if (vtxFallback_) {
0111     return ((std::abs(t.vertex().z()) < 15.9) && (t.vertex().Rho() < 0.2));
0112   }
0113   return ok;
0114 }
0115
0116 bool TrackWithVertexSelector::operator()(const reco::Track &t) const {
0117   if (!testTrack(t))
0118     return false;
0119   if (nVertices_ == 0)
0120     return true;
0121   return testVertices(t, *vcoll_);
0122 }
0123
0124 bool TrackWithVertexSelector::operator()(const reco::TrackRef &tref) const {
0125   if (!testTrack(tref))
0126     return false;
0127   if (nVertices_ == 0)
0128     return true;
0129   return testVertices(tref, *vcoll_);
0130 }