Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:28:07 

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])
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   float lostMidFrac = tk.numberOfLostHits() / (tk.numberOfValidHits() + tk.numberOfLostHits());
0426   if (lostMidFrac > max_lostHitFraction_[tsNum])
0427     return false;
0428 
0429   //other track parameters
0430   float d0 = -tk.dxy(vertexBeamSpot.position()), d0E = tk.d0Error(), dz = tk.dz(vertexBeamSpot.position()),
0431         dzE = tk.dzError();
0432 
0433   // parametrized d0 resolution for the track pt
0434   float nomd0E = sqrt(res_par_[tsNum][0] * res_par_[tsNum][0] + (res_par_[tsNum][1] / pt) * (res_par_[tsNum][1] / pt));
0435   // parametrized z0 resolution for the track pt and eta
0436   float nomdzE = nomd0E * (std::cosh(eta));
0437 
0438   float dzCut = std::min(powN(dz_par1_[tsNum][0] * nlayers, int(dz_par1_[tsNum][1] + 0.5)) * nomdzE,
0439                          powN(dz_par2_[tsNum][0] * nlayers, int(dz_par2_[tsNum][1] + 0.5)) * dzE);
0440   float d0Cut = std::min(powN(d0_par1_[tsNum][0] * nlayers, int(d0_par1_[tsNum][1] + 0.5)) * nomd0E,
0441                          powN(d0_par2_[tsNum][0] * nlayers, int(d0_par2_[tsNum][1] + 0.5)) * d0E);
0442 
0443   // ---- PrimaryVertex compatibility cut
0444   bool primaryVertexZCompatibility(false);
0445   bool primaryVertexD0Compatibility(false);
0446 
0447   if (points.empty()) {  //If not primaryVertices are reconstructed, check just the compatibility with the BS
0448     //z0 within (n sigma + dzCut) of the beam spot z, if no good vertex is found
0449     if (abs(dz) < hypot(vertexBeamSpot.sigmaZ() * nSigmaZ_[tsNum], dzCut))
0450       primaryVertexZCompatibility = true;
0451     // d0 compatibility with beam line
0452     if (abs(d0) < d0Cut)
0453       primaryVertexD0Compatibility = true;
0454   }
0455 
0456   int iv = 0;
0457   for (std::vector<Point>::const_iterator point = points.begin(), end = points.end(); point != end; ++point) {
0458     LogTrace("TrackSelection") << "Test track w.r.t. vertex with z position " << point->z();
0459     if (primaryVertexZCompatibility && primaryVertexD0Compatibility)
0460       break;
0461     float dzPV = tk.dz(*point);   //re-evaluate the dz with respect to the vertex position
0462     float d0PV = tk.dxy(*point);  //re-evaluate the dxy with respect to the vertex position
0463     if (useVtxError_) {
0464       float dzErrPV = std::sqrt(dzE * dzE + vzerr[iv] * vzerr[iv]);  // include vertex error in z
0465       float d0ErrPV = std::sqrt(d0E * d0E + vterr[iv] * vterr[iv]);  // include vertex error in xy
0466       iv++;
0467       if (abs(dzPV) < dz_par1_[tsNum][0] * pow(nlayers, dz_par1_[tsNum][1]) * nomdzE &&
0468           abs(dzPV) < dz_par2_[tsNum][0] * pow(nlayers, dz_par2_[tsNum][1]) * dzErrPV && abs(dzPV) < max_z0_[tsNum])
0469         primaryVertexZCompatibility = true;
0470       if (abs(d0PV) < d0_par1_[tsNum][0] * pow(nlayers, d0_par1_[tsNum][1]) * nomd0E &&
0471           abs(d0PV) < d0_par2_[tsNum][0] * pow(nlayers, d0_par2_[tsNum][1]) * d0ErrPV && abs(d0PV) < max_d0_[tsNum])
0472         primaryVertexD0Compatibility = true;
0473     } else {
0474       if (abs(dzPV) < dzCut)
0475         primaryVertexZCompatibility = true;
0476       if (abs(d0PV) < d0Cut)
0477         primaryVertexD0Compatibility = true;
0478     }
0479     LogTrace("TrackSelection") << "distances " << dzPV << " " << d0PV << " vs " << dzCut << " " << d0Cut;
0480   }
0481 
0482   if (points.empty() && applyAbsCutsIfNoPV_[tsNum]) {
0483     if (abs(dz) > max_z0NoPV_[tsNum] || abs(d0) > max_d0NoPV_[tsNum])
0484       return false;
0485   } else {
0486     // Absolute cuts on all tracks impact parameters with respect to beam-spot.
0487     // If BS is not compatible, verify if at least the reco-vertex is compatible (useful for incorrect BS settings)
0488     if (abs(d0) > max_d0_[tsNum] && !primaryVertexD0Compatibility)
0489       return false;
0490     LogTrace("TrackSelection") << "absolute cuts on d0 passed";
0491     if (abs(dz) > max_z0_[tsNum] && !primaryVertexZCompatibility)
0492       return false;
0493     LogTrace("TrackSelection") << "absolute cuts on dz passed";
0494   }
0495 
0496   LogTrace("TrackSelection") << "cuts on PV: apply adapted PV cuts? " << applyAdaptedPVCuts_[tsNum]
0497                              << " d0 compatibility? " << primaryVertexD0Compatibility << " z compatibility? "
0498                              << primaryVertexZCompatibility;
0499 
0500   if (applyAdaptedPVCuts_[tsNum]) {
0501     return (primaryVertexD0Compatibility && primaryVertexZCompatibility);
0502   } else {
0503     return true;
0504   }
0505 }
0506 
0507 void MultiTrackSelector::selectVertices(unsigned int tsNum,
0508                                         const reco::VertexCollection& vtxs,
0509                                         std::vector<Point>& points,
0510                                         std::vector<float>& vterr,
0511                                         std::vector<float>& vzerr) const {
0512   // Select good primary vertices
0513   using namespace reco;
0514   int32_t toTake = vtxNumber_[tsNum];
0515   for (VertexCollection::const_iterator it = vtxs.begin(), ed = vtxs.end(); it != ed; ++it) {
0516     LogDebug("SelectVertex") << " select vertex with z position " << it->z() << " " << it->chi2() << " " << it->ndof()
0517                              << " " << TMath::Prob(it->chi2(), static_cast<int32_t>(it->ndof()));
0518     const Vertex& vtx = *it;
0519     bool pass = vertexCut_[tsNum](vtx);
0520     if (pass) {
0521       points.push_back(it->position());
0522       vterr.push_back(sqrt(it->yError() * it->xError()));
0523       vzerr.push_back(it->zError());
0524       LogTrace("SelectVertex") << " SELECTED vertex with z position " << it->z();
0525       toTake--;
0526       if (toTake == 0)
0527         break;
0528     }
0529   }
0530 }
0531 
0532 void MultiTrackSelector::processMVA(edm::Event& evt,
0533                                     const edm::EventSetup& es,
0534                                     const reco::BeamSpot& beamspot,
0535                                     const reco::VertexCollection& vertices,
0536                                     int selIndex,
0537                                     std::vector<float>& mvaVals_,
0538                                     bool writeIt) const {
0539   using namespace std;
0540   using namespace edm;
0541   using namespace reco;
0542 
0543   // Get tracks
0544   Handle<TrackCollection> hSrcTrack;
0545   evt.getByToken(src_, hSrcTrack);
0546   const TrackCollection& srcTracks(*hSrcTrack);
0547   RefToBaseProd<Track> rtbpTrackCollection(hSrcTrack);
0548   assert(mvaVals_.size() == srcTracks.size());
0549 
0550   // get hits in track..
0551   Handle<TrackingRecHitCollection> hSrcHits;
0552   evt.getByToken(hSrc_, hSrcHits);
0553   const TrackingRecHitCollection& srcHits(*hSrcHits);
0554 
0555   auto mvaValValueMap = std::make_unique<edm::ValueMap<float>>();
0556   edm::ValueMap<float>::Filler mvaFiller(*mvaValValueMap);
0557 
0558   if (!useAnyMVA_ && writeIt) {
0559     // mvaVals_ already initalized...
0560     mvaFiller.insert(hSrcTrack, mvaVals_.begin(), mvaVals_.end());
0561     mvaFiller.fill();
0562     evt.put(std::move(mvaValValueMap), "MVAVals");
0563     auto mvas = std::make_unique<MVACollection>(mvaVals_.begin(), mvaVals_.end());
0564     evt.put(std::move(mvas), "MVAValues");
0565     return;
0566   }
0567 
0568   if (!useMVA_[selIndex] && !writeIt)
0569     return;
0570 
0571   size_t current = 0;
0572   for (TrackCollection::const_iterator it = srcTracks.begin(), ed = srcTracks.end(); it != ed; ++it, ++current) {
0573     const Track& trk = *it;
0574     RefToBase<Track> trackRef(rtbpTrackCollection, current);
0575     auto tmva_ndof_ = trk.ndof();
0576     auto tmva_nlayers_ = trk.hitPattern().trackerLayersWithMeasurement();
0577     auto tmva_nlayers3D_ =
0578         trk.hitPattern().pixelLayersWithMeasurement() + trk.hitPattern().numberOfValidStripLayersWithMonoAndStereo();
0579     auto tmva_nlayerslost_ = trk.hitPattern().trackerLayersWithoutMeasurement(reco::HitPattern::TRACK_HITS);
0580     float chi2n = trk.normalizedChi2();
0581     float chi2n_no1Dmod = chi2n;
0582 
0583     int count1dhits = 0;
0584     auto ith = trk.extra()->firstRecHit();
0585     auto edh = ith + trk.recHitsSize();
0586     for (; ith < edh; ++ith) {
0587       const TrackingRecHit& hit = srcHits[ith];
0588       if (hit.dimension() == 1)
0589         ++count1dhits;
0590     }
0591     if (count1dhits > 0) {
0592       float chi2 = trk.chi2();
0593       float ndof = trk.ndof();
0594       chi2n = (chi2 + count1dhits) / float(ndof + count1dhits);
0595     }
0596     auto tmva_chi2n_ = chi2n;
0597     auto tmva_chi2n_no1dmod_ = chi2n_no1Dmod;
0598     auto tmva_eta_ = trk.eta();
0599     auto tmva_relpterr_ = float(trk.ptError()) / std::max(float(trk.pt()), 0.000001f);
0600     auto tmva_nhits_ = trk.numberOfValidHits();
0601     int lostIn = trk.hitPattern().numberOfLostTrackerHits(reco::HitPattern::MISSING_INNER_HITS);
0602     int lostOut = trk.hitPattern().numberOfLostTrackerHits(reco::HitPattern::MISSING_OUTER_HITS);
0603     auto tmva_minlost_ = std::min(lostIn, lostOut);
0604     auto tmva_lostmidfrac_ = trk.numberOfLostHits() / (trk.numberOfValidHits() + trk.numberOfLostHits());
0605     auto tmva_absd0_ = fabs(-trk.dxy(beamspot.position()));
0606     auto tmva_absdz_ = fabs(trk.dz(beamspot.position()));
0607     Point bestVertex = getBestVertex(trackRef, vertices);
0608     auto tmva_absd0PV_ = fabs(trk.dxy(bestVertex));
0609     auto tmva_absdzPV_ = fabs(trk.dz(bestVertex));
0610     auto tmva_pt_ = trk.pt();
0611 
0612     GBRForest const* forest = forest_[selIndex];
0613     if (useForestFromDB_) {
0614       forest = &es.getData(forestToken_[selIndex]);
0615     }
0616 
0617     float gbrVals_[16];
0618     gbrVals_[0] = tmva_pt_;
0619     gbrVals_[1] = tmva_lostmidfrac_;
0620     gbrVals_[2] = tmva_minlost_;
0621     gbrVals_[3] = tmva_nhits_;
0622     gbrVals_[4] = tmva_relpterr_;
0623     gbrVals_[5] = tmva_eta_;
0624     gbrVals_[6] = tmva_chi2n_no1dmod_;
0625     gbrVals_[7] = tmva_chi2n_;
0626     gbrVals_[8] = tmva_nlayerslost_;
0627     gbrVals_[9] = tmva_nlayers3D_;
0628     gbrVals_[10] = tmva_nlayers_;
0629     gbrVals_[11] = tmva_ndof_;
0630     gbrVals_[12] = tmva_absd0PV_;
0631     gbrVals_[13] = tmva_absdzPV_;
0632     gbrVals_[14] = tmva_absdz_;
0633     gbrVals_[15] = tmva_absd0_;
0634 
0635     if (mvaType_[selIndex] == "Prompt") {
0636       auto gbrVal = forest->GetClassifier(gbrVals_);
0637       mvaVals_[current] = gbrVal;
0638     } else {
0639       float detachedGbrVals_[12];
0640       for (int jjj = 0; jjj < 12; jjj++)
0641         detachedGbrVals_[jjj] = gbrVals_[jjj];
0642       auto gbrVal = forest->GetClassifier(detachedGbrVals_);
0643       mvaVals_[current] = gbrVal;
0644     }
0645   }
0646 
0647   if (writeIt) {
0648     mvaFiller.insert(hSrcTrack, mvaVals_.begin(), mvaVals_.end());
0649     mvaFiller.fill();
0650     evt.put(std::move(mvaValValueMap), "MVAVals");
0651     auto mvas = std::make_unique<MVACollection>(mvaVals_.begin(), mvaVals_.end());
0652     evt.put(std::move(mvas), "MVAValues");
0653   }
0654 }
0655 
0656 MultiTrackSelector::Point MultiTrackSelector::getBestVertex(TrackBaseRef track, VertexCollection vertices) const {
0657   Point p(0, 0, -99999);
0658   Point p_dz(0, 0, -99999);
0659   float bestWeight = 0;
0660   float dzmin = 10000;
0661   bool weightMatch = false;
0662   for (auto const& vertex : vertices) {
0663     float w = vertex.trackWeight(track);
0664     const Point& v_pos = vertex.position();
0665     if (w > bestWeight) {
0666       p = v_pos;
0667       bestWeight = w;
0668       weightMatch = true;
0669     }
0670     float dz = fabs(track.get()->dz(v_pos));
0671     if (dz < dzmin) {
0672       p_dz = v_pos;
0673       dzmin = dz;
0674     }
0675   }
0676   if (weightMatch)
0677     return p;
0678   else
0679     return p_dz;
0680 }
0681 
0682 #include "FWCore/PluginManager/interface/ModuleDef.h"
0683 #include "FWCore/Framework/interface/MakerMacros.h"
0684 
0685 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!