File indexing completed on 2025-01-08 03:36:32
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 #include <utility>
0012 #include <vector>
0013 #include <memory>
0014 #include <algorithm>
0015 #include <map>
0016
0017 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0018 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0019 #include "DataFormats/TrackReco/interface/Track.h"
0020 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0021 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0022 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h"
0023 #include "DataFormats/VertexReco/interface/Vertex.h"
0024 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0025 #include "FWCore/Framework/interface/Event.h"
0026 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0027 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0028 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0029 #include "FWCore/Utilities/interface/InputTag.h"
0030 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0031 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0032
0033 #include "MultiTrackSelector.h"
0034
0035 using namespace reco;
0036
0037 class dso_hidden AnalyticalTrackSelector final : public MultiTrackSelector {
0038 private:
0039 public:
0040
0041 explicit AnalyticalTrackSelector(const edm::ParameterSet& cfg);
0042
0043 ~AnalyticalTrackSelector() override = default;
0044
0045 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0046
0047 private:
0048 typedef math::XYZPoint Point;
0049
0050 void run(edm::Event& evt, const edm::EventSetup& es) const override;
0051
0052
0053 bool copyExtras_;
0054
0055 bool copyTrajectories_;
0056
0057 double minEta_;
0058 double maxEta_;
0059
0060 edm::EDGetTokenT<std::vector<Trajectory>> srcTraj_;
0061 edm::EDGetTokenT<TrajTrackAssociationCollection> srcTass_;
0062 };
0063
0064 #include "DataFormats/Common/interface/ValueMap.h"
0065
0066 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0067 #include <Math/DistFunc.h>
0068 #include "TMath.h"
0069
0070 AnalyticalTrackSelector::AnalyticalTrackSelector(const edm::ParameterSet& cfg) : MultiTrackSelector() {
0071
0072
0073 qualityToSet_.reserve(1);
0074 vtxNumber_.reserve(1);
0075 vertexCut_.reserve(1);
0076 res_par_.reserve(1);
0077 chi2n_par_.reserve(1);
0078 chi2n_no1Dmod_par_.reserve(1);
0079 d0_par1_.reserve(1);
0080 dz_par1_.reserve(1);
0081 d0_par2_.reserve(1);
0082 dz_par2_.reserve(1);
0083 applyAdaptedPVCuts_.reserve(1);
0084 max_d0_.reserve(1);
0085 max_z0_.reserve(1);
0086 nSigmaZ_.reserve(1);
0087 min_layers_.reserve(1);
0088 min_3Dlayers_.reserve(1);
0089 max_lostLayers_.reserve(1);
0090 min_hits_bypass_.reserve(1);
0091 applyAbsCutsIfNoPV_.reserve(1);
0092 max_d0NoPV_.reserve(1);
0093 max_z0NoPV_.reserve(1);
0094 preFilter_.reserve(1);
0095 max_relpterr_.reserve(1);
0096 min_nhits_.reserve(1);
0097 max_minMissHitOutOrIn_.reserve(1);
0098 max_lostHitFraction_.reserve(1);
0099 min_eta_.reserve(1);
0100 max_eta_.reserve(1);
0101 forest_.reserve(1);
0102 mvaType_.reserve(1);
0103 useMVA_.reserve(1);
0104
0105 produces<edm::ValueMap<float>>("MVAVals");
0106
0107 produces<MVACollection>("MVAValues");
0108 forest_[0] = nullptr;
0109 useAnyMVA_ = cfg.getParameter<bool>("useAnyMVA");
0110
0111 src_ = consumes<reco::TrackCollection>(cfg.getParameter<edm::InputTag>("src"));
0112 hSrc_ = consumes<TrackingRecHitCollection>(cfg.getParameter<edm::InputTag>("src"));
0113 beamspot_ = consumes<reco::BeamSpot>(cfg.getParameter<edm::InputTag>("beamspot"));
0114 useVertices_ = cfg.getParameter<bool>("useVertices");
0115 useVtxError_ = cfg.getParameter<bool>("useVtxError");
0116 if (useVertices_)
0117 vertices_ = consumes<reco::VertexCollection>(cfg.getParameter<edm::InputTag>("vertices"));
0118 copyExtras_ = cfg.getUntrackedParameter<bool>("copyExtras", false);
0119 copyTrajectories_ = cfg.getUntrackedParameter<bool>("copyTrajectories", false);
0120 if (copyTrajectories_) {
0121 srcTraj_ = consumes<std::vector<Trajectory>>(cfg.getParameter<edm::InputTag>("src"));
0122 srcTass_ = consumes<TrajTrackAssociationCollection>(cfg.getParameter<edm::InputTag>("src"));
0123 }
0124
0125 qualityToSet_.push_back(TrackBase::undefQuality);
0126
0127 vtxNumber_.push_back(useVertices_ ? cfg.getParameter<int32_t>("vtxNumber") : 0);
0128 vertexCut_.push_back(useVertices_ ? cfg.getParameter<std::string>("vertexCut") : "");
0129
0130 res_par_.push_back(cfg.getParameter<std::vector<double>>("res_par"));
0131 chi2n_par_.push_back(cfg.getParameter<double>("chi2n_par"));
0132 chi2n_no1Dmod_par_.push_back(cfg.getParameter<double>("chi2n_no1Dmod_par"));
0133 d0_par1_.push_back(cfg.getParameter<std::vector<double>>("d0_par1"));
0134 dz_par1_.push_back(cfg.getParameter<std::vector<double>>("dz_par1"));
0135 d0_par2_.push_back(cfg.getParameter<std::vector<double>>("d0_par2"));
0136 dz_par2_.push_back(cfg.getParameter<std::vector<double>>("dz_par2"));
0137
0138
0139 applyAdaptedPVCuts_.push_back(cfg.getParameter<bool>("applyAdaptedPVCuts"));
0140
0141 max_d0_.push_back(cfg.getParameter<double>("max_d0"));
0142 max_z0_.push_back(cfg.getParameter<double>("max_z0"));
0143 nSigmaZ_.push_back(cfg.getParameter<double>("nSigmaZ"));
0144
0145 min_layers_.push_back(cfg.getParameter<uint32_t>("minNumberLayers"));
0146 min_3Dlayers_.push_back(cfg.getParameter<uint32_t>("minNumber3DLayers"));
0147 max_lostLayers_.push_back(cfg.getParameter<uint32_t>("maxNumberLostLayers"));
0148 min_hits_bypass_.push_back(cfg.getParameter<uint32_t>("minHitsToBypassChecks"));
0149 max_relpterr_.push_back(cfg.getParameter<double>("max_relpterr"));
0150 min_nhits_.push_back(cfg.getParameter<uint32_t>("min_nhits"));
0151 max_minMissHitOutOrIn_.push_back(cfg.getParameter<int32_t>("max_minMissHitOutOrIn"));
0152 max_lostHitFraction_.push_back(cfg.getParameter<double>("max_lostHitFraction"));
0153 min_eta_.push_back(cfg.getParameter<double>("min_eta"));
0154 max_eta_.push_back(cfg.getParameter<double>("max_eta"));
0155
0156
0157 applyAbsCutsIfNoPV_.push_back(cfg.getParameter<bool>("applyAbsCutsIfNoPV"));
0158 keepAllTracks_.push_back(cfg.getParameter<bool>("keepAllTracks"));
0159
0160 setQualityBit_.push_back(false);
0161 std::string qualityStr = cfg.getParameter<std::string>("qualityBit");
0162
0163 if (d0_par1_[0].size() != 2 || dz_par1_[0].size() != 2 || d0_par2_[0].size() != 2 || dz_par2_[0].size() != 2) {
0164 edm::LogError("MisConfiguration") << "vector of size less then 2";
0165 throw;
0166 }
0167
0168 if (!qualityStr.empty()) {
0169 setQualityBit_[0] = true;
0170 qualityToSet_[0] = TrackBase::qualityByName(cfg.getParameter<std::string>("qualityBit"));
0171 }
0172
0173 if (keepAllTracks_[0] && !setQualityBit_[0])
0174 throw cms::Exception("Configuration")
0175 << "If you set 'keepAllTracks' to true, you must specify which qualityBit to set.\n";
0176 if (setQualityBit_[0] && (qualityToSet_[0] == TrackBase::undefQuality))
0177 throw cms::Exception("Configuration")
0178 << "You can't set the quality bit " << cfg.getParameter<std::string>("qualityBit")
0179 << " as it is 'undefQuality' or unknown.\n";
0180 if (applyAbsCutsIfNoPV_[0]) {
0181 max_d0NoPV_.push_back(cfg.getParameter<double>("max_d0NoPV"));
0182 max_z0NoPV_.push_back(cfg.getParameter<double>("max_z0NoPV"));
0183 } else {
0184 max_d0NoPV_.push_back(0.);
0185 max_z0NoPV_.push_back(0.);
0186 }
0187
0188 if (useAnyMVA_) {
0189 bool thisMVA = cfg.getParameter<bool>("useMVA");
0190 useMVA_.push_back(thisMVA);
0191 if (thisMVA) {
0192 double minVal = cfg.getParameter<double>("minMVA");
0193 min_MVA_.push_back(minVal);
0194 mvaType_.push_back(cfg.getParameter<std::string>("mvaType"));
0195 forestLabel_.push_back(cfg.getParameter<std::string>("GBRForestLabel"));
0196 useMVAonly_.push_back(cfg.getParameter<bool>("useMVAonly"));
0197 } else {
0198 min_MVA_.push_back(-9999.0);
0199 useMVAonly_.push_back(false);
0200 mvaType_.push_back("Detached");
0201 forestLabel_.push_back("MVASelectorIter0");
0202 }
0203 } else {
0204 useMVA_.push_back(false);
0205 useMVAonly_.push_back(false);
0206 min_MVA_.push_back(-9999.0);
0207 mvaType_.push_back("Detached");
0208 forestLabel_.push_back("MVASelectorIter0");
0209 }
0210
0211 std::string alias(cfg.getParameter<std::string>("@module_label"));
0212 if (copyExtras_) {
0213 produces<reco::TrackExtraCollection>().setBranchAlias(alias + "TrackExtras");
0214 produces<TrackingRecHitCollection>().setBranchAlias(alias + "RecHits");
0215 }
0216 if (copyTrajectories_) {
0217 produces<std::vector<Trajectory>>().setBranchAlias(alias + "Trajectories");
0218 produces<TrajTrackAssociationCollection>().setBranchAlias(alias + "TrajectoryTrackAssociations");
0219 }
0220
0221
0222
0223 produces<reco::TrackCollection>().setBranchAlias(alias + "Tracks");
0224 }
0225
0226 void AnalyticalTrackSelector::run(edm::Event& evt, const edm::EventSetup& es) const {
0227
0228 std::unique_ptr<reco::TrackCollection> selTracks_;
0229 std::unique_ptr<reco::TrackExtraCollection> selTrackExtras_;
0230 std::unique_ptr<TrackingRecHitCollection> selHits_;
0231 std::unique_ptr<std::vector<Trajectory>> selTrajs_;
0232 std::unique_ptr<std::vector<const Trajectory*>> selTrajPtrs_;
0233 std::unique_ptr<TrajTrackAssociationCollection> selTTAss_;
0234 reco::TrackRefProd rTracks_;
0235 reco::TrackExtraRefProd rTrackExtras_;
0236 TrackingRecHitRefProd rHits_;
0237 edm::RefProd<std::vector<Trajectory>> rTrajectories_;
0238 std::vector<reco::TrackRef> trackRefs_;
0239
0240 using namespace std;
0241 using namespace edm;
0242 using namespace reco;
0243
0244 Handle<TrackCollection> hSrcTrack;
0245 Handle<vector<Trajectory>> hTraj;
0246 Handle<vector<Trajectory>> hTrajP;
0247 Handle<TrajTrackAssociationCollection> hTTAss;
0248
0249
0250 edm::Handle<reco::BeamSpot> hBsp;
0251 evt.getByToken(beamspot_, hBsp);
0252 reco::BeamSpot vertexBeamSpot;
0253 vertexBeamSpot = *hBsp;
0254
0255
0256 const reco::VertexCollection dummyVtx;
0257 const reco::VertexCollection* vtxPtr = &dummyVtx;
0258 std::vector<Point> points;
0259 std::vector<float> vterr, vzerr;
0260 if (useVertices_) {
0261 vtxPtr = &evt.get(vertices_);
0262 selectVertices(0, *vtxPtr, points, vterr, vzerr);
0263
0264 LogDebug("SelectVertex") << points.size() << " good pixel vertices";
0265 }
0266
0267
0268 evt.getByToken(src_, hSrcTrack);
0269
0270 Handle<TrackingRecHitCollection> hSrcHits;
0271 evt.getByToken(hSrc_, hSrcHits);
0272 const TrackingRecHitCollection& srcHits(*hSrcHits);
0273
0274 selTracks_ = std::make_unique<TrackCollection>();
0275 rTracks_ = evt.getRefBeforePut<TrackCollection>();
0276 if (copyExtras_) {
0277 selTrackExtras_ = std::make_unique<TrackExtraCollection>();
0278 selHits_ = std::make_unique<TrackingRecHitCollection>();
0279 rHits_ = evt.getRefBeforePut<TrackingRecHitCollection>();
0280 rTrackExtras_ = evt.getRefBeforePut<TrackExtraCollection>();
0281 }
0282
0283 if (copyTrajectories_)
0284 trackRefs_.resize(hSrcTrack->size());
0285
0286 std::vector<float> mvaVals_(hSrcTrack->size(), -99.f);
0287 processMVA(evt, es, vertexBeamSpot, *vtxPtr, 0, mvaVals_, true);
0288
0289
0290 size_t current = 0;
0291 for (TrackCollection::const_iterator it = hSrcTrack->begin(), ed = hSrcTrack->end(); it != ed; ++it, ++current) {
0292 const Track& trk = *it;
0293
0294
0295 LogTrace("TrackSelection") << "ready to check track with pt=" << trk.pt();
0296
0297 float mvaVal = 0;
0298 if (useAnyMVA_)
0299 mvaVal = mvaVals_[current];
0300 bool ok = select(0, vertexBeamSpot, srcHits, trk, points, vterr, vzerr, mvaVal);
0301 if (!ok) {
0302 LogTrace("TrackSelection") << "track with pt=" << trk.pt() << " NOT selected";
0303
0304 if (copyTrajectories_)
0305 trackRefs_[current] = reco::TrackRef();
0306 if (!keepAllTracks_[0])
0307 continue;
0308 }
0309 LogTrace("TrackSelection") << "track with pt=" << trk.pt() << " selected";
0310 selTracks_->push_back(Track(trk));
0311 if (ok && setQualityBit_[0]) {
0312 selTracks_->back().setQuality(qualityToSet_[0]);
0313 if (qualityToSet_[0] == TrackBase::tight) {
0314 selTracks_->back().setQuality(TrackBase::loose);
0315 } else if (qualityToSet_[0] == TrackBase::highPurity) {
0316 selTracks_->back().setQuality(TrackBase::loose);
0317 selTracks_->back().setQuality(TrackBase::tight);
0318 }
0319 if (!points.empty()) {
0320 if (qualityToSet_[0] == TrackBase::loose) {
0321 selTracks_->back().setQuality(TrackBase::looseSetWithPV);
0322 } else if (qualityToSet_[0] == TrackBase::highPurity) {
0323 selTracks_->back().setQuality(TrackBase::looseSetWithPV);
0324 selTracks_->back().setQuality(TrackBase::highPuritySetWithPV);
0325 }
0326 }
0327 }
0328 if (copyExtras_) {
0329
0330 selTrackExtras_->push_back(TrackExtra(trk.outerPosition(),
0331 trk.outerMomentum(),
0332 trk.outerOk(),
0333 trk.innerPosition(),
0334 trk.innerMomentum(),
0335 trk.innerOk(),
0336 trk.outerStateCovariance(),
0337 trk.outerDetId(),
0338 trk.innerStateCovariance(),
0339 trk.innerDetId(),
0340 trk.seedDirection(),
0341 trk.seedRef()));
0342 selTracks_->back().setExtra(TrackExtraRef(rTrackExtras_, selTrackExtras_->size() - 1));
0343 TrackExtra& tx = selTrackExtras_->back();
0344 tx.setResiduals(trk.residuals());
0345
0346 auto const firstHitIndex = selHits_->size();
0347 for (trackingRecHit_iterator hit = trk.recHitsBegin(); hit != trk.recHitsEnd(); ++hit) {
0348 selHits_->push_back((*hit)->clone());
0349 }
0350 tx.setHits(rHits_, firstHitIndex, selHits_->size() - firstHitIndex);
0351 tx.setTrajParams(trk.extra()->trajParams(), trk.extra()->chi2sX5());
0352 }
0353 if (copyTrajectories_) {
0354 trackRefs_[current] = TrackRef(rTracks_, selTracks_->size() - 1);
0355 }
0356 }
0357 if (copyTrajectories_) {
0358 Handle<vector<Trajectory>> hTraj;
0359 Handle<TrajTrackAssociationCollection> hTTAss;
0360 evt.getByToken(srcTass_, hTTAss);
0361 evt.getByToken(srcTraj_, hTraj);
0362 selTrajs_ = std::make_unique<std::vector<Trajectory>>();
0363 rTrajectories_ = evt.getRefBeforePut<vector<Trajectory>>();
0364 selTTAss_ = std::make_unique<TrajTrackAssociationCollection>(rTrajectories_, rTracks_);
0365 for (size_t i = 0, n = hTraj->size(); i < n; ++i) {
0366 Ref<vector<Trajectory>> trajRef(hTraj, i);
0367 TrajTrackAssociationCollection::const_iterator match = hTTAss->find(trajRef);
0368 if (match != hTTAss->end()) {
0369 const Ref<TrackCollection>& trkRef = match->val;
0370 short oldKey = static_cast<short>(trkRef.key());
0371 if (trackRefs_[oldKey].isNonnull()) {
0372 selTrajs_->push_back(Trajectory(*trajRef));
0373 selTTAss_->insert(Ref<vector<Trajectory>>(rTrajectories_, selTrajs_->size() - 1), trackRefs_[oldKey]);
0374 }
0375 }
0376 }
0377 }
0378
0379 evt.put(std::move(selTracks_));
0380 if (copyExtras_) {
0381 evt.put(std::move(selTrackExtras_));
0382 evt.put(std::move(selHits_));
0383 }
0384 if (copyTrajectories_) {
0385 evt.put(std::move(selTrajs_));
0386 evt.put(std::move(selTTAss_));
0387 }
0388 }
0389
0390 void AnalyticalTrackSelector::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0391 edm::ParameterSetDescription desc;
0392
0393 desc.add<edm::InputTag>("src", edm::InputTag("generalTracks"));
0394 desc.add<bool>("keepAllTracks", false)
0395 ->setComment("if set to true tracks failing this filter are kept in the output");
0396 desc.add<edm::InputTag>("beamspot", edm::InputTag("offlineBeamSpot"));
0397
0398
0399 desc.add<bool>("useVertices", true);
0400 desc.add<bool>("useVtxError", false);
0401 desc.add<edm::InputTag>("vertices", edm::InputTag("firstStepPrimaryVertices"));
0402 desc.add<int32_t>("vtxNumber", -1);
0403 desc.add<std::string>("vertexCut", "ndof>=2&!isFake");
0404
0405 desc.addUntracked<bool>("copyExtras", false);
0406 desc.addUntracked<bool>("copyTrajectories", false);
0407 desc.add<std::string>("qualityBit", std::string(""))->setComment("set to ''if you don't want to set the bit");
0408
0409
0410 desc.add<double>("chi2n_no1Dmod_par", 9999.)
0411 ->setComment("parameter for adapted optimal cuts on chi2 and primary vertex compatibility");
0412 desc.add<double>("chi2n_par", 1.6)
0413 ->setComment("parameter for adapted optimal cuts on chi2 and primary vertex compatibility");
0414 desc.add<std::vector<double>>("res_par", {0.003, 0.01})->setComment("default: Loose");
0415 desc.add<std::vector<double>>("d0_par1", {0.55, 4.0})->setComment("default: Loose");
0416 desc.add<std::vector<double>>("d0_par2", {0.65, 4.0})->setComment("default: Loose");
0417 desc.add<std::vector<double>>("dz_par1", {0.55, 4.0})->setComment("default: Loose");
0418 desc.add<std::vector<double>>("dz_par2", {0.45, 4.0})->setComment("default: Loose");
0419 desc.add<bool>("applyAdaptedPVCuts", true)
0420 ->setComment("Boolean indicating if adapted primary vertex compatibility cuts are to be applied.");
0421
0422
0423 desc.add<double>("max_d0", 100.)->setComment("transverse impact parameter absolute cut");
0424 desc.add<double>("max_z0", 100.)->setComment("longitudinal impact parameter absolute cut");
0425 desc.add<double>("nSigmaZ", 4.);
0426
0427
0428 desc.add<uint32_t>("minNumberLayers", 0);
0429 desc.add<uint32_t>("minNumber3DLayers", 0);
0430 desc.add<uint32_t>("minHitsToBypassChecks", 20);
0431 desc.add<uint32_t>("maxNumberLostLayers", 999);
0432
0433
0434 desc.add<bool>("applyAbsCutsIfNoPV", false);
0435 desc.add<double>("max_d0NoPV", 100.);
0436 desc.add<double>("max_z0NoPV", 100.);
0437
0438
0439 desc.add<double>("max_relpterr", 9999)->setComment("parameter for cutting on pterror/pt");
0440 desc.add<uint32_t>("min_nhits", 0)->setComment("parameter for cutting on number of valid hits");
0441
0442 desc.add<double>("max_lostHitFraction", 1.0);
0443 desc.add<int32_t>("max_minMissHitOutOrIn", 99);
0444
0445
0446 desc.add<double>("max_eta", 9999.);
0447 desc.add<double>("min_eta", -9999.);
0448
0449
0450 desc.add<bool>("useMVA", false);
0451 desc.add<bool>("useAnyMVA", false);
0452 desc.add<bool>("useMVAonly", false);
0453 desc.add<double>("minMVA", -1.);
0454 desc.add<std::string>("GBRForestLabel", "MVASelectorIter0");
0455 desc.add<std::string>("mvaType", "Detached");
0456
0457 descriptions.addWithDefaultLabel(desc);
0458 }
0459
0460 #include "FWCore/PluginManager/interface/ModuleDef.h"
0461 #include "FWCore/Framework/interface/MakerMacros.h"
0462
0463 DEFINE_FWK_MODULE(AnalyticalTrackSelector);