FinalTrackSelectors/plugins/TrackCutClassifier.cc

0001 #include "RecoTracker/FinalTrackSelectors/interface/TrackMVAClassifier.h"
0002
0003 #include "DataFormats/TrackReco/interface/Track.h"
0004
0005 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0006
0007 #include "RecoTracker/FinalTrackSelectors/interface/getBestVertex.h"
0008
0009 #include <cassert>
0010
0011 #include "powN.h"
0012
0013 namespace {
0014
0015   void fillArrayF(float* x, const edm::ParameterSet& cfg, const char* name) {
0016     auto v = cfg.getParameter<std::vector<double>>(name);
0017     assert(v.size() == 3);
0018     std::copy(std::begin(v), std::end(v), x);
0019   }
0020
0021   void fillArrayI(int* x, const edm::ParameterSet& cfg, const char* name) {
0022     auto v = cfg.getParameter<std::vector<int>>(name);
0023     assert(v.size() == 3);
0024     std::copy(std::begin(v), std::end(v), x);
0025   }
0026
0027   // fake mva value to return for loose,tight,hp
0028   constexpr float mvaVal[3] = {-.5, .5, 1.};
0029
0030   template <typename T, typename Comp>
0031   inline float cut(T val, const T* cuts, Comp comp) {
0032     for (int i = 2; i >= 0; --i) {
0033       if (comp(val, cuts[i]))
0034         return mvaVal[i];
0035     }
0036     return -1.f;
0037   }
0038
0039   inline float chi2n(reco::Track const& tk) { return tk.normalizedChi2(); }
0040
0041   inline float relPtErr(reco::Track const& tk) {
0042     return (tk.pt() != 0. ? float(tk.ptError()) / float(tk.pt()) : 9999999.);
0043   }
0044
0045   inline int lostLayers(reco::Track const& tk) {
0046     return tk.hitPattern().trackerLayersWithoutMeasurement(reco::HitPattern::TRACK_HITS);
0047   }
0048
0049   inline int n3DLayers(reco::Track const& tk, bool isHLT) {
0050     uint32_t nlayers3D = tk.hitPattern().pixelLayersWithMeasurement();
0051     if (!isHLT)
0052       nlayers3D += tk.hitPattern().numberOfValidStripLayersWithMonoAndStereo();
0053     else {
0054       size_t count3D = 0;
0055       for (auto it = tk.recHitsBegin(), et = tk.recHitsEnd(); it != et; ++it) {
0056         const TrackingRecHit* hit = (*it);
0057         if (!trackerHitRTTI::isFromDetOrFast(*hit))
0058           continue;
0059
0060         if (hit->dimension() == 2) {
0061           auto const& thit = static_cast<BaseTrackerRecHit const&>(*hit);
0062           if (thit.isMatched())
0063             count3D++;
0064         }
0065       }
0066       nlayers3D += count3D;
0067     }
0068     return nlayers3D;
0069   }
0070
0071   inline int nHits(reco::Track const& tk) { return tk.numberOfValidHits(); }
0072
0073   inline int nPixelHits(reco::Track const& tk) { return tk.hitPattern().numberOfValidPixelHits(); }
0074
0075   inline float dz(reco::Track const& trk, Point const& bestVertex) { return std::abs(trk.dz(bestVertex)); }
0076   inline float dr(reco::Track const& trk, Point const& bestVertex) { return std::abs(trk.dxy(bestVertex)); }
0077
0078   inline void dzCut_par1(reco::Track const& trk, int& nLayers, const float* par, const int* exp, float dzCut[]) {
0079     float dzE = trk.dzError();
0080     for (int i = 2; i >= 0; --i) {
0081       dzCut[i] = powN(par[i] * nLayers, exp[i]) * dzE;
0082     }
0083   }
0084   inline void drCut_par1(reco::Track const& trk, int& nLayers, const float* par, const int* exp, float drCut[]) {
0085     float drE = trk.d0Error();
0086     for (int i = 2; i >= 0; --i) {
0087       drCut[i] = powN(par[i] * nLayers, exp[i]) * drE;
0088     }
0089   }
0090
0091   inline void dzCut_par2(reco::Track const& trk,
0092                          int& nLayers,
0093                          const float* par,
0094                          const int* exp,
0095                          const float* d0err,
0096                          const float* d0err_par,
0097                          float dzCut[]) {
0098     float pt = float(trk.pt());
0099     float p = float(trk.p());
0100
0101     for (int i = 2; i >= 0; --i) {
0102       // parametrized d0 resolution for the track pt
0103       float nomd0E = sqrt(d0err[i] * d0err[i] + (d0err_par[i] / pt) * (d0err_par[i] / pt));
0104       // parametrized z0 resolution for the track pt and eta
0105       float nomdzE = nomd0E * (p / pt);  // cosh(eta):=abs(p)/pt
0106
0107       dzCut[i] = powN(par[i] * nLayers, exp[i]) * nomdzE;
0108     }
0109   }
0110   inline void drCut_par2(reco::Track const& trk,
0111                          int& nLayers,
0112                          const float* par,
0113                          const int* exp,
0114                          const float* d0err,
0115                          const float* d0err_par,
0116                          float drCut[]) {
0117     float pt = trk.pt();
0118
0119     for (int i = 2; i >= 0; --i) {
0120       // parametrized d0 resolution for the track pt
0121       float nomd0E = sqrt(d0err[i] * d0err[i] + (d0err_par[i] / pt) * (d0err_par[i] / pt));
0122
0123       drCut[i] = powN(par[i] * nLayers, exp[i]) * nomd0E;
0124     }
0125   }
0126
0127   inline void dzCut_wPVerror_par(reco::Track const& trk,
0128                                  int& nLayers,
0129                                  const float* par,
0130                                  const int* exp,
0131                                  Point const& bestVertexError,
0132                                  float dzCut[]) {
0133     float dzE = trk.dzError();
0134     float zPVerr = bestVertexError.z();
0135
0136     float dzErrPV = std::sqrt(dzE * dzE + zPVerr * zPVerr);
0137     for (int i = 2; i >= 0; --i) {
0138       dzCut[i] = par[i] * dzErrPV;
0139       if (exp[i] != 0)
0140         dzCut[i] *= pow(nLayers, exp[i]);
0141     }
0142   }
0143
0144   inline void drCut_wPVerror_par(reco::Track const& trk,
0145                                  int& nLayers,
0146                                  const float* par,
0147                                  const int* exp,
0148                                  Point const& bestVertexError,
0149                                  float drCut[]) {
0150     float drE = trk.d0Error();
0151     // shouldn't it be bestVertex.xError()*bestVertex.xError()+bestVertex.yError()*bestVertex.yError() ?!?!?
0152     float rPVerr = sqrt(bestVertexError.x() * bestVertexError.y());
0153
0154     float drErrPV = std::sqrt(drE * drE + rPVerr * rPVerr);
0155     for (int i = 2; i >= 0; --i) {
0156       drCut[i] = par[i] * drErrPV;
0157       if (exp[i] != 0)
0158         drCut[i] *= pow(nLayers, exp[i]);
0159     }
0160   }
0161
0162   struct Cuts {
0163     Cuts(const edm::ParameterSet& cfg, edm::ConsumesCollector) {
0164       isHLT = cfg.getParameter<bool>("isHLT");
0165       fillArrayF(minNdof, cfg, "minNdof");
0166       fillArrayF(maxChi2, cfg, "maxChi2");
0167       fillArrayF(maxChi2n, cfg, "maxChi2n");
0168       fillArrayI(minHits4pass, cfg, "minHits4pass");
0169       fillArrayI(minHits, cfg, "minHits");
0170       fillArrayI(minPixelHits, cfg, "minPixelHits");
0171       fillArrayI(min3DLayers, cfg, "min3DLayers");
0172       fillArrayI(minLayers, cfg, "minLayers");
0173       fillArrayI(maxLostLayers, cfg, "maxLostLayers");
0174       fillArrayF(maxRelPtErr, cfg, "maxRelPtErr");
0175       minNVtxTrk = cfg.getParameter<int>("minNVtxTrk");
0176       fillArrayF(maxDz, cfg, "maxDz");
0177       fillArrayF(maxDzWrtBS, cfg, "maxDzWrtBS");
0178       fillArrayF(maxDr, cfg, "maxDr");
0179       edm::ParameterSet dz_par = cfg.getParameter<edm::ParameterSet>("dz_par");
0180       fillArrayI(dz_exp, dz_par, "dz_exp");
0181       fillArrayF(dz_par1, dz_par, "dz_par1");
0182       fillArrayF(dz_par2, dz_par, "dz_par2");
0183       fillArrayF(dzWPVerr_par, dz_par, "dzWPVerr_par");
0184       edm::ParameterSet dr_par = cfg.getParameter<edm::ParameterSet>("dr_par");
0185       fillArrayI(dr_exp, dr_par, "dr_exp");
0186       fillArrayF(dr_par1, dr_par, "dr_par1");
0187       fillArrayF(dr_par2, dr_par, "dr_par2");
0188       fillArrayF(d0err, dr_par, "d0err");
0189       fillArrayF(d0err_par, dr_par, "d0err_par");
0190       fillArrayF(drWPVerr_par, dr_par, "drWPVerr_par");
0191     }
0192
0193     void beginStream() {}
0194     void initEvent(const edm::EventSetup&) {}
0195
0196     float operator()(reco::Track const& trk,
0197                      reco::BeamSpot const& beamSpot,
0198                      reco::VertexCollection const& vertices) const {
0199       float ret = 1.f;
0200       // minimum number of hits for by-passing the other checks
0201       if (minHits4pass[0] < std::numeric_limits<int>::max()) {
0202         ret = std::min(ret, cut(nHits(trk), minHits4pass, std::greater_equal<int>()));
0203         if (ret == 1.f)
0204           return ret;
0205       }
0206
0207       if (maxRelPtErr[2] < std::numeric_limits<float>::max()) {
0208         ret = std::min(ret, cut(relPtErr(trk), maxRelPtErr, std::less_equal<float>()));
0209         if (ret == -1.f)
0210           return ret;
0211       }
0212
0213       ret = std::min(ret, cut(float(trk.ndof()), minNdof, std::greater_equal<float>()));
0214       if (ret == -1.f)
0215         return ret;
0216
0217       auto nLayers = trk.hitPattern().trackerLayersWithMeasurement();
0218       ret = std::min(ret, cut(nLayers, minLayers, std::greater_equal<int>()));
0219       if (ret == -1.f)
0220         return ret;
0221
0222       ret = std::min(ret, cut(chi2n(trk) / float(nLayers), maxChi2n, std::less_equal<float>()));
0223       if (ret == -1.f)
0224         return ret;
0225
0226       ret = std::min(ret, cut(chi2n(trk), maxChi2, std::less_equal<float>()));
0227       if (ret == -1.f)
0228         return ret;
0229
0230       ret = std::min(ret, cut(n3DLayers(trk, isHLT), min3DLayers, std::greater_equal<int>()));
0231       if (ret == -1.f)
0232         return ret;
0233
0234       ret = std::min(ret, cut(nHits(trk), minHits, std::greater_equal<int>()));
0235       if (ret == -1.f)
0236         return ret;
0237
0238       ret = std::min(ret, cut(nPixelHits(trk), minPixelHits, std::greater_equal<int>()));
0239       if (ret == -1.f)
0240         return ret;
0241
0242       ret = std::min(ret, cut(lostLayers(trk), maxLostLayers, std::less_equal<int>()));
0243       if (ret == -1.f)
0244         return ret;
0245
0246       // original dz and dr cut
0247       if (maxDz[2] < std::numeric_limits<float>::max() || maxDr[2] < std::numeric_limits<float>::max()) {
0248         // if not primaryVertices are reconstructed, check compatibility w.r.t. beam spot
0249         // min number of tracks [2 (=default) for offline, 3 for HLT]
0250         Point bestVertex = getBestVertex(trk, vertices, minNVtxTrk);
0251         float maxDzcut[3];
0252         std::copy(std::begin(maxDz), std::end(maxDz), std::begin(maxDzcut));
0253         if (bestVertex.z() < -99998.) {
0254           bestVertex = beamSpot.position();
0255           std::copy(std::begin(maxDzWrtBS), std::end(maxDzWrtBS), std::begin(maxDzcut));
0256         }
0257         ret = std::min(ret, cut(dr(trk, bestVertex), maxDr, std::less<float>()));
0258         if (ret == -1.f)
0259           return ret;
0260
0261         ret = std::min(ret, cut(dz(trk, bestVertex), maxDzcut, std::less<float>()));
0262         if (ret == -1.f)
0263           return ret;
0264       }
0265
0266       // parametrized dz and dr cut by using PV error
0267       if (dzWPVerr_par[2] < std::numeric_limits<float>::max() || drWPVerr_par[2] < std::numeric_limits<float>::max()) {
0268         Point bestVertexError(-1., -1., -1.);
0269         // min number of tracks [2 (=default) for offline, 3 for HLT]
0270         Point bestVertex = getBestVertex_withError(trk, vertices, bestVertexError, minNVtxTrk);
0271
0272         float maxDz_par[3];
0273         float maxDr_par[3];
0274         dzCut_wPVerror_par(trk, nLayers, dzWPVerr_par, dz_exp, bestVertexError, maxDz_par);
0275         drCut_wPVerror_par(trk, nLayers, drWPVerr_par, dr_exp, bestVertexError, maxDr_par);
0276
0277         ret = std::min(ret, cut(dr(trk, bestVertex), maxDr_par, std::less<float>()));
0278         if (ret == -1.f)
0279           return ret;
0280
0281         ret = std::min(ret, cut(dz(trk, bestVertex), maxDr_par, std::less<float>()));
0282         if (ret == -1.f)
0283           return ret;
0284       }
0285
0286       // parametrized dz and dr cut by using their error
0287       if (dz_par1[2] < std::numeric_limits<float>::max() || dr_par1[2] < std::numeric_limits<float>::max()) {
0288         float maxDz_par1[3];
0289         float maxDr_par1[3];
0290         dzCut_par1(trk, nLayers, dz_par1, dz_exp, maxDz_par1);
0291         drCut_par1(trk, nLayers, dr_par1, dr_exp, maxDr_par1);
0292
0293         float maxDz_par[3];
0294         float maxDr_par[3];
0295         std::copy(std::begin(maxDz_par1), std::end(maxDz_par1), std::begin(maxDz_par));
0296         std::copy(std::begin(maxDr_par1), std::end(maxDr_par1), std::begin(maxDr_par));
0297
0298         // parametrized dz and dr cut by using d0 and z0 resolution
0299         if (dz_par2[2] < std::numeric_limits<float>::max() || dr_par2[2] < std::numeric_limits<float>::max()) {
0300           float maxDz_par2[3];
0301           float maxDr_par2[3];
0302           dzCut_par2(trk, nLayers, dz_par2, dz_exp, d0err, d0err_par, maxDz_par2);
0303           drCut_par2(trk, nLayers, dr_par2, dr_exp, d0err, d0err_par, maxDr_par2);
0304
0305           for (int i = 2; i >= 0; --i) {
0306             if (maxDr_par2[i] < maxDr_par[i])
0307               maxDr_par[i] = maxDr_par2[i];
0308             if (maxDz_par2[i] < maxDz_par[i])
0309               maxDz_par[i] = maxDz_par2[i];
0310           }
0311         }
0312
0313         Point bestVertex = getBestVertex(trk, vertices, minNVtxTrk);  // min number of tracks 3 @HLT
0314         if (bestVertex.z() < -99998.) {
0315           bestVertex = beamSpot.position();
0316         }
0317
0318         ret = std::min(ret, cut(dz(trk, bestVertex), maxDz_par, std::less<float>()));
0319         if (ret == -1.f)
0320           return ret;
0321         ret = std::min(ret, cut(dr(trk, bestVertex), maxDr_par, std::less<float>()));
0322         if (ret == -1.f)
0323           return ret;
0324       }
0325       if (ret == -1.f)
0326         return ret;
0327
0328       return ret;
0329     }
0330
0331     static const char* name() { return "TrackCutClassifier"; }
0332
0333     static void fillDescriptions(edm::ParameterSetDescription& desc) {
0334       desc.add<bool>("isHLT", false);
0335       desc.add<std::vector<int>>(
0336           "minHits4pass",
0337           {std::numeric_limits<int>::max(), std::numeric_limits<int>::max(), std::numeric_limits<int>::max()});
0338       desc.add<std::vector<int>>("minHits", {0, 0, 1});
0339       desc.add<std::vector<int>>("minPixelHits", {0, 0, 1});
0340       desc.add<std::vector<int>>("minLayers", {3, 4, 5});
0341       desc.add<std::vector<int>>("min3DLayers", {1, 2, 3});
0342       desc.add<std::vector<int>>("maxLostLayers", {99, 3, 3});
0343       desc.add<std::vector<double>>(
0344           "maxRelPtErr",
0345           {std::numeric_limits<float>::max(), std::numeric_limits<float>::max(), std::numeric_limits<float>::max()});
0346       desc.add<std::vector<double>>("minNdof", {-1., -1., -1.});
0347       desc.add<std::vector<double>>("maxChi2", {9999., 25., 16.});
0348       desc.add<std::vector<double>>("maxChi2n", {9999., 1.0, 0.4});
0349
0350       desc.add<int>("minNVtxTrk", 2);
0351
0352       desc.add<std::vector<double>>(
0353           "maxDz",
0354           {std::numeric_limits<float>::max(), std::numeric_limits<float>::max(), std::numeric_limits<float>::max()});
0355       desc.add<std::vector<double>>("maxDzWrtBS", {std::numeric_limits<float>::max(), 24., 15.});
0356       desc.add<std::vector<double>>(
0357           "maxDr",
0358           {std::numeric_limits<float>::max(), std::numeric_limits<float>::max(), std::numeric_limits<float>::max()});
0359
0360       edm::ParameterSetDescription dz_par;
0361       dz_par.add<std::vector<int>>("dz_exp",
0362                                    {std::numeric_limits<int>::max(),
0363                                     std::numeric_limits<int>::max(),
0364                                     std::numeric_limits<int>::max()});  // par = 4
0365       dz_par.add<std::vector<double>>("dz_par1",
0366                                       {std::numeric_limits<float>::max(),
0367                                        std::numeric_limits<float>::max(),
0368                                        std::numeric_limits<float>::max()});  // par = 0.4
0369       dz_par.add<std::vector<double>>("dz_par2",
0370                                       {std::numeric_limits<float>::max(),
0371                                        std::numeric_limits<float>::max(),
0372                                        std::numeric_limits<float>::max()});  // par = 0.35
0373       dz_par.add<std::vector<double>>("dzWPVerr_par",
0374                                       {std::numeric_limits<float>::max(),
0375                                        std::numeric_limits<float>::max(),
0376                                        std::numeric_limits<float>::max()});  // par = 3.
0377       desc.add<edm::ParameterSetDescription>("dz_par", dz_par);
0378
0379       edm::ParameterSetDescription dr_par;
0380       dr_par.add<std::vector<int>>("dr_exp",
0381                                    {std::numeric_limits<int>::max(),
0382                                     std::numeric_limits<int>::max(),
0383                                     std::numeric_limits<int>::max()});  // par = 4
0384       dr_par.add<std::vector<double>>("dr_par1",
0385                                       {std::numeric_limits<float>::max(),
0386                                        std::numeric_limits<float>::max(),
0387                                        std::numeric_limits<float>::max()});  // par = 0.4
0388       dr_par.add<std::vector<double>>("dr_par2",
0389                                       {std::numeric_limits<float>::max(),
0390                                        std::numeric_limits<float>::max(),
0391                                        std::numeric_limits<float>::max()});  // par = 0.3
0392       dr_par.add<std::vector<double>>("d0err", {0.003, 0.003, 0.003});
0393       dr_par.add<std::vector<double>>("d0err_par", {0.001, 0.001, 0.001});
0394       dr_par.add<std::vector<double>>("drWPVerr_par",
0395                                       {std::numeric_limits<float>::max(),
0396                                        std::numeric_limits<float>::max(),
0397                                        std::numeric_limits<float>::max()});  // par = 3.
0398       desc.add<edm::ParameterSetDescription>("dr_par", dr_par);
0399     }
0400
0401     bool isHLT;
0402     float maxRelPtErr[3];
0403     float minNdof[3];
0404     float maxChi2[3];
0405     float maxChi2n[3];
0406     int minLayers[3];
0407     int min3DLayers[3];
0408     int minHits4pass[3];
0409     int minHits[3];
0410     int minPixelHits[3];
0411     int maxLostLayers[3];
0412     int minNVtxTrk;
0413     float maxDz[3];
0414     float maxDzWrtBS[3];
0415     float maxDr[3];
0416     int dz_exp[3];
0417     float dz_par1[3];
0418     float dz_par2[3];
0419     float dzWPVerr_par[3];
0420     int dr_exp[3];
0421     float dr_par1[3];
0422     float dr_par2[3];
0423     float d0err[3];
0424     float d0err_par[3];
0425     float drWPVerr_par[3];
0426   };
0427
0428   using TrackCutClassifier = TrackMVAClassifier<Cuts>;
0429
0430 }  // namespace
0431
0432 #include "FWCore/PluginManager/interface/ModuleDef.h"
0433 #include "FWCore/Framework/interface/MakerMacros.h"
0434
0435 DEFINE_FWK_MODULE(TrackCutClassifier);