Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:18:14 

0001 #include "HIMultiTrackSelector.h"
0002 
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "DataFormats/Common/interface/ValueMap.h"
0005 #include "CondFormats/DataRecord/interface/GBRWrapperRcd.h"
0006 #include "FWCore/Framework/interface/EventSetup.h"
0007 
0008 #include <Math/DistFunc.h>
0009 #include <TMath.h>
0010 #include <TFile.h>
0011 
0012 namespace {
0013   // not a generic solution (wrong for N negative for instance)
0014   template <int N>
0015   struct PowN {
0016     template <typename T>
0017     static T op(T t) {
0018       return PowN<N / 2>::op(t) * PowN<(N + 1) / 2>::op(t);
0019     }
0020   };
0021   template <>
0022   struct PowN<0> {
0023     template <typename T>
0024     static T op(T t) {
0025       return T(1);
0026     }
0027   };
0028   template <>
0029   struct PowN<1> {
0030     template <typename T>
0031     static T op(T t) {
0032       return t;
0033     }
0034   };
0035   template <>
0036   struct PowN<2> {
0037     template <typename T>
0038     static T op(T t) {
0039       return t * t;
0040     }
0041   };
0042 
0043   template <typename T>
0044   T powN(T t, int n) {
0045     switch (n) {
0046       case 4:
0047         return PowN<4>::op(t);  // the only one that matters
0048       case 3:
0049         return PowN<3>::op(t);  // and this
0050       case 8:
0051         return PowN<8>::op(t);  // used in conversion????
0052       case 2:
0053         return PowN<2>::op(t);
0054       case 5:
0055         return PowN<5>::op(t);
0056       case 6:
0057         return PowN<6>::op(t);
0058       case 7:
0059         return PowN<7>::op(t);
0060       case 0:
0061         return PowN<0>::op(t);
0062       case 1:
0063         return PowN<1>::op(t);
0064       default:
0065         return std::pow(t, T(n));
0066     }
0067   }
0068 
0069 }  // namespace
0070 
0071 using namespace reco;
0072 
0073 HIMultiTrackSelector::HIMultiTrackSelector() {
0074   useForestFromDB_ = true;
0075   forest_ = nullptr;
0076 }
0077 
0078 void HIMultiTrackSelector::ParseForestVars() {
0079   mvavars_indices.clear();
0080   for (unsigned i = 0; i < forestVars_.size(); i++) {
0081     std::string v = forestVars_[i];
0082     int ind = -1;
0083     if (v == "chi2perdofperlayer")
0084       ind = chi2perdofperlayer;
0085     if (v == "dxyperdxyerror")
0086       ind = dxyperdxyerror;
0087     if (v == "dzperdzerror")
0088       ind = dzperdzerror;
0089     if (v == "relpterr")
0090       ind = relpterr;
0091     if (v == "lostmidfrac")
0092       ind = lostmidfrac;
0093     if (v == "minlost")
0094       ind = minlost;
0095     if (v == "nhits")
0096       ind = nhits;
0097     if (v == "eta")
0098       ind = eta;
0099     if (v == "chi2n_no1dmod")
0100       ind = chi2n_no1dmod;
0101     if (v == "chi2n")
0102       ind = chi2n;
0103     if (v == "nlayerslost")
0104       ind = nlayerslost;
0105     if (v == "nlayers3d")
0106       ind = nlayers3d;
0107     if (v == "nlayers")
0108       ind = nlayers;
0109     if (v == "ndof")
0110       ind = ndof;
0111     if (v == "etaerror")
0112       ind = etaerror;
0113 
0114     if (ind == -1)
0115       edm::LogWarning("HIMultiTrackSelector")
0116           << "Unknown forest variable " << v << ". Please make sure it's in the list of supported variables\n";
0117 
0118     mvavars_indices.push_back(ind);
0119   }
0120 }
0121 
0122 HIMultiTrackSelector::HIMultiTrackSelector(const edm::ParameterSet &cfg)
0123     : src_(consumes<reco::TrackCollection>(cfg.getParameter<edm::InputTag>("src"))),
0124       hSrc_(consumes<TrackingRecHitCollection>(cfg.getParameter<edm::InputTag>("src"))),
0125       beamspot_(consumes<reco::BeamSpot>(cfg.getParameter<edm::InputTag>("beamspot"))),
0126       useVertices_(cfg.getParameter<bool>("useVertices")),
0127       useVtxError_(cfg.getParameter<bool>("useVtxError"))
0128 // now get the pset for each selector
0129 {
0130   if (useVertices_)
0131     vertices_ = consumes<reco::VertexCollection>(cfg.getParameter<edm::InputTag>("vertices"));
0132 
0133   applyPixelMergingCuts_ = false;
0134   if (cfg.exists("applyPixelMergingCuts"))
0135     applyPixelMergingCuts_ = cfg.getParameter<bool>("applyPixelMergingCuts");
0136 
0137   useAnyMVA_ = false;
0138   forestLabel_ = "MVASelectorIter0";
0139   std::string type = "BDTG";
0140   useForestFromDB_ = true;
0141   dbFileName_ = "";
0142 
0143   forest_ = nullptr;
0144 
0145   if (cfg.exists("useAnyMVA"))
0146     useAnyMVA_ = cfg.getParameter<bool>("useAnyMVA");
0147   if (useAnyMVA_) {
0148     if (cfg.exists("mvaType"))
0149       type = cfg.getParameter<std::string>("mvaType");
0150     if (cfg.exists("GBRForestLabel"))
0151       forestLabel_ = cfg.getParameter<std::string>("GBRForestLabel");
0152     if (cfg.exists("GBRForestVars")) {
0153       forestVars_ = cfg.getParameter<std::vector<std::string>>("GBRForestVars");
0154       ParseForestVars();
0155     }
0156     if (cfg.exists("GBRForestFileName")) {
0157       dbFileName_ = cfg.getParameter<std::string>("GBRForestFileName");
0158       useForestFromDB_ = false;
0159     }
0160     mvaType_ = type;
0161   }
0162   if (useForestFromDB_) {
0163     forestToken_ = esConsumes(edm::ESInputTag("", forestLabel_));
0164   }
0165   std::vector<edm::ParameterSet> trkSelectors(cfg.getParameter<std::vector<edm::ParameterSet>>("trackSelectors"));
0166   qualityToSet_.reserve(trkSelectors.size());
0167   vtxNumber_.reserve(trkSelectors.size());
0168   vertexCut_.reserve(trkSelectors.size());
0169   res_par_.reserve(trkSelectors.size());
0170   chi2n_par_.reserve(trkSelectors.size());
0171   chi2n_no1Dmod_par_.reserve(trkSelectors.size());
0172   d0_par1_.reserve(trkSelectors.size());
0173   dz_par1_.reserve(trkSelectors.size());
0174   d0_par2_.reserve(trkSelectors.size());
0175   dz_par2_.reserve(trkSelectors.size());
0176   applyAdaptedPVCuts_.reserve(trkSelectors.size());
0177   max_d0_.reserve(trkSelectors.size());
0178   max_z0_.reserve(trkSelectors.size());
0179   nSigmaZ_.reserve(trkSelectors.size());
0180   pixel_pTMinCut_.reserve(trkSelectors.size());
0181   pixel_pTMaxCut_.reserve(trkSelectors.size());
0182   min_layers_.reserve(trkSelectors.size());
0183   min_3Dlayers_.reserve(trkSelectors.size());
0184   max_lostLayers_.reserve(trkSelectors.size());
0185   min_hits_bypass_.reserve(trkSelectors.size());
0186   applyAbsCutsIfNoPV_.reserve(trkSelectors.size());
0187   max_d0NoPV_.reserve(trkSelectors.size());
0188   max_z0NoPV_.reserve(trkSelectors.size());
0189   preFilter_.reserve(trkSelectors.size());
0190   max_relpterr_.reserve(trkSelectors.size());
0191   min_nhits_.reserve(trkSelectors.size());
0192   max_minMissHitOutOrIn_.reserve(trkSelectors.size());
0193   max_lostHitFraction_.reserve(trkSelectors.size());
0194   min_eta_.reserve(trkSelectors.size());
0195   max_eta_.reserve(trkSelectors.size());
0196   useMVA_.reserve(trkSelectors.size());
0197   //mvaReaders_.reserve(trkSelectors.size());
0198   min_MVA_.reserve(trkSelectors.size());
0199   //mvaType_.reserve(trkSelectors.size());
0200 
0201   produces<edm::ValueMap<float>>("MVAVals");
0202 
0203   for (unsigned int i = 0; i < trkSelectors.size(); i++) {
0204     qualityToSet_.push_back(TrackBase::undefQuality);
0205     // parameters for vertex selection
0206     vtxNumber_.push_back(useVertices_ ? trkSelectors[i].getParameter<int32_t>("vtxNumber") : 0);
0207     vertexCut_.push_back(useVertices_ ? trkSelectors[i].getParameter<std::string>("vertexCut") : nullptr);
0208     //  parameters for adapted optimal cuts on chi2 and primary vertex compatibility
0209     res_par_.push_back(trkSelectors[i].getParameter<std::vector<double>>("res_par"));
0210     chi2n_par_.push_back(trkSelectors[i].getParameter<double>("chi2n_par"));
0211     chi2n_no1Dmod_par_.push_back(trkSelectors[i].getParameter<double>("chi2n_no1Dmod_par"));
0212     d0_par1_.push_back(trkSelectors[i].getParameter<std::vector<double>>("d0_par1"));
0213     dz_par1_.push_back(trkSelectors[i].getParameter<std::vector<double>>("dz_par1"));
0214     d0_par2_.push_back(trkSelectors[i].getParameter<std::vector<double>>("d0_par2"));
0215     dz_par2_.push_back(trkSelectors[i].getParameter<std::vector<double>>("dz_par2"));
0216     // Boolean indicating if adapted primary vertex compatibility cuts are to be applied.
0217     applyAdaptedPVCuts_.push_back(trkSelectors[i].getParameter<bool>("applyAdaptedPVCuts"));
0218     // Impact parameter absolute cuts.
0219     max_d0_.push_back(trkSelectors[i].getParameter<double>("max_d0"));
0220     max_z0_.push_back(trkSelectors[i].getParameter<double>("max_z0"));
0221     nSigmaZ_.push_back(trkSelectors[i].getParameter<double>("nSigmaZ"));
0222     // Cuts on numbers of layers with hits/3D hits/lost hits.
0223     min_layers_.push_back(trkSelectors[i].getParameter<uint32_t>("minNumberLayers"));
0224     min_3Dlayers_.push_back(trkSelectors[i].getParameter<uint32_t>("minNumber3DLayers"));
0225     max_lostLayers_.push_back(trkSelectors[i].getParameter<uint32_t>("maxNumberLostLayers"));
0226     min_hits_bypass_.push_back(trkSelectors[i].getParameter<uint32_t>("minHitsToBypassChecks"));
0227     // Flag to apply absolute cuts if no PV passes the selection
0228     applyAbsCutsIfNoPV_.push_back(trkSelectors[i].getParameter<bool>("applyAbsCutsIfNoPV"));
0229     keepAllTracks_.push_back(trkSelectors[i].getParameter<bool>("keepAllTracks"));
0230     max_relpterr_.push_back(trkSelectors[i].getParameter<double>("max_relpterr"));
0231     min_nhits_.push_back(trkSelectors[i].getParameter<uint32_t>("min_nhits"));
0232     max_minMissHitOutOrIn_.push_back(trkSelectors[i].existsAs<int32_t>("max_minMissHitOutOrIn")
0233                                          ? trkSelectors[i].getParameter<int32_t>("max_minMissHitOutOrIn")
0234                                          : 99);
0235     max_lostHitFraction_.push_back(trkSelectors[i].existsAs<double>("max_lostHitFraction")
0236                                        ? trkSelectors[i].getParameter<double>("max_lostHitFraction")
0237                                        : 1.0);
0238     min_eta_.push_back(trkSelectors[i].existsAs<double>("min_eta") ? trkSelectors[i].getParameter<double>("min_eta")
0239                                                                    : -9999);
0240     max_eta_.push_back(trkSelectors[i].existsAs<double>("max_eta") ? trkSelectors[i].getParameter<double>("max_eta")
0241                                                                    : 9999);
0242 
0243     setQualityBit_.push_back(false);
0244     std::string qualityStr = trkSelectors[i].getParameter<std::string>("qualityBit");
0245     if (!qualityStr.empty()) {
0246       setQualityBit_[i] = true;
0247       qualityToSet_[i] = TrackBase::qualityByName(trkSelectors[i].getParameter<std::string>("qualityBit"));
0248     }
0249 
0250     if (setQualityBit_[i] && (qualityToSet_[i] == TrackBase::undefQuality))
0251       throw cms::Exception("Configuration")
0252           << "You can't set the quality bit " << trkSelectors[i].getParameter<std::string>("qualityBit")
0253           << " as it is 'undefQuality' or unknown.\n";
0254 
0255     if (applyAbsCutsIfNoPV_[i]) {
0256       max_d0NoPV_.push_back(trkSelectors[i].getParameter<double>("max_d0NoPV"));
0257       max_z0NoPV_.push_back(trkSelectors[i].getParameter<double>("max_z0NoPV"));
0258     } else {  //dummy values
0259       max_d0NoPV_.push_back(0.);
0260       max_z0NoPV_.push_back(0.);
0261     }
0262 
0263     name_.push_back(trkSelectors[i].getParameter<std::string>("name"));
0264 
0265     preFilter_[i] = trkSelectors.size();  // no prefilter
0266 
0267     std::string pfName = trkSelectors[i].getParameter<std::string>("preFilterName");
0268     if (!pfName.empty()) {
0269       bool foundPF = false;
0270       for (unsigned int j = 0; j < i; j++)
0271         if (name_[j] == pfName) {
0272           foundPF = true;
0273           preFilter_[i] = j;
0274         }
0275       if (!foundPF)
0276         throw cms::Exception("Configuration") << "Invalid prefilter name in HIMultiTrackSelector "
0277                                               << trkSelectors[i].getParameter<std::string>("preFilterName");
0278     }
0279 
0280     if (applyPixelMergingCuts_) {
0281       pixel_pTMinCut_.push_back(trkSelectors[i].getParameter<std::vector<double>>("pixel_pTMinCut"));
0282       pixel_pTMaxCut_.push_back(trkSelectors[i].getParameter<std::vector<double>>("pixel_pTMaxCut"));
0283     }
0284 
0285     //    produces<std::vector<int> >(name_[i]).setBranchAlias( name_[i] + "TrackQuals");
0286     produces<edm::ValueMap<int>>(name_[i]).setBranchAlias(name_[i] + "TrackQuals");
0287     if (useAnyMVA_) {
0288       bool thisMVA = false;
0289       if (trkSelectors[i].exists("useMVA"))
0290         thisMVA = trkSelectors[i].getParameter<bool>("useMVA");
0291       useMVA_.push_back(thisMVA);
0292       if (thisMVA) {
0293         double minVal = -1;
0294         if (trkSelectors[i].exists("minMVA"))
0295           minVal = trkSelectors[i].getParameter<double>("minMVA");
0296         min_MVA_.push_back(minVal);
0297 
0298       } else {
0299         min_MVA_.push_back(-9999.0);
0300       }
0301     } else {
0302       min_MVA_.push_back(-9999.0);
0303     }
0304   }
0305 }
0306 
0307 HIMultiTrackSelector::~HIMultiTrackSelector() { delete forest_; }
0308 
0309 void HIMultiTrackSelector::beginStream(edm::StreamID) {
0310   if (!useForestFromDB_) {
0311     TFile gbrfile(dbFileName_.c_str());
0312     forest_ = (GBRForest *)gbrfile.Get(forestLabel_.c_str());
0313   }
0314 }
0315 
0316 void HIMultiTrackSelector::run(edm::Event &evt, const edm::EventSetup &es) const {
0317   using namespace std;
0318   using namespace edm;
0319   using namespace reco;
0320 
0321   // Get tracks
0322   Handle<TrackCollection> hSrcTrack;
0323   evt.getByToken(src_, hSrcTrack);
0324   const TrackCollection &srcTracks(*hSrcTrack);
0325 
0326   // get hits in track..
0327   Handle<TrackingRecHitCollection> hSrcHits;
0328   evt.getByToken(hSrc_, hSrcHits);
0329   const TrackingRecHitCollection &srcHits(*hSrcHits);
0330 
0331   // looking for the beam spot
0332   edm::Handle<reco::BeamSpot> hBsp;
0333   evt.getByToken(beamspot_, hBsp);
0334   const reco::BeamSpot &vertexBeamSpot(*hBsp);
0335 
0336   // Select good primary vertices for use in subsequent track selection
0337   edm::Handle<reco::VertexCollection> hVtx;
0338   if (useVertices_)
0339     evt.getByToken(vertices_, hVtx);
0340 
0341   unsigned int trkSize = srcTracks.size();
0342   std::vector<int> selTracksSave(qualityToSet_.size() * trkSize, 0);
0343 
0344   std::vector<float> mvaVals_(srcTracks.size(), -99.f);
0345   processMVA(evt, es, mvaVals_, *hVtx);
0346 
0347   for (unsigned int i = 0; i < qualityToSet_.size(); i++) {
0348     std::vector<int> selTracks(trkSize, 0);
0349     auto selTracksValueMap = std::make_unique<edm::ValueMap<int>>();
0350     edm::ValueMap<int>::Filler filler(*selTracksValueMap);
0351 
0352     std::vector<Point> points;
0353     std::vector<float> vterr, vzerr;
0354     if (useVertices_)
0355       selectVertices(i, *hVtx, points, vterr, vzerr);
0356 
0357     // Loop over tracks
0358     size_t current = 0;
0359     for (TrackCollection::const_iterator it = srcTracks.begin(), ed = srcTracks.end(); it != ed; ++it, ++current) {
0360       const Track &trk = *it;
0361       // Check if this track passes cuts
0362 
0363       LogTrace("TrackSelection") << "ready to check track with pt=" << trk.pt();
0364 
0365       //already removed
0366       bool ok = true;
0367       float mvaVal = 0;
0368       if (preFilter_[i] < i && selTracksSave[preFilter_[i] * trkSize + current] < 0) {
0369         selTracks[current] = -1;
0370         ok = false;
0371         if (!keepAllTracks_[i])
0372           continue;
0373       } else {
0374         if (useAnyMVA_)
0375           mvaVal = mvaVals_[current];
0376         ok = select(i, vertexBeamSpot, srcHits, trk, points, vterr, vzerr, mvaVal);
0377         if (!ok) {
0378           LogTrace("TrackSelection") << "track with pt=" << trk.pt() << " NOT selected";
0379           if (!keepAllTracks_[i]) {
0380             selTracks[current] = -1;
0381             continue;
0382           }
0383         } else
0384           LogTrace("TrackSelection") << "track with pt=" << trk.pt() << " selected";
0385       }
0386 
0387       if (preFilter_[i] < i) {
0388         selTracks[current] = selTracksSave[preFilter_[i] * trkSize + current];
0389       } else {
0390         selTracks[current] = trk.qualityMask();
0391       }
0392       if (ok && setQualityBit_[i]) {
0393         selTracks[current] = (selTracks[current] | (1 << qualityToSet_[i]));
0394         if (qualityToSet_[i] == TrackBase::tight) {
0395           selTracks[current] = (selTracks[current] | (1 << TrackBase::loose));
0396         } else if (qualityToSet_[i] == TrackBase::highPurity) {
0397           selTracks[current] = (selTracks[current] | (1 << TrackBase::loose));
0398           selTracks[current] = (selTracks[current] | (1 << TrackBase::tight));
0399         }
0400 
0401         if (!points.empty()) {
0402           if (qualityToSet_[i] == TrackBase::loose) {
0403             selTracks[current] = (selTracks[current] | (1 << TrackBase::looseSetWithPV));
0404           } else if (qualityToSet_[i] == TrackBase::highPurity) {
0405             selTracks[current] = (selTracks[current] | (1 << TrackBase::looseSetWithPV));
0406             selTracks[current] = (selTracks[current] | (1 << TrackBase::highPuritySetWithPV));
0407           }
0408         }
0409       }
0410     }
0411     for (unsigned int j = 0; j < trkSize; j++)
0412       selTracksSave[j + i * trkSize] = selTracks[j];
0413     filler.insert(hSrcTrack, selTracks.begin(), selTracks.end());
0414     filler.fill();
0415 
0416     //    evt.put(std::move(selTracks),name_[i]);
0417     evt.put(std::move(selTracksValueMap), name_[i]);
0418   }
0419 }
0420 
0421 bool HIMultiTrackSelector::select(unsigned int tsNum,
0422                                   const reco::BeamSpot &vertexBeamSpot,
0423                                   const TrackingRecHitCollection &recHits,
0424                                   const reco::Track &tk,
0425                                   const std::vector<Point> &points,
0426                                   std::vector<float> &vterr,
0427                                   std::vector<float> &vzerr,
0428                                   double mvaVal) const {
0429   // Decide if the given track passes selection cuts.
0430 
0431   using namespace std;
0432 
0433   //cuts on number of valid hits
0434   auto nhits = tk.numberOfValidHits();
0435   if (nhits >= min_hits_bypass_[tsNum])
0436     return true;
0437   if (nhits < min_nhits_[tsNum])
0438     return false;
0439 
0440   if (tk.ndof() < 1E-5)
0441     return false;
0442 
0443   //////////////////////////////////////////////////
0444   //Adding the MVA selection before any other cut//
0445   ////////////////////////////////////////////////
0446   if (useAnyMVA_ && useMVA_[tsNum]) {
0447     if (mvaVal < min_MVA_[tsNum])
0448       return false;
0449     else
0450       return true;
0451   }
0452   /////////////////////////////////
0453   //End of MVA selection section//
0454   ///////////////////////////////
0455 
0456   // Cuts on numbers of layers with hits/3D hits/lost hits.
0457   uint32_t nlayers = tk.hitPattern().trackerLayersWithMeasurement();
0458   uint32_t nlayers3D =
0459       tk.hitPattern().pixelLayersWithMeasurement() + tk.hitPattern().numberOfValidStripLayersWithMonoAndStereo();
0460   uint32_t nlayersLost = tk.hitPattern().trackerLayersWithoutMeasurement(reco::HitPattern::TRACK_HITS);
0461   LogDebug("TrackSelection") << "cuts on nlayers: " << nlayers << " " << nlayers3D << " " << nlayersLost << " vs "
0462                              << min_layers_[tsNum] << " " << min_3Dlayers_[tsNum] << " " << max_lostLayers_[tsNum];
0463   if (nlayers < min_layers_[tsNum])
0464     return false;
0465   if (nlayers3D < min_3Dlayers_[tsNum])
0466     return false;
0467   if (nlayersLost > max_lostLayers_[tsNum])
0468     return false;
0469   LogTrace("TrackSelection") << "cuts on nlayers passed";
0470 
0471   float chi2n = tk.normalizedChi2();
0472   float chi2n_no1Dmod = chi2n;
0473 
0474   int count1dhits = 0;
0475   auto ith = tk.extra()->firstRecHit();
0476   auto edh = ith + tk.recHitsSize();
0477   for (; ith < edh; ++ith) {
0478     const TrackingRecHit &hit = recHits[ith];
0479     if (hit.dimension() == 1)
0480       ++count1dhits;
0481   }
0482   if (count1dhits > 0) {
0483     float chi2 = tk.chi2();
0484     float ndof = tk.ndof();
0485     chi2n = (chi2 + count1dhits) / float(ndof + count1dhits);
0486   }
0487   // For each 1D rechit, the chi^2 and ndof is increased by one.  This is a way of retaining approximately
0488   // the same normalized chi^2 distribution as with 2D rechits.
0489   if (chi2n > chi2n_par_[tsNum] * nlayers)
0490     return false;
0491 
0492   if (chi2n_no1Dmod > chi2n_no1Dmod_par_[tsNum] * nlayers)
0493     return false;
0494 
0495   // Get track parameters
0496   float pt = std::max(float(tk.pt()), 0.000001f);
0497   float eta = tk.eta();
0498   if (eta < min_eta_[tsNum] || eta > max_eta_[tsNum])
0499     return false;
0500 
0501   //cuts on relative error on pt
0502   float relpterr = float(tk.ptError()) / pt;
0503   if (relpterr > max_relpterr_[tsNum])
0504     return false;
0505 
0506   int lostIn = tk.hitPattern().numberOfLostTrackerHits(reco::HitPattern::MISSING_INNER_HITS);
0507   int lostOut = tk.hitPattern().numberOfLostTrackerHits(reco::HitPattern::MISSING_OUTER_HITS);
0508   int minLost = std::min(lostIn, lostOut);
0509   if (minLost > max_minMissHitOutOrIn_[tsNum])
0510     return false;
0511   float lostMidFrac =
0512       tk.numberOfLostHits() == 0 ? 0. : tk.numberOfLostHits() / (tk.numberOfValidHits() + tk.numberOfLostHits());
0513   if (lostMidFrac > max_lostHitFraction_[tsNum])
0514     return false;
0515 
0516   // Pixel Track Merging pT dependent cuts
0517   if (applyPixelMergingCuts_) {
0518     // hard cut at absolute min/max pt
0519     if (pt < pixel_pTMinCut_[tsNum][0])
0520       return false;
0521     if (pt > pixel_pTMaxCut_[tsNum][0])
0522       return false;
0523     // tapering cuts with chi2n_no1Dmod
0524     double pTMaxCutPos = (pixel_pTMaxCut_[tsNum][0] - pt) / (pixel_pTMaxCut_[tsNum][0] - pixel_pTMaxCut_[tsNum][1]);
0525     double pTMinCutPos = (pt - pixel_pTMinCut_[tsNum][0]) / (pixel_pTMinCut_[tsNum][1] - pixel_pTMinCut_[tsNum][0]);
0526     if (pt > pixel_pTMaxCut_[tsNum][1] &&
0527         chi2n_no1Dmod > pixel_pTMaxCut_[tsNum][2] * nlayers * pow(pTMaxCutPos, pixel_pTMaxCut_[tsNum][3]))
0528       return false;
0529     if (pt < pixel_pTMinCut_[tsNum][1] &&
0530         chi2n_no1Dmod > pixel_pTMinCut_[tsNum][2] * nlayers * pow(pTMinCutPos, pixel_pTMinCut_[tsNum][3]))
0531       return false;
0532   }
0533 
0534   //other track parameters
0535   float d0 = -tk.dxy(vertexBeamSpot.position()), d0E = tk.d0Error(), dz = tk.dz(vertexBeamSpot.position()),
0536         dzE = tk.dzError();
0537 
0538   // parametrized d0 resolution for the track pt
0539   float nomd0E = sqrt(res_par_[tsNum][0] * res_par_[tsNum][0] + (res_par_[tsNum][1] / pt) * (res_par_[tsNum][1] / pt));
0540   // parametrized z0 resolution for the track pt and eta
0541   float nomdzE = nomd0E * (std::cosh(eta));
0542 
0543   float dzCut = std::min(powN(dz_par1_[tsNum][0] * nlayers, int(dz_par1_[tsNum][1] + 0.5)) * nomdzE,
0544                          powN(dz_par2_[tsNum][0] * nlayers, int(dz_par2_[tsNum][1] + 0.5)) * dzE);
0545   float d0Cut = std::min(powN(d0_par1_[tsNum][0] * nlayers, int(d0_par1_[tsNum][1] + 0.5)) * nomd0E,
0546                          powN(d0_par2_[tsNum][0] * nlayers, int(d0_par2_[tsNum][1] + 0.5)) * d0E);
0547 
0548   // ---- PrimaryVertex compatibility cut
0549   bool primaryVertexZCompatibility(false);
0550   bool primaryVertexD0Compatibility(false);
0551 
0552   if (points.empty()) {  //If not primaryVertices are reconstructed, check just the compatibility with the BS
0553     //z0 within (n sigma + dzCut) of the beam spot z, if no good vertex is found
0554     if (abs(dz) < hypot(vertexBeamSpot.sigmaZ() * nSigmaZ_[tsNum], dzCut))
0555       primaryVertexZCompatibility = true;
0556     // d0 compatibility with beam line
0557     if (abs(d0) < d0Cut)
0558       primaryVertexD0Compatibility = true;
0559   }
0560 
0561   int iv = 0;
0562   for (std::vector<Point>::const_iterator point = points.begin(), end = points.end(); point != end; ++point) {
0563     LogTrace("TrackSelection") << "Test track w.r.t. vertex with z position " << point->z();
0564     if (primaryVertexZCompatibility && primaryVertexD0Compatibility)
0565       break;
0566     float dzPV = tk.dz(*point);   //re-evaluate the dz with respect to the vertex position
0567     float d0PV = tk.dxy(*point);  //re-evaluate the dxy with respect to the vertex position
0568     if (useVtxError_) {
0569       float dzErrPV = std::sqrt(dzE * dzE + vzerr[iv] * vzerr[iv]);  // include vertex error in z
0570       float d0ErrPV = std::sqrt(d0E * d0E + vterr[iv] * vterr[iv]);  // include vertex error in xy
0571       iv++;
0572       if (abs(dzPV) < dz_par1_[tsNum][0] * pow(nlayers, dz_par1_[tsNum][1]) * nomdzE &&
0573           abs(dzPV) < dz_par2_[tsNum][0] * pow(nlayers, dz_par2_[tsNum][1]) * dzErrPV && abs(dzPV) < max_z0_[tsNum])
0574         primaryVertexZCompatibility = true;
0575       if (abs(d0PV) < d0_par1_[tsNum][0] * pow(nlayers, d0_par1_[tsNum][1]) * nomd0E &&
0576           abs(d0PV) < d0_par2_[tsNum][0] * pow(nlayers, d0_par2_[tsNum][1]) * d0ErrPV && abs(d0PV) < max_d0_[tsNum])
0577         primaryVertexD0Compatibility = true;
0578     } else {
0579       if (abs(dzPV) < dzCut)
0580         primaryVertexZCompatibility = true;
0581       if (abs(d0PV) < d0Cut)
0582         primaryVertexD0Compatibility = true;
0583     }
0584     LogTrace("TrackSelection") << "distances " << dzPV << " " << d0PV << " vs " << dzCut << " " << d0Cut;
0585   }
0586 
0587   if (points.empty() && applyAbsCutsIfNoPV_[tsNum]) {
0588     if (abs(dz) > max_z0NoPV_[tsNum] || abs(d0) > max_d0NoPV_[tsNum])
0589       return false;
0590   } else {
0591     // Absolute cuts on all tracks impact parameters with respect to beam-spot.
0592     // If BS is not compatible, verify if at least the reco-vertex is compatible (useful for incorrect BS settings)
0593     if (abs(d0) > max_d0_[tsNum] && !primaryVertexD0Compatibility)
0594       return false;
0595     LogTrace("TrackSelection") << "absolute cuts on d0 passed";
0596     if (abs(dz) > max_z0_[tsNum] && !primaryVertexZCompatibility)
0597       return false;
0598     LogTrace("TrackSelection") << "absolute cuts on dz passed";
0599   }
0600 
0601   LogTrace("TrackSelection") << "cuts on PV: apply adapted PV cuts? " << applyAdaptedPVCuts_[tsNum]
0602                              << " d0 compatibility? " << primaryVertexD0Compatibility << " z compatibility? "
0603                              << primaryVertexZCompatibility;
0604 
0605   if (applyAdaptedPVCuts_[tsNum]) {
0606     return (primaryVertexD0Compatibility && primaryVertexZCompatibility);
0607   } else {
0608     return true;
0609   }
0610 }
0611 
0612 void HIMultiTrackSelector::selectVertices(unsigned int tsNum,
0613                                           const reco::VertexCollection &vtxs,
0614                                           std::vector<Point> &points,
0615                                           std::vector<float> &vterr,
0616                                           std::vector<float> &vzerr) const {
0617   // Select good primary vertices
0618   using namespace reco;
0619   int32_t toTake = vtxNumber_[tsNum];
0620   for (VertexCollection::const_iterator it = vtxs.begin(), ed = vtxs.end(); it != ed; ++it) {
0621     LogDebug("SelectVertex") << " select vertex with z position " << it->z() << " " << it->chi2() << " " << it->ndof()
0622                              << " " << TMath::Prob(it->chi2(), static_cast<int32_t>(it->ndof()));
0623     const Vertex &vtx = *it;
0624     bool pass = vertexCut_[tsNum](vtx);
0625     if (pass) {
0626       points.push_back(it->position());
0627       vterr.push_back(sqrt(it->yError() * it->xError()));
0628       vzerr.push_back(it->zError());
0629       LogTrace("SelectVertex") << " SELECTED vertex with z position " << it->z();
0630       toTake--;
0631       if (toTake == 0)
0632         break;
0633     }
0634   }
0635 }
0636 
0637 void HIMultiTrackSelector::processMVA(edm::Event &evt,
0638                                       const edm::EventSetup &es,
0639                                       std::vector<float> &mvaVals_,
0640                                       const reco::VertexCollection &vertices) const {
0641   using namespace std;
0642   using namespace edm;
0643   using namespace reco;
0644 
0645   // Get tracks
0646   Handle<TrackCollection> hSrcTrack;
0647   evt.getByToken(src_, hSrcTrack);
0648   const TrackCollection &srcTracks(*hSrcTrack);
0649   assert(mvaVals_.size() == srcTracks.size());
0650 
0651   // get hits in track..
0652   Handle<TrackingRecHitCollection> hSrcHits;
0653   evt.getByToken(hSrc_, hSrcHits);
0654   const TrackingRecHitCollection &srcHits(*hSrcHits);
0655 
0656   auto mvaValValueMap = std::make_unique<edm::ValueMap<float>>();
0657   edm::ValueMap<float>::Filler mvaFiller(*mvaValValueMap);
0658 
0659   if (!useAnyMVA_) {
0660     // mvaVals_ already initalized...
0661     mvaFiller.insert(hSrcTrack, mvaVals_.begin(), mvaVals_.end());
0662     mvaFiller.fill();
0663     evt.put(std::move(mvaValValueMap), "MVAVals");
0664     return;
0665   }
0666 
0667   bool checkvertex =
0668       std::find(mvavars_indices.begin(), mvavars_indices.end(), dxyperdxyerror) != mvavars_indices.end() ||
0669       std::find(mvavars_indices.begin(), mvavars_indices.end(), dzperdzerror) != mvavars_indices.end();
0670 
0671   size_t current = 0;
0672   for (TrackCollection::const_iterator it = srcTracks.begin(), ed = srcTracks.end(); it != ed; ++it, ++current) {
0673     const Track &trk = *it;
0674 
0675     float mvavalues[15];
0676     mvavalues[ndof] = trk.ndof();
0677     mvavalues[nlayers] = trk.hitPattern().trackerLayersWithMeasurement();
0678     mvavalues[nlayers3d] =
0679         trk.hitPattern().pixelLayersWithMeasurement() + trk.hitPattern().numberOfValidStripLayersWithMonoAndStereo();
0680     mvavalues[nlayerslost] = trk.hitPattern().trackerLayersWithoutMeasurement(reco::HitPattern::TRACK_HITS);
0681     mvavalues[chi2n_no1dmod] = trk.normalizedChi2();
0682     mvavalues[chi2perdofperlayer] = mvavalues[chi2n_no1dmod] / mvavalues[nlayers];
0683 
0684     float chi2n1d = trk.normalizedChi2();
0685     int count1dhits = 0;
0686     auto ith = trk.extra()->firstRecHit();
0687     auto edh = ith + trk.recHitsSize();
0688     for (; ith < edh; ++ith) {
0689       const TrackingRecHit &hit = srcHits[ith];
0690       if (hit.dimension() == 1)
0691         ++count1dhits;
0692     }
0693     if (count1dhits > 0) {
0694       float chi2 = trk.chi2();
0695       float ndof = trk.ndof();
0696       chi2n1d = (chi2 + count1dhits) / float(ndof + count1dhits);
0697     }
0698 
0699     mvavalues[chi2n] = chi2n1d;  //chi2 and 1d modes
0700 
0701     mvavalues[eta] = trk.eta();
0702     mvavalues[relpterr] = float(trk.ptError()) / std::max(float(trk.pt()), 0.000001f);
0703     mvavalues[nhits] = trk.numberOfValidHits();
0704 
0705     int lostIn = trk.hitPattern().numberOfLostTrackerHits(reco::HitPattern::MISSING_INNER_HITS);
0706     int lostOut = trk.hitPattern().numberOfLostTrackerHits(reco::HitPattern::MISSING_OUTER_HITS);
0707     mvavalues[minlost] = std::min(lostIn, lostOut);
0708     mvavalues[lostmidfrac] = trk.numberOfLostHits() / (trk.numberOfValidHits() + trk.numberOfLostHits());
0709 
0710     mvavalues[etaerror] = trk.etaError();
0711 
0712     float reldz = 0;
0713     float reldxy = 0;
0714     if (checkvertex) {
0715       int vtxind = 0;  // only first vertex is taken into account for the speed purposes
0716       float dxy = trk.dxy(vertices[vtxind].position()),
0717             dxyE = sqrt(trk.dxyError() * trk.dxyError() + vertices[vtxind].xError() * vertices[vtxind].yError());
0718       float dz = trk.dz(vertices[vtxind].position()),
0719             dzE = sqrt(trk.dzError() * trk.dzError() + vertices[vtxind].zError() * vertices[vtxind].zError());
0720       reldz = dz / dzE;
0721       reldxy = dxy / dxyE;
0722     }
0723     mvavalues[dxyperdxyerror] = reldxy;
0724     mvavalues[dzperdzerror] = reldz;
0725 
0726     std::vector<float> gbrValues;
0727 
0728     //fill in the gbrValues vector with the necessary variables
0729     for (unsigned i = 0; i < mvavars_indices.size(); i++) {
0730       gbrValues.push_back(mvavalues[mvavars_indices[i]]);
0731     }
0732 
0733     GBRForest const *forest = forest_;
0734     if (useForestFromDB_) {
0735       forest = &es.getData(forestToken_);
0736     }
0737 
0738     auto gbrVal = forest->GetClassifier(&gbrValues[0]);
0739     mvaVals_[current] = gbrVal;
0740   }
0741   mvaFiller.insert(hSrcTrack, mvaVals_.begin(), mvaVals_.end());
0742   mvaFiller.fill();
0743   evt.put(std::move(mvaValValueMap), "MVAVals");
0744 }
0745 
0746 #include "FWCore/PluginManager/interface/ModuleDef.h"
0747 #include "FWCore/Framework/interface/MakerMacros.h"
0748 
0749 DEFINE_FWK_MODULE(HIMultiTrackSelector);
[ 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!