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