Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-10-08 23:10:00 

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