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