File indexing completed on 2024-09-07 04:37:58
0001 #include "FWCore/Framework/interface/stream/EDProducer.h"
0002 #include "FWCore/Framework/interface/Event.h"
0003 #include "FWCore/Framework/interface/EventSetup.h"
0004 #include "FWCore/Framework/interface/ESHandle.h"
0005 #include "FWCore/Utilities/interface/TypeDemangler.h"
0006
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0010 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0011 #include "FWCore/Utilities/interface/InputTag.h"
0012 #include "FWCore/Utilities/interface/ESGetToken.h"
0013
0014 #include "DataFormats/Common/interface/View.h"
0015 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0016 #include "DataFormats/TrackReco/interface/Track.h"
0017 #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
0018 #include "DataFormats/TrackingRecHit/interface/InvalidTrackingRecHit.h"
0019
0020 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0021 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0022 #include "MagneticField/Engine/interface/MagneticField.h"
0023 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0024 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0025 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0026
0027 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0028 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0029 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0030
0031
0032 #include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
0033 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0034 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0035 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h"
0036 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0037 #include "RecoLocalTracker/SiStripClusterizer/interface/SiStripClusterInfo.h"
0038
0039 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0040 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0041 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0042 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0043 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0044 #include "TMath.h"
0045
0046 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0047 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0048 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
0049
0050 #include "DataFormats/Provenance/interface/RunLumiEventNumber.h"
0051
0052 #include <boost/regex.hpp>
0053 #include <map>
0054 #include <optional>
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081 namespace reco {
0082
0083 namespace modules {
0084 class TrackerTrackHitFilter : public edm::stream::EDProducer<> {
0085 public:
0086 TrackerTrackHitFilter(const edm::ParameterSet &iConfig);
0087 void produce(edm::Event &iEvent, const edm::EventSetup &iSetup) override;
0088 int checkHit(const edm::EventSetup &iSetup, const DetId &detid, const TrackingRecHit *hit);
0089 void produceFromTrajectory(const edm::EventSetup &iSetup,
0090 const Trajectory *itt,
0091 std::vector<TrackingRecHit *> &hits);
0092 void produceFromTrack(const edm::EventSetup &iSetup, const Track *itt, std::vector<TrackingRecHit *> &hits);
0093 unsigned int getSequLayer(const reco::Track &tk, unsigned int prevSequLayers, std::vector<bool> isNotValidVec);
0094 bool isFirstValidHitInLayer(const reco::Track &tk, std::vector<bool> isNotValidVec);
0095
0096 static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0097
0098 private:
0099 class Rule {
0100 public:
0101
0102 Rule(const std::string &str);
0103
0104 void apply(DetId detid, const TrackerTopology *tTopo, bool &verdict) const {
0105
0106 if (detid.subdetId() == subdet_) {
0107
0108 if ((layer_ == 0) || (layer_ == layer(detid, tTopo))) {
0109
0110 verdict = keep_;
0111
0112 }
0113
0114 }
0115
0116 }
0117
0118 private:
0119 int subdet_;
0120
0121 int layer_;
0122 bool keep_;
0123 int layer(DetId detid, const TrackerTopology *tTopo) const;
0124 };
0125
0126 edm::InputTag src_;
0127
0128 edm::RunNumber_t iRun;
0129 edm::EventNumber_t iEvt;
0130
0131 size_t minimumHits_;
0132
0133 unsigned int layersRemaining_;
0134 bool truncateTracks_;
0135
0136 bool replaceWithInactiveHits_;
0137 bool stripFrontInvalidHits_;
0138 bool stripBackInvalidHits_;
0139 bool stripAllInvalidHits_;
0140
0141 bool isPhase2_;
0142 bool rejectBadStoNHits_;
0143 std::string CMNSubtractionMode_;
0144 std::vector<bool> subdetStoN_;
0145 std::vector<double> subdetStoNlowcut_;
0146 std::vector<double> subdetStoNhighcut_;
0147 bool checkStoN(const DetId &id, const TrackingRecHit *therechit);
0148 void parseStoN(const std::string &str);
0149
0150 std::vector<uint32_t> detsToIgnore_;
0151 std::vector<Rule> rules_;
0152
0153 bool useTrajectories_;
0154 bool rejectLowAngleHits_;
0155 double TrackAngleCut_;
0156 bool checkHitAngle(const TrajectoryMeasurement &meas);
0157 bool checkPXLQuality_;
0158 double pxlTPLProbXY_;
0159 double pxlTPLProbXYQ_;
0160 std::vector<int32_t> pxlTPLqBin_;
0161
0162 bool checkPXLCorrClustCharge(const TrajectoryMeasurement &meas);
0163 double PXLcorrClusChargeCut_;
0164
0165 edm::ESHandle<TrackerGeometry> theGeometry;
0166 edm::ESHandle<MagneticField> theMagField;
0167
0168 edm::EDGetTokenT<reco::TrackCollection> tokenTracks;
0169 edm::EDGetTokenT<TrajTrackAssociationCollection> tokenTrajTrack;
0170 edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tokenGeometry;
0171 edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tokenMagField;
0172 edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tokenTrackerTopo;
0173
0174 std::optional<SiStripClusterInfo> siStripClusterInfo_;
0175
0176 bool tagOverlaps_;
0177 int nOverlaps;
0178 int layerFromId(const DetId &id, const TrackerTopology *tTopo) const;
0179 int sideFromId(const DetId &id, const TrackerTopology *tTopo) const;
0180
0181
0182 TrackCandidate makeCandidate(const reco::Track &tk,
0183 std::vector<TrackingRecHit *>::iterator hitsBegin,
0184 std::vector<TrackingRecHit *>::iterator hitsEnd);
0185
0186 };
0187
0188 TrackerTrackHitFilter::Rule::Rule(const std::string &str) {
0189 static const boost::regex rule("(keep|drop)\\s+([A-Z]+)(\\s+(\\d+))?");
0190 boost::cmatch match;
0191 std::string match_1;
0192 std::string match_2;
0193 std::string match_3;
0194
0195 if (!regex_match(str.c_str(), match, rule)) {
0196 throw cms::Exception("Configuration") << "Rule '" << str << "' not understood.\n";
0197 } else {
0198 edm::LogInfo("TrackerTrackHitFilter") << "*** Rule Command given to TrackerTrackHitFilter:\t" << str;
0199 }
0200
0201
0202 keep_ = (strncmp(match[1].first, "keep", 4) == 0);
0203
0204
0205 subdet_ = -1;
0206 if (strncmp(match[2].first, "PXB", 3) == 0)
0207 subdet_ = PixelSubdetector::PixelBarrel;
0208 else if (strncmp(match[2].first, "PXE", 3) == 0)
0209 subdet_ = PixelSubdetector::PixelEndcap;
0210 else if (strncmp(match[2].first, "TIB", 3) == 0)
0211 subdet_ = StripSubdetector::TIB;
0212 else if (strncmp(match[2].first, "TID", 3) == 0)
0213 subdet_ = StripSubdetector::TID;
0214 else if (strncmp(match[2].first, "TOB", 3) == 0)
0215 subdet_ = StripSubdetector::TOB;
0216 else if (strncmp(match[2].first, "TEC", 3) == 0)
0217 subdet_ = StripSubdetector::TEC;
0218 if (subdet_ == -1) {
0219 throw cms::Exception("Configuration")
0220 << "Detector '" << match[2].first << "' not understood. Should be PXB, PXE, TIB, TID, TOB, TEC.\n";
0221 }
0222
0223 if (match[4].first != match[4].second) {
0224 layer_ = atoi(match[4].first);
0225 } else {
0226 layer_ = 0;
0227 }
0228 }
0229
0230 int TrackerTrackHitFilter::Rule::layer(DetId detid, const TrackerTopology *tTopo) const {
0231 return tTopo->layer(detid);
0232 }
0233
0234 void TrackerTrackHitFilter::parseStoN(const std::string &str) {
0235
0236
0237
0238
0239 static const boost::regex rule(
0240 "([A-Z]+)"
0241 "\\s*(\\d+\\.*\\d*)?"
0242 "\\s*(\\d+\\.*\\d*)?");
0243
0244 boost::cmatch match;
0245 std::string match_1;
0246 std::string match_2;
0247 std::string match_3;
0248
0249 if (!regex_match(str.c_str(), match, rule)) {
0250 throw cms::Exception("Configuration") << "Rule for S to N cut '" << str << "' not understood.\n";
0251 } else {
0252 std::string match_0 = match[0].second;
0253 match_1 = match[1].second;
0254 match_2 = match[2].second;
0255 match_3 = match[3].second;
0256 }
0257
0258 int cnt = 0;
0259 float subdet_ind[6];
0260 for (cnt = 0; cnt < 6; cnt++) {
0261 subdet_ind[cnt] = -1.0;
0262 }
0263
0264 bool doALL = false;
0265 std::string match_1a(match[1].first, match[1].second);
0266 if (strncmp(match[1].first, "ALL", 3) == 0)
0267 doALL = true;
0268 if (doALL || strncmp(match[1].first, "PXB", 3) == 0)
0269 subdet_ind[0] = +1.0;
0270 if (doALL || strncmp(match[1].first, "PXE", 3) == 0)
0271 subdet_ind[1] = +1.0;
0272 if (doALL || strncmp(match[1].first, "TIB", 3) == 0)
0273 subdet_ind[2] = +1.0;
0274 if (doALL || strncmp(match[1].first, "TID", 3) == 0)
0275 subdet_ind[3] = +1.0;
0276 if (doALL || strncmp(match[1].first, "TOB", 3) == 0)
0277 subdet_ind[4] = +1.0;
0278 if (doALL || strncmp(match[1].first, "TEC", 3) == 0)
0279 subdet_ind[5] = +1.0;
0280
0281 for (cnt = 0; cnt < 6; cnt++) {
0282 if (subdet_ind[cnt] > 0.0) {
0283 subdetStoN_[cnt] = true;
0284 if (match[2].first != match[2].second) {
0285 subdetStoNlowcut_[cnt] = atof(match[2].first);
0286 }
0287 if (match[3].first != match[3].second) {
0288 subdetStoNhighcut_[cnt] = atof(match[3].first);
0289 }
0290 edm::LogInfo("TrackerTrackHitFilter") << "Setting thresholds*&^ for subdet #" << cnt + 1 << " = "
0291 << subdetStoNlowcut_[cnt] << " - " << subdetStoNhighcut_[cnt];
0292 }
0293 }
0294
0295 bool correct_regex = false;
0296 for (cnt = 0; cnt < 6; cnt++) {
0297 if (subdetStoN_[cnt])
0298 correct_regex = true;
0299 }
0300
0301 if (!correct_regex) {
0302 throw cms::Exception("Configuration")
0303 << "Detector '" << match_1a << "' not understood in parseStoN. Should be PXB, PXE, TIB, TID, TOB, TEC.\n";
0304 }
0305
0306
0307 }
0308
0309 TrackerTrackHitFilter::TrackerTrackHitFilter(const edm::ParameterSet &iConfig)
0310 : src_(iConfig.getParameter<edm::InputTag>("src")),
0311 minimumHits_(iConfig.getParameter<uint32_t>("minimumHits")),
0312 layersRemaining_(iConfig.getParameter<uint32_t>("layersRemaining")),
0313 truncateTracks_(iConfig.getParameter<bool>("truncateTracks")),
0314 replaceWithInactiveHits_(iConfig.getParameter<bool>("replaceWithInactiveHits")),
0315 stripFrontInvalidHits_(iConfig.getParameter<bool>("stripFrontInvalidHits")),
0316 stripBackInvalidHits_(iConfig.getParameter<bool>("stripBackInvalidHits")),
0317 stripAllInvalidHits_(iConfig.getParameter<bool>("stripAllInvalidHits")),
0318 isPhase2_(iConfig.getParameter<bool>("isPhase2")),
0319 rejectBadStoNHits_(iConfig.getParameter<bool>("rejectBadStoNHits")),
0320 CMNSubtractionMode_(iConfig.getParameter<std::string>("CMNSubtractionMode")),
0321 detsToIgnore_(iConfig.getParameter<std::vector<uint32_t> >("detsToIgnore")),
0322 useTrajectories_(iConfig.getParameter<bool>("useTrajectories")),
0323 rejectLowAngleHits_(iConfig.getParameter<bool>("rejectLowAngleHits")),
0324 TrackAngleCut_(iConfig.getParameter<double>("TrackAngleCut")),
0325 checkPXLQuality_(iConfig.getParameter<bool>("usePixelQualityFlag")),
0326 pxlTPLProbXY_(iConfig.getParameter<double>("PxlTemplateProbXYCut")),
0327 pxlTPLProbXYQ_(iConfig.getParameter<double>("PxlTemplateProbXYChargeCut")),
0328 pxlTPLqBin_(iConfig.getParameter<std::vector<int32_t> >("PxlTemplateqBinCut")),
0329 PXLcorrClusChargeCut_(iConfig.getParameter<double>("PxlCorrClusterChargeCut")),
0330 tagOverlaps_(iConfig.getParameter<bool>("tagOverlaps")) {
0331
0332
0333 if (!isPhase2_) {
0334 siStripClusterInfo_ = consumesCollector();
0335 }
0336
0337 tokenGeometry = esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>();
0338 tokenMagField = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0339 tokenTrackerTopo = esConsumes<TrackerTopology, TrackerTopologyRcd>();
0340 if (useTrajectories_)
0341 tokenTrajTrack = consumes<TrajTrackAssociationCollection>(src_);
0342 else
0343 tokenTracks = consumes<reco::TrackCollection>(src_);
0344
0345
0346 if (stripAllInvalidHits_ && replaceWithInactiveHits_) {
0347 throw cms::Exception("Configuration") << "Inconsistent Configuration: you can't set both 'stripAllInvalidHits' "
0348 "and 'replaceWithInactiveHits' to true\n";
0349 }
0350 if (rejectLowAngleHits_ && !useTrajectories_) {
0351 throw cms::Exception("Configuration") << "Wrong configuration of TrackerTrackHitFilter. You cannot apply the "
0352 "cut on the track angle without using Trajectories!\n";
0353 }
0354 if (!useTrajectories_ && PXLcorrClusChargeCut_ > 0) {
0355 throw cms::Exception("Configuration")
0356 << "Wrong configuration of TrackerTrackHitFilter. You cannot apply the cut on the corrected pixel cluster "
0357 "charge without using Trajectories!\n";
0358 }
0359
0360 if (pxlTPLqBin_.size() > 2) {
0361 edm::LogInfo("TrackerTrackHitFIlter") << "Warning from TrackerTrackHitFilter: vector with qBin cuts has size > "
0362 "2. Additional items will be ignored.";
0363 }
0364
0365
0366 std::vector<std::string> str_rules = iConfig.getParameter<std::vector<std::string> >("commands");
0367 rules_.reserve(str_rules.size());
0368 for (std::vector<std::string>::const_iterator it = str_rules.begin(), ed = str_rules.end(); it != ed; ++it) {
0369 rules_.push_back(Rule(*it));
0370 }
0371
0372 if (rejectBadStoNHits_) {
0373
0374 subdetStoN_.reserve(6);
0375 subdetStoNlowcut_.reserve(6);
0376 subdetStoNhighcut_.reserve(6);
0377 int cnt = 0;
0378 for (cnt = 0; cnt < 6; cnt++) {
0379 subdetStoN_[cnt] = false;
0380 subdetStoNlowcut_[cnt] = -1.0;
0381 subdetStoNhighcut_[cnt] = -1.0;
0382 }
0383
0384 std::vector<std::string> str_StoNrules = iConfig.getParameter<std::vector<std::string> >("StoNcommands");
0385 for (std::vector<std::string>::const_iterator str_StoN = str_StoNrules.begin(); str_StoN != str_StoNrules.end();
0386 ++str_StoN) {
0387 parseStoN(*str_StoN);
0388 }
0389
0390 edm::LogInfo("TrackerTrackHitFilter") << "Finished parsing S/N. Applying following cuts to subdets:";
0391 for (cnt = 0; cnt < 6; cnt++) {
0392
0393 edm::LogVerbatim("TrackerTrackHitFilter")
0394 << "Subdet #" << cnt + 1 << " -> " << subdetStoNlowcut_[cnt] << " , " << subdetStoNhighcut_[cnt];
0395 }
0396 }
0397
0398 if (rejectLowAngleHits_)
0399 edm::LogInfo("TrackerTrackHitFilter") << "\nApplying cut on angle track = " << TrackAngleCut_;
0400
0401
0402 std::sort(detsToIgnore_.begin(), detsToIgnore_.end());
0403
0404
0405 produces<TrackCandidateCollection>();
0406 }
0407
0408 void TrackerTrackHitFilter::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0409
0410 iRun = iEvent.id().run();
0411 iEvt = iEvent.id().event();
0412
0413
0414 edm::Handle<std::vector<reco::Track> > tracks;
0415 edm::Handle<TrajTrackAssociationCollection> assoMap;
0416
0417 if (useTrajectories_)
0418 iEvent.getByToken(tokenTrajTrack, assoMap);
0419 else
0420 iEvent.getByToken(tokenTracks, tracks);
0421
0422
0423 theGeometry = iSetup.getHandle(tokenGeometry);
0424 theMagField = iSetup.getHandle(tokenMagField);
0425 if (!isPhase2_)
0426 siStripClusterInfo_->initEvent(iSetup);
0427
0428
0429 size_t candcollsize;
0430 if (useTrajectories_)
0431 candcollsize = assoMap->size();
0432 else
0433 candcollsize = tracks->size();
0434 auto output = std::make_unique<TrackCandidateCollection>();
0435
0436 output->reserve(candcollsize);
0437
0438
0439 std::vector<TrackingRecHit *> hits;
0440
0441 if (useTrajectories_) {
0442 for (TrajTrackAssociationCollection::const_iterator itass = assoMap->begin(); itass != assoMap->end();
0443 ++itass) {
0444 const edm::Ref<std::vector<Trajectory> > traj = itass->key;
0445 const reco::TrackRef tkref = itass->val;
0446
0447
0448 const Track *trk = &(*tkref);
0449 const Trajectory *myTrajectory = &(*traj);
0450 produceFromTrajectory(iSetup, myTrajectory, hits);
0451
0452 std::vector<TrackingRecHit *>::iterator begin = hits.begin(), end = hits.end();
0453
0454
0455 if (stripFrontInvalidHits_) {
0456 while ((begin != end) && ((*begin)->isValid() == false))
0457 ++begin;
0458 }
0459
0460 if (stripBackInvalidHits_ && (begin != end)) {
0461 --end;
0462 while ((begin != end) && ((*end)->isValid() == false))
0463 --end;
0464 ++end;
0465 }
0466
0467
0468 if (replaceWithInactiveHits_) {
0469 int nvalidhits = 0;
0470 for (std::vector<TrackingRecHit *>::iterator ithit = begin; ithit != end; ++ithit) {
0471 if ((*ithit)->isValid())
0472 nvalidhits++;
0473 }
0474 if (nvalidhits >= int(minimumHits_)) {
0475 output->push_back(makeCandidate(*trk, begin, end));
0476 }
0477 } else {
0478 if ((end - begin) >= int(minimumHits_)) {
0479 output->push_back(makeCandidate(*trk, begin, end));
0480 }
0481 }
0482
0483
0484
0485
0486
0487
0488
0489 for (begin = hits.begin(), end = hits.end(); begin != end; ++begin) {
0490 if (*begin)
0491 delete *begin;
0492 }
0493 hits.clear();
0494 }
0495
0496 } else {
0497
0498
0499 for (std::vector<reco::Track>::const_iterator ittrk = tracks->begin(), edtrk = tracks->end(); ittrk != edtrk;
0500 ++ittrk) {
0501
0502
0503 const Track *trk = &(*ittrk);
0504
0505 produceFromTrack(iSetup, trk, hits);
0506
0507
0508
0509
0510
0511
0512
0513
0514
0515
0516
0517
0518
0519
0520
0521
0522
0523
0524
0525
0526
0527
0528 std::vector<TrackingRecHit *>::iterator begin = hits.begin(), end = hits.end();
0529
0530
0531 if (stripFrontInvalidHits_) {
0532 while ((begin != end) && ((*begin)->isValid() == false))
0533 ++begin;
0534 }
0535
0536 if (stripBackInvalidHits_ && (begin != end)) {
0537 --end;
0538 while ((begin != end) && ((*end)->isValid() == false))
0539 --end;
0540 ++end;
0541 }
0542
0543
0544 if (replaceWithInactiveHits_) {
0545 int nvalidhits = 0;
0546 for (std::vector<TrackingRecHit *>::iterator ithit = begin; ithit != end; ++ithit) {
0547 if ((*ithit)->isValid())
0548 nvalidhits++;
0549 }
0550 if (nvalidhits >= int(minimumHits_)) {
0551 output->push_back(makeCandidate(*ittrk, begin, end));
0552 }
0553
0554 } else {
0555 if ((end - begin) >= int(minimumHits_)) {
0556 output->push_back(makeCandidate(*ittrk, begin, end));
0557 }
0558 }
0559
0560
0561 for (begin = hits.begin(), end = hits.end(); begin != end; ++begin) {
0562 if (*begin)
0563 delete *begin;
0564 }
0565 hits.clear();
0566 }
0567 }
0568
0569
0570
0571 iEvent.put(std::move(output));
0572 }
0573
0574 bool TrackerTrackHitFilter::isFirstValidHitInLayer(const reco::Track &tk, std::vector<bool> isNotValidVec) {
0575 reco::HitPattern hp = tk.hitPattern();
0576
0577 int vecSize = static_cast<int>(isNotValidVec.size());
0578
0579 if (isNotValidVec[vecSize - 1])
0580 return false;
0581
0582
0583 if (vecSize == 1)
0584 return true;
0585
0586 uint32_t pHit = hp.getHitPattern(reco::HitPattern::TRACK_HITS, vecSize - 1);
0587 uint32_t thisLayer = hp.getLayer(pHit);
0588 uint32_t thisSubStruct = hp.getSubStructure(pHit);
0589
0590
0591
0592
0593
0594 for (int j = 0; j < vecSize; ++j) {
0595 if (vecSize > (j + 1)) {
0596 uint32_t nHit = hp.getHitPattern(reco::HitPattern::TRACK_HITS, vecSize - (j + 2));
0597 uint32_t prevLayer = hp.getLayer(nHit);
0598 uint32_t prevSubStruct = hp.getSubStructure(nHit);
0599 if ((thisLayer == prevLayer) && (thisSubStruct == prevSubStruct)) {
0600 if (isNotValidVec[vecSize - (j + 2)] == false) {
0601 return false;
0602 }
0603 } else {
0604 return true;
0605 }
0606 } else {
0607 return true;
0608 }
0609 }
0610 return false;
0611 }
0612
0613 unsigned int TrackerTrackHitFilter::getSequLayer(const reco::Track &tk,
0614 unsigned int prevSequLayers,
0615 std::vector<bool> isNotValidVec) {
0616 unsigned int sequLayers = 0;
0617 sequLayers = isFirstValidHitInLayer(tk, isNotValidVec) ? prevSequLayers + 1 : prevSequLayers;
0618
0619 return sequLayers;
0620 }
0621
0622 TrackCandidate TrackerTrackHitFilter::makeCandidate(const reco::Track &tk,
0623 std::vector<TrackingRecHit *>::iterator hitsBegin,
0624 std::vector<TrackingRecHit *>::iterator hitsEnd) {
0625 PropagationDirection pdir = tk.seedDirection();
0626 PTrajectoryStateOnDet state;
0627 if (pdir == anyDirection)
0628 throw cms::Exception("UnimplementedFeature") << "Cannot work with tracks that have 'anyDirecton' \n";
0629
0630
0631
0632
0633 if ((pdir == alongMomentum) == ((tk.outerPosition() - tk.innerPosition()).Dot(tk.momentum()) >= 0)) {
0634
0635 TrajectoryStateOnSurface originalTsosIn(
0636 trajectoryStateTransform::innerStateOnSurface(tk, *theGeometry, &*theMagField));
0637 state = trajectoryStateTransform::persistentState(originalTsosIn, DetId(tk.innerDetId()));
0638 } else {
0639
0640 TrajectoryStateOnSurface originalTsosOut(
0641 trajectoryStateTransform::outerStateOnSurface(tk, *theGeometry, &*theMagField));
0642 state = trajectoryStateTransform::persistentState(originalTsosOut, DetId(tk.outerDetId()));
0643 }
0644 TrajectorySeed seed(state, TrackCandidate::RecHitContainer(), pdir);
0645 TrackCandidate::RecHitContainer ownHits;
0646 ownHits.reserve(hitsEnd - hitsBegin);
0647 const reco::HitPattern &hp = tk.hitPattern();
0648 unsigned int sequLayers = 0;
0649 std::vector<bool> isNotValidVec;
0650 isNotValidVec.clear();
0651 bool breakHitLoop = false;
0652 if (int(int(hp.numberOfValidHits()) - int(hp.numberOfAllHits(reco::HitPattern::TRACK_HITS))) != 0) {
0653 breakHitLoop = true;
0654 }
0655 for (; hitsBegin != hitsEnd; ++hitsBegin) {
0656
0657 if (truncateTracks_) {
0658 if (breakHitLoop)
0659 break;
0660
0661 if (!(*hitsBegin)->isValid()) {
0662 isNotValidVec.push_back(true);
0663 } else {
0664 isNotValidVec.push_back(false);
0665 }
0666 sequLayers = getSequLayer(tk, sequLayers, isNotValidVec);
0667 if (sequLayers > layersRemaining_)
0668 break;
0669 }
0670 ownHits.push_back(*hitsBegin);
0671 }
0672
0673 TrackCandidate cand(ownHits, seed, state, tk.seedRef());
0674
0675 return cand;
0676 }
0677
0678 void TrackerTrackHitFilter::produceFromTrack(const edm::EventSetup &iSetup,
0679 const Track *itt,
0680 std::vector<TrackingRecHit *> &hits) {
0681
0682 hits.clear();
0683
0684 for (trackingRecHit_iterator ith = itt->recHitsBegin(), edh = itt->recHitsEnd(); ith != edh; ++ith) {
0685 const TrackingRecHit *hit = (*ith);
0686
0687 DetId detid = hit->geographicalId();
0688
0689
0690 if (hit->isValid() && hit == nullptr && detid.rawId() == 0)
0691 continue;
0692
0693 int verdict = checkHit(iSetup, detid, hit);
0694 if (verdict == 0) {
0695
0696 hits.push_back(hit->clone());
0697 } else if (verdict < -2) {
0698
0699 if (replaceWithInactiveHits_) {
0700 hits.push_back(new InvalidTrackingRecHit(*(hit->det()), TrackingRecHit::inactive));
0701 }
0702 } else if (verdict == -2)
0703 hits.push_back(hit->clone());
0704 else if (verdict == -1) {
0705 if (!stripAllInvalidHits_) {
0706 hits.push_back(hit->clone());
0707 }
0708 }
0709 }
0710
0711 }
0712
0713 void TrackerTrackHitFilter::produceFromTrajectory(const edm::EventSetup &iSetup,
0714 const Trajectory *itt,
0715 std::vector<TrackingRecHit *> &hits) {
0716 hits.clear();
0717 nOverlaps = 0;
0718
0719
0720 edm::ESHandle<TrackerTopology> tTopoHand = iSetup.getHandle(tokenTrackerTopo);
0721 const TrackerTopology *tTopo = tTopoHand.product();
0722
0723 std::vector<TrajectoryMeasurement> tmColl = itt->measurements();
0724
0725
0726 const TrajectoryMeasurement *previousTM(nullptr);
0727 DetId previousId(0);
0728
0729
0730
0731 for (std::vector<TrajectoryMeasurement>::const_iterator itTrajMeas = tmColl.begin(); itTrajMeas != tmColl.end();
0732 itTrajMeas++) {
0733 TransientTrackingRecHit::ConstRecHitPointer hitpointer = itTrajMeas->recHit();
0734
0735
0736 if (hitpointer->isValid() && hitpointer->hit() == nullptr) {
0737 continue;
0738 }
0739
0740 const TrackingRecHit *hit = ((*hitpointer).hit());
0741 DetId detid = hit->geographicalId();
0742 int verdict = checkHit(iSetup, detid, hit);
0743
0744 if (verdict == 0) {
0745 if (rejectLowAngleHits_ && !checkHitAngle(*itTrajMeas)) {
0746 verdict = -6;
0747 }
0748 }
0749
0750
0751
0752
0753
0754
0755
0756
0757
0758
0759 if (verdict == 0) {
0760
0761 if (tagOverlaps_) {
0762
0763
0764 int side(sideFromId(detid, tTopo));
0765 int layer(layerFromId(detid, tTopo));
0766 int subDet = detid.subdetId();
0767
0768
0769 if ((previousTM != nullptr) && (layer != -1)) {
0770
0771 for (std::vector<TrajectoryMeasurement>::const_iterator itmCompare = itTrajMeas - 1;
0772 itmCompare >= tmColl.begin() && itmCompare > itTrajMeas - 4;
0773 --itmCompare) {
0774 DetId compareId = itmCompare->recHit()->geographicalId();
0775 if (subDet != compareId.subdetId() || side != sideFromId(compareId, tTopo) ||
0776 layer != layerFromId(compareId, tTopo))
0777 break;
0778 if (!itmCompare->recHit()->isValid())
0779 continue;
0780 if (GeomDetEnumerators::isTrackerPixel(theGeometry->geomDetSubDetector(detid.subdetId())) ||
0781 (GeomDetEnumerators::isTrackerStrip(theGeometry->geomDetSubDetector(detid.subdetId())) &&
0782 SiStripDetId(detid).stereo() ==
0783 SiStripDetId(compareId).stereo())) {
0784
0785
0786
0787
0788
0789
0790
0791
0792
0793
0794 nOverlaps++;
0795 break;
0796 }
0797 }
0798
0799 }
0800
0801 previousTM = &(*itTrajMeas);
0802 previousId = detid;
0803
0804 }
0805
0806
0807 hits.push_back(hit->clone());
0808 }
0809 else if (verdict < -2) {
0810
0811 if (replaceWithInactiveHits_) {
0812 hits.push_back(new InvalidTrackingRecHit(*hit->det(), TrackingRecHit::inactive));
0813 }
0814 } else if (verdict == -2)
0815 hits.push_back(hit->clone());
0816 else if (verdict == -1) {
0817 if (!stripAllInvalidHits_) {
0818 hits.push_back(hit->clone());
0819 }
0820 }
0821 }
0822
0823 std::reverse(hits.begin(), hits.end());
0824 }
0825
0826 int TrackerTrackHitFilter::checkHit(const edm::EventSetup &iSetup, const DetId &detid, const TrackingRecHit *hit) {
0827
0828 edm::ESHandle<TrackerTopology> tTopoHand = iSetup.getHandle(tokenTrackerTopo);
0829 const TrackerTopology *tTopo = tTopoHand.product();
0830
0831 int hitresult = 0;
0832 if (hit->isValid()) {
0833 if (detid.det() == DetId::Tracker) {
0834 bool verdict = true;
0835
0836 for (std::vector<Rule>::const_iterator itr = rules_.begin(), edr = rules_.end(); itr != edr; ++itr) {
0837 itr->apply(detid, tTopo, verdict);
0838 }
0839
0840
0841 if (verdict) {
0842 if (std::binary_search(detsToIgnore_.begin(), detsToIgnore_.end(), detid.rawId())) {
0843 hitresult = -4;
0844 }
0845 } else
0846 hitresult = -3;
0847
0848 if (hitresult == 0 && rejectBadStoNHits_) {
0849 if (!checkStoN(detid, hit))
0850 hitresult = -5;
0851 }
0852 }
0853 else
0854 hitresult = -2;
0855 }
0856 else
0857 hitresult = -1;
0858 return hitresult;
0859 }
0860
0861 bool TrackerTrackHitFilter::checkStoN(const DetId &id, const TrackingRecHit *therechit) {
0862 bool keepthishit = true;
0863
0864
0865
0866
0867 int subdet_cnt = id.subdetId();
0868
0869
0870
0871
0872
0873
0874 if (GeomDetEnumerators::isOuterTracker(theGeometry->geomDetSubDetector(id.subdetId())) &&
0875 !GeomDetEnumerators::isTrackerStrip(theGeometry->geomDetSubDetector(id.subdetId()))) {
0876 return true;
0877 }
0878
0879 if (GeomDetEnumerators::isTrackerStrip(theGeometry->geomDetSubDetector(id.subdetId()))) {
0880 if (subdetStoN_[subdet_cnt - 1]) {
0881
0882 const std::type_info &type = typeid(*therechit);
0883 const SiStripCluster *cluster;
0884 std::optional<const SiStripCluster *> stereo_cluster;
0885 if (type == typeid(SiStripRecHit2D)) {
0886 const SiStripRecHit2D *hit = dynamic_cast<const SiStripRecHit2D *>(therechit);
0887 if (hit != nullptr)
0888 cluster = &*(hit->cluster());
0889 else {
0890 edm::LogError("TrackerTrackHitFilter")
0891 << "TrackerTrackHitFilter::checkStoN : Unknown valid tracker hit in subdet " << id.subdetId()
0892 << "(detID=" << id.rawId() << ")\n ";
0893 keepthishit = false;
0894 }
0895 } else if (type == typeid(SiStripRecHit1D)) {
0896 const SiStripRecHit1D *hit = dynamic_cast<const SiStripRecHit1D *>(therechit);
0897 if (hit != nullptr)
0898 cluster = &*(hit->cluster());
0899 else {
0900 edm::LogError("TrackerTrackHitFilter")
0901 << "TrackerTrackHitFilter::checkStoN : Unknown valid tracker hit in subdet " << id.subdetId()
0902 << "(detID=" << id.rawId() << ")\n ";
0903 keepthishit = false;
0904 }
0905 } else if (type == typeid(SiStripMatchedRecHit2D)) {
0906 const SiStripMatchedRecHit2D *hit = dynamic_cast<const SiStripMatchedRecHit2D *>(therechit);
0907 if (hit != nullptr) {
0908 cluster = &(hit->monoCluster());
0909 stereo_cluster = &(hit->stereoCluster());
0910 } else {
0911 edm::LogError("TrackerTrackHitFilter")
0912 << "TrackerTrackHitFilter::checkStoN : Unknown valid tracker hit in subdet " << id.subdetId()
0913 << "(detID=" << id.rawId() << ")\n ";
0914 keepthishit = false;
0915 }
0916 } else if (type == typeid(ProjectedSiStripRecHit2D)) {
0917 const ProjectedSiStripRecHit2D *hit = dynamic_cast<const ProjectedSiStripRecHit2D *>(therechit);
0918 if (hit != nullptr)
0919 cluster = &*(hit->cluster());
0920 else {
0921 edm::LogError("TrackerTrackHitFilter")
0922 << "TrackerTrackHitFilter::checkStoN : Unknown valid tracker hit in subdet " << id.subdetId()
0923 << "(detID=" << id.rawId() << ")\n ";
0924 keepthishit = false;
0925 }
0926 } else {
0927 throw cms::Exception("Unknown RecHit Type")
0928 << "RecHit of type " << edm::typeDemangle(type.name()) << " not supported.";
0929 }
0930
0931 if (keepthishit) {
0932 double sToNlowCut = subdetStoNlowcut_[subdet_cnt - 1];
0933 double sToNhighCut = subdetStoNhighcut_[subdet_cnt - 1];
0934
0935 float sOverN{-1.f};
0936
0937 if UNLIKELY (stereo_cluster.has_value()) {
0938 const auto &matched = dynamic_cast<const SiStripMatchedRecHit2D *>(therechit);
0939 if (!matched) {
0940 throw cms::Exception("LogicError")
0941 << "expected a SiStripMatchedRecHit2D but could not cast!" << std::endl;
0942 }
0943 siStripClusterInfo_->setCluster(*cluster, matched->monoId());
0944 sOverN = siStripClusterInfo_->signalOverNoise();
0945
0946 siStripClusterInfo_->setCluster(**stereo_cluster, matched->stereoId());
0947 sOverN += siStripClusterInfo_->signalOverNoise();
0948 sOverN *= 0.5f;
0949 } else {
0950 siStripClusterInfo_->setCluster(*cluster, id.rawId());
0951 sOverN = siStripClusterInfo_->signalOverNoise();
0952 }
0953
0954 if ((sToNlowCut > 0) && (sOverN < sToNlowCut)) {
0955 keepthishit = false;
0956 }
0957
0958 if ((sToNhighCut > 0) && (sOverN > sToNhighCut)) {
0959 keepthishit = false;
0960 }
0961
0962 }
0963 }
0964
0965 }
0966 else if (GeomDetEnumerators::isTrackerPixel(theGeometry->geomDetSubDetector(id.subdetId()))) {
0967
0968
0969 keepthishit = true;
0970
0971
0972
0973
0974
0975
0976 if (checkPXLQuality_) {
0977 const SiPixelRecHit *pixelhit = dynamic_cast<const SiPixelRecHit *>(therechit);
0978 if (pixelhit != nullptr) {
0979
0980 float xyprob = pixelhit->clusterProbability(0);
0981
0982 float xychargeprob = pixelhit->clusterProbability(1);
0983
0984 bool haspassed_tplreco = pixelhit->hasFilledProb();
0985 int qbin = pixelhit->qBin();
0986
0987
0988
0989
0990
0991 keepthishit = false;
0992
0993 if (haspassed_tplreco && xyprob > pxlTPLProbXY_ && xychargeprob > pxlTPLProbXYQ_ && qbin > pxlTPLqBin_[0] &&
0994 qbin <= pxlTPLqBin_[1])
0995 keepthishit = true;
0996
0997 } else {
0998 edm::LogInfo("TrackerTrackHitFilter") << "HIT IN PIXEL (" << subdet_cnt << ") but PixelRecHit is EMPTY!!!";
0999 }
1000 }
1001 }
1002
1003
1004
1005
1006 return keepthishit;
1007 }
1008
1009 bool TrackerTrackHitFilter::checkHitAngle(const TrajectoryMeasurement &meas) {
1010 bool angle_ok = false;
1011 bool corrcharge_ok = true;
1012 const TrajectoryStateOnSurface &tsos = meas.updatedState();
1013
1014
1015
1016
1017
1018
1019
1020 if (tsos.isValid()) {
1021
1022 float mom_x = tsos.localDirection().x();
1023 float mom_y = tsos.localDirection().y();
1024 float mom_z = tsos.localDirection().z();
1025
1026 float angle = TMath::ASin(TMath::Abs(mom_z) / sqrt(pow(mom_x, 2) + pow(mom_y, 2) + pow(mom_z, 2)));
1027 if (!rejectLowAngleHits_ || angle >= TrackAngleCut_)
1028 angle_ok = true;
1029
1030
1031 if (angle_ok && PXLcorrClusChargeCut_ > 0.0) {
1032
1033
1034 const TransientTrackingRecHit::ConstRecHitPointer &hitpointer = meas.recHit();
1035 if (hitpointer->isValid()) {
1036 const TrackingRecHit *hit = (*hitpointer).hit();
1037 if (GeomDetEnumerators::isTrackerPixel(
1038 theGeometry->geomDetSubDetector(hit->geographicalId().subdetId()))) {
1039 corrcharge_ok = false;
1040 float clust_alpha = atan2(mom_z, mom_x);
1041 float clust_beta = atan2(mom_z, mom_y);
1042
1043
1044
1045 const SiPixelRecHit *pixelhit = dynamic_cast<const SiPixelRecHit *>(hit);
1046 float clust_charge = pixelhit->cluster()->charge();
1047 float corr_clust_charge =
1048 clust_charge * sqrt(1.0 / (1.0 / pow(tan(clust_alpha), 2) + 1.0 / pow(tan(clust_beta), 2) + 1.0));
1049
1050 if (corr_clust_charge > PXLcorrClusChargeCut_) {
1051 corrcharge_ok = true;
1052 }
1053 }
1054 }
1055
1056 }
1057
1058 }
1059 else {
1060 edm::LogWarning("TrackerTrackHitFilter") << "TSOS not valid ! Impossible to calculate track angle.";
1061 }
1062
1063 return angle_ok && corrcharge_ok;
1064 }
1065
1066 bool TrackerTrackHitFilter::checkPXLCorrClustCharge(const TrajectoryMeasurement &meas) {
1067
1068
1069
1070
1071 bool corrcharge_ok = false;
1072
1073 const TransientTrackingRecHit::ConstRecHitPointer &hitpointer = meas.recHit();
1074 if (!hitpointer->isValid())
1075 return corrcharge_ok;
1076 const TrackingRecHit *hit = (*hitpointer).hit();
1077 if (GeomDetEnumerators::isTrackerStrip(
1078 theGeometry->geomDetSubDetector(hit->geographicalId().subdetId()))) {
1079 return corrcharge_ok;
1080 }
1081
1082 const TrajectoryStateOnSurface &tsos = meas.updatedState();
1083 if (tsos.isValid()) {
1084 float mom_x = tsos.localDirection().x();
1085 float mom_y = tsos.localDirection().y();
1086 float mom_z = tsos.localDirection().z();
1087 float clust_alpha = atan2(mom_z, mom_x);
1088 float clust_beta = atan2(mom_z, mom_y);
1089
1090
1091
1092 const SiPixelRecHit *pixelhit = dynamic_cast<const SiPixelRecHit *>(hit);
1093 float clust_charge = pixelhit->cluster()->charge();
1094 float corr_clust_charge =
1095 clust_charge * sqrt(1.0 / (1.0 / pow(tan(clust_alpha), 2) + 1.0 / pow(tan(clust_beta), 2) + 1.0));
1096 if (corr_clust_charge > PXLcorrClusChargeCut_)
1097 corrcharge_ok = true;
1098
1099 }
1100 return corrcharge_ok;
1101
1102 }
1103
1104 int TrackerTrackHitFilter::layerFromId(const DetId &id, const TrackerTopology *tTopo) const {
1105 return tTopo->layer(id);
1106 }
1107
1108 int TrackerTrackHitFilter::sideFromId(const DetId &id, const TrackerTopology *tTopo) const {
1109 return tTopo->side(id);
1110 }
1111
1112
1113 void TrackerTrackHitFilter::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
1114 edm::ParameterSetDescription desc;
1115 desc.setComment("");
1116 desc.add<edm::InputTag>("src", edm::InputTag("generalTracks"));
1117 desc.add<uint32_t>("minimumHits", 3)->setComment("number of hits for refit");
1118 desc.add<bool>("replaceWithInactiveHits", false)
1119 ->setComment(
1120 " instead of removing hits replace them with inactive hits, so you still consider the multiple "
1121 "scattering");
1122 desc.add<bool>("truncateTracks", false)
1123 ->setComment(
1124 "determines if recHits collection is to be truncated to provide tracks with layersRemaining number of "
1125 "layers after refitting");
1126 desc.add<uint32_t>("layersRemaining", 8)
1127 ->setComment("number of tracker layers with measurement remaining after truncating track");
1128 desc.add<bool>("stripFrontInvalidHits", false)
1129 ->setComment("strip invalid & inactive hits from any end of the track");
1130 desc.add<bool>("stripBackInvalidHits", false)
1131 ->setComment("strip invalid & inactive hits from any end of the track");
1132 desc.add<bool>("stripAllInvalidHits", false)->setComment("dangerous to turn on, you might forget about MS");
1133 desc.add<bool>("isPhase2", false);
1134 desc.add<bool>("rejectBadStoNHits", false);
1135 desc.add<std::string>("CMNSubtractionMode", std::string("Median"))->setComment("TT6");
1136 desc.add<std::vector<uint32_t> >("detsToIgnore", {});
1137 desc.add<bool>("useTrajectories", false);
1138 desc.add<bool>("rejectLowAngleHits", false);
1139 desc.add<double>("TrackAngleCut", 0.25)->setComment("rad");
1140 desc.add<bool>("usePixelQualityFlag", false);
1141 desc.add<double>("PxlTemplateProbXYCut", 0.000125);
1142 desc.add<double>("PxlTemplateProbXYChargeCut", -99.);
1143 desc.add<std::vector<int32_t> >("PxlTemplateqBinCut", {0, 3});
1144 desc.add<double>("PxlCorrClusterChargeCut", -999.0);
1145 desc.add<bool>("tagOverlaps", false);
1146 desc.add<std::vector<std::string> >("commands", {})->setComment("layers to remove");
1147 desc.add<std::vector<std::string> >("StoNcommands", {})->setComment("S/N cut per layer");
1148 descriptions.addWithDefaultLabel(desc);
1149 }
1150
1151 }
1152 }
1153
1154
1155 #include "FWCore/PluginManager/interface/ModuleDef.h"
1156 #include "FWCore/Framework/interface/MakerMacros.h"
1157
1158 using reco::modules::TrackerTrackHitFilter;
1159 DEFINE_FWK_MODULE(TrackerTrackHitFilter);