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
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);
0048 case 3:
0049 return PowN<3>::op(t);
0050 case 8:
0051 return PowN<8>::op(t);
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 }
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
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
0198 min_MVA_.reserve(trkSelectors.size());
0199
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
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
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
0217 applyAdaptedPVCuts_.push_back(trkSelectors[i].getParameter<bool>("applyAdaptedPVCuts"));
0218
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
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
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 {
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();
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
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
0322 Handle<TrackCollection> hSrcTrack;
0323 evt.getByToken(src_, hSrcTrack);
0324 const TrackCollection &srcTracks(*hSrcTrack);
0325
0326
0327 Handle<TrackingRecHitCollection> hSrcHits;
0328 evt.getByToken(hSrc_, hSrcHits);
0329 const TrackingRecHitCollection &srcHits(*hSrcHits);
0330
0331
0332 edm::Handle<reco::BeamSpot> hBsp;
0333 evt.getByToken(beamspot_, hBsp);
0334 const reco::BeamSpot &vertexBeamSpot(*hBsp);
0335
0336
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
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
0362
0363 LogTrace("TrackSelection") << "ready to check track with pt=" << trk.pt();
0364
0365
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
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
0430
0431 using namespace std;
0432
0433
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
0445
0446 if (useAnyMVA_ && useMVA_[tsNum]) {
0447 if (mvaVal < min_MVA_[tsNum])
0448 return false;
0449 else
0450 return true;
0451 }
0452
0453
0454
0455
0456
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
0488
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
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
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
0517 if (applyPixelMergingCuts_) {
0518
0519 if (pt < pixel_pTMinCut_[tsNum][0])
0520 return false;
0521 if (pt > pixel_pTMaxCut_[tsNum][0])
0522 return false;
0523
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
0535 float d0 = -tk.dxy(vertexBeamSpot.position()), d0E = tk.d0Error(), dz = tk.dz(vertexBeamSpot.position()),
0536 dzE = tk.dzError();
0537
0538
0539 float nomd0E = sqrt(res_par_[tsNum][0] * res_par_[tsNum][0] + (res_par_[tsNum][1] / pt) * (res_par_[tsNum][1] / pt));
0540
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
0549 bool primaryVertexZCompatibility(false);
0550 bool primaryVertexD0Compatibility(false);
0551
0552 if (points.empty()) {
0553
0554 if (abs(dz) < hypot(vertexBeamSpot.sigmaZ() * nSigmaZ_[tsNum], dzCut))
0555 primaryVertexZCompatibility = true;
0556
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);
0567 float d0PV = tk.dxy(*point);
0568 if (useVtxError_) {
0569 float dzErrPV = std::sqrt(dzE * dzE + vzerr[iv] * vzerr[iv]);
0570 float d0ErrPV = std::sqrt(d0E * d0E + vterr[iv] * vterr[iv]);
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
0592
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
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
0646 Handle<TrackCollection> hSrcTrack;
0647 evt.getByToken(src_, hSrcTrack);
0648 const TrackCollection &srcTracks(*hSrcTrack);
0649 assert(mvaVals_.size() == srcTracks.size());
0650
0651
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
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;
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;
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
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);