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
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
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
0084 produces<MVACollection>("MVAValues");
0085
0086 for (unsigned int i = 0; i < trkSelectors.size(); i++) {
0087 qualityToSet_.push_back(TrackBase::undefQuality);
0088
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
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
0100 applyAdaptedPVCuts_.push_back(trkSelectors[i].getParameter<bool>("applyAdaptedPVCuts"));
0101
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
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
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 {
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();
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
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
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
0233 Handle<TrackingRecHitCollection> hSrcHits;
0234 evt.getByToken(hSrc_, hSrcHits);
0235 const TrackingRecHitCollection& srcHits(*hSrcHits);
0236
0237
0238 edm::Handle<reco::BeamSpot> hBsp;
0239 evt.getByToken(beamspot_, hBsp);
0240 const reco::BeamSpot& vertexBeamSpot(*hBsp);
0241
0242
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
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
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
0273
0274 LogTrace("TrackSelection") << "ready to check track with pt=" << trk.pt();
0275
0276
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
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
0344
0345 using namespace std;
0346
0347
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
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
0368
0369
0370
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
0402
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
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
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
0426 [[clang::suppress]]
0427 float lostMidFrac = tk.numberOfLostHits() / (tk.numberOfValidHits() + tk.numberOfLostHits());
0428 if (lostMidFrac > max_lostHitFraction_[tsNum])
0429 return false;
0430
0431
0432 float d0 = -tk.dxy(vertexBeamSpot.position()), d0E = tk.d0Error(), dz = tk.dz(vertexBeamSpot.position()),
0433 dzE = tk.dzError();
0434
0435
0436 float nomd0E = sqrt(res_par_[tsNum][0] * res_par_[tsNum][0] + (res_par_[tsNum][1] / pt) * (res_par_[tsNum][1] / pt));
0437
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
0446 bool primaryVertexZCompatibility(false);
0447 bool primaryVertexD0Compatibility(false);
0448
0449 if (points.empty()) {
0450
0451 if (abs(dz) < hypot(vertexBeamSpot.sigmaZ() * nSigmaZ_[tsNum], dzCut))
0452 primaryVertexZCompatibility = true;
0453
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);
0464 float d0PV = tk.dxy(*point);
0465 if (useVtxError_) {
0466 float dzErrPV = std::sqrt(dzE * dzE + vzerr[iv] * vzerr[iv]);
0467 float d0ErrPV = std::sqrt(d0E * d0E + vterr[iv] * vterr[iv]);
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
0489
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
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
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
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
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);