FinalTrackSelectors/plugins/TrackerTrackHitFilter.cc

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 //for S/N cut
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 //for angle cut
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 //#include <math>
0056
0057 /**
0058  * Configurables:
0059  *
0060  *   Generic:
0061  *     tracks                  = InputTag of a collection of tracks to read from
0062  *     minimumHits             = Minimum hits that the output TrackCandidate must have to be saved
0063  *     replaceWithInactiveHits = instead of discarding hits, replace them with a invalid "inactive" hits,
0064  *                               so multiple scattering is accounted for correctly.
0065  *     truncateTracks          = determines if recHits collection is to be truncated to provide tracks with
0066  *                               layersRemaining number of layers after refitting
0067  *     layersRemaining         = number of tracker layers with measurement remaining after truncating track
0068  *     stripFrontInvalidHits   = strip invalid hits at the beginning of the track
0069  *     stripBackInvalidHits    = strip invalid hits at the end of the track
0070  *     stripAllInvalidHits     = remove ALL invald hits (might be a problem for multiple scattering, use with care!)
0071  *
0072  *   Per structure:
0073  *      commands = list of commands, to be applied in order as they are written
0074  *      commands can be:
0075  *          "keep XYZ"  , "drop XYZ"    (XYZ = PXB, PXE, TIB, TID, TOB, TEC)
0076  *          "keep XYZ l", "drop XYZ n"  (XYZ as above, n = layer, wheel or disk = 1 .. 6 ; positive and negative are the same )
0077  *
0078  *   Individual modules:
0079  *     detsToIgnore        = individual list of detids on which hits must be discarded
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         // parse a rule from a string
0102         Rule(const std::string &str);
0103         // check this DetId, update the verdict if the rule has anything to say about it
0104         void apply(DetId detid, const TrackerTopology *tTopo, bool &verdict) const {
0105           // check detector
0106           if (detid.subdetId() == subdet_) {
0107             // check layer
0108             if ((layer_ == 0) || (layer_ == layer(detid, tTopo))) {
0109               // override verdict
0110               verdict = keep_;
0111               //  std::cout<<"Verdict is "<<verdict<<" for subdet "<<subdet_<<", layer "<<layer_<<"! "<<std::endl;
0112             }
0113             // else std::cout<<"No, sorry, wrong layer.Retry..."<<std::endl;
0114           }
0115           //    else{ std::cout<<"No, sorry, wrong subdet.Retry..."<<std::endl;}
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_;           //(6); //,std::bool(false));
0145       std::vector<double> subdetStoNlowcut_;   //(6,-1.0);
0146       std::vector<double> subdetStoNhighcut_;  //(6,-1.0);
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       // bool checkOverlapHit();
0181
0182       TrackCandidate makeCandidate(const reco::Track &tk,
0183                                    std::vector<TrackingRecHit *>::iterator hitsBegin,
0184                                    std::vector<TrackingRecHit *>::iterator hitsEnd);
0185       //const TransientTrackingRecHitBuilder *RHBuilder;
0186     };  // class
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       // match and check it works
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       // Set up fields:
0201       //  rule type
0202       keep_ = (strncmp(match[1].first, "keep", 4) == 0);
0203
0204       //  subdet
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       //   layer (if present)
0223       if (match[4].first != match[4].second) {
0224         layer_ = atoi(match[4].first);
0225       } else {
0226         layer_ = 0;
0227       }
0228     }  //end Rule::Rule
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       // match a set of capital case chars (preceded by an arbitrary number of leading blanks),
0236       //followed b an arbitrary number of blanks, one or more digits (not necessary, they cannot also be,
0237       // another set of blank spaces and, again another *eventual* digit
0238       // static boost::regex rule("\\s+([A-Z]+)(\\s+(\\d+)(\\.)?(\\d+))?(\\s+(\\d+)(\\.)?(\\d+))?");
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       // match and check it works
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++) {  //loop on subdets
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++) {  //check that the regex was correct
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       //std::cout<<"Reached end of parseStoN"<<std::endl;
0307     }  //end parseStoN
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       // construct the SiStripClusterInfo object only for Phase-0 / Phase-1
0332       // no Strip modules in Phase-2
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       // sanity check
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       // read and parse commands
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_) {  //commands for S/N cut
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         ////edm::LogDebug("TrackerTrackHitFilter")
0390         edm::LogInfo("TrackerTrackHitFilter") << "Finished parsing S/N. Applying following cuts to subdets:";
0391         for (cnt = 0; cnt < 6; cnt++) {
0392           ////edm::LogDebug("TrackerTrackHitFilter")
0393           edm::LogVerbatim("TrackerTrackHitFilter")
0394               << "Subdet #" << cnt + 1 << " -> " << subdetStoNlowcut_[cnt] << " , " << subdetStoNhighcut_[cnt];
0395         }
0396       }  //end if rejectBadStoNHits_
0397
0398       if (rejectLowAngleHits_)
0399         edm::LogInfo("TrackerTrackHitFilter") << "\nApplying cut on angle track = " << TrackAngleCut_;
0400
0401       // sort detids to ignore
0402       std::sort(detsToIgnore_.begin(), detsToIgnore_.end());
0403
0404       // issue the produce<>
0405       produces<TrackCandidateCollection>();
0406     }
0407
0408     void TrackerTrackHitFilter::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0409       //Dump Run and Event
0410       iRun = iEvent.id().run();
0411       iEvt = iEvent.id().event();
0412
0413       // read with View, so we can read also a TrackRefVector
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       // read from EventSetup
0423       theGeometry = iSetup.getHandle(tokenGeometry);
0424       theMagField = iSetup.getHandle(tokenMagField);
0425       if (!isPhase2_)
0426         siStripClusterInfo_->initEvent(iSetup);
0427
0428       // prepare output collection
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       // working area and tools
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;  //trajectory in the collection
0445           const reco::TrackRef tkref = itass->val;                     //associated track track in the collection
0446           //std::cout<<"The hit collection has size "<<hits.size()<<" (should be 0) while the track contains initially "<< tkref->recHitsEnd() - tkref->recHitsBegin()<<std::endl;
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           // strip invalid hits at the beginning
0455           if (stripFrontInvalidHits_) {
0456             while ((begin != end) && ((*begin)->isValid() == false))
0457               ++begin;
0458           }
0459           // strip invalid hits at the end
0460           if (stripBackInvalidHits_ && (begin != end)) {
0461             --end;
0462             while ((begin != end) && ((*end)->isValid() == false))
0463               --end;
0464             ++end;
0465           }
0466
0467           // if we still have some hits build the track candidate
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 {  //all invalid hits have been already kicked out
0478             if ((end - begin) >= int(minimumHits_)) {
0479               output->push_back(makeCandidate(*trk, begin, end));
0480             }
0481           }
0482
0483           // if we still have some hits build the track candidate
0484           //if ((end - begin) >= int(minimumHits_)) {
0485           //    output->push_back( makeCandidate ( *trk, begin, end ) );
0486           //}
0487
0488           // now delete the hits not used by the candidate
0489           for (begin = hits.begin(), end = hits.end(); begin != end; ++begin) {
0490             if (*begin)
0491               delete *begin;
0492           }
0493           hits.clear();
0494         }  // loop on trajectories
0495
0496       } else {  //use plain tracks
0497
0498         // loop on tracks
0499         for (std::vector<reco::Track>::const_iterator ittrk = tracks->begin(), edtrk = tracks->end(); ittrk != edtrk;
0500              ++ittrk) {
0501           //    std::cout<<"The hit collection has size "<<hits.size()<<" (should be 0) while the track contains initially "<< ittrk->recHitsEnd() - ittrk->recHitsBegin()<<std::endl;
0502
0503           const Track *trk = &(*ittrk);
0504
0505           produceFromTrack(iSetup, trk, hits);
0506           //-----------------------
0507           /*
0508       std::cout<<"Hit collection in output has size "<<hits.size()<<". Dumping hit positions..."<<std::endl;
0509         for (std::vector<TrackingRecHit *>::iterator ith = hits.begin(), edh = hits.end(); ith != edh; ++ith) {
0510       const TrackingRecHit *myhit = *(ith);
0511         TransientTrackingRecHit::RecHitPointer ttrh;
0512         float radius=0.0;
0513         float xx=-999.0,yy=-999.0,zz=-999.0;
0514         unsigned int myid=0;
0515         if(myhit->isValid()){
0516           ttrh = RHBuilder->build(myhit);
0517           xx=ttrh->globalPosition().x();
0518           yy=ttrh->globalPosition().y();
0519           zz=ttrh->globalPosition().z();
0520           radius = sqrt( pow(xx,2)+pow(yy,2) );
0521           myid=myhit->geographicalId().rawId();
0522         }
0523         std::cout<<"-$-$ OUTPUT Hit position: ( "<<xx<<" , " <<yy<<" , " <<zz<<" ) , RADIUS = "  <<radius<<"  on DetID= "<< myid<<std::endl;
0524     }//end loop on hits
0525       */
0526           //-----------------------
0527
0528           std::vector<TrackingRecHit *>::iterator begin = hits.begin(), end = hits.end();
0529           // std::cout << "Back in the main producer (TRK), the final hit collection has size " << hits.size() << std::endl;
0530           // strip invalid hits at the beginning
0531           if (stripFrontInvalidHits_) {
0532             while ((begin != end) && ((*begin)->isValid() == false))
0533               ++begin;
0534           }
0535           // strip invalid hits at the end
0536           if (stripBackInvalidHits_ && (begin != end)) {
0537             --end;
0538             while ((begin != end) && ((*end)->isValid() == false))
0539               --end;
0540             ++end;
0541           }
0542
0543           // if we still have some hits build the track candidate
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 {  //all invalid hits have been already kicked out
0555             if ((end - begin) >= int(minimumHits_)) {
0556               output->push_back(makeCandidate(*ittrk, begin, end));
0557             }
0558           }
0559
0560           // now delete the hits not used by the candidate
0561           for (begin = hits.begin(), end = hits.end(); begin != end; ++begin) {
0562             if (*begin)
0563               delete *begin;
0564           }
0565           hits.clear();
0566         }  // loop on tracks
0567       }  //end else useTracks
0568
0569       // std::cout<<"OUTPUT SIZE: "<<output->size()<<std::endl;
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       // If hit is not valid, it will not count as a tracker layer with measurement -> don't increase sequLayers
0579       if (isNotValidVec[vecSize - 1])
0580         return false;
0581
0582       // If very first valid layer -> increase sequLayers
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       // This loop compares the previous hits substructure and layer with current hit. If hits in the same layer
0591       // and substructure and previous hit is valid, skip layer. If previous hit is not valid, even if in same layer
0592       // and substructure, check the previous previous hit. Repeat process for every previous hit until reaching
0593       // a valid hit or different layer/substructure
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       //  double innerP=sqrt( pow(tk.innerMomentum().X(),2)+pow(tk.innerMomentum().Y(),2)+pow(tk.innerMomentum().Z(),2) );
0631       //  if ( (pdir == alongMomentum) == ( innerP >= tk.outerP() ) ) {
0632
0633       if ((pdir == alongMomentum) == ((tk.outerPosition() - tk.innerPosition()).Dot(tk.momentum()) >= 0)) {
0634         // use inner state
0635         TrajectoryStateOnSurface originalTsosIn(
0636             trajectoryStateTransform::innerStateOnSurface(tk, *theGeometry, &*theMagField));
0637         state = trajectoryStateTransform::persistentState(originalTsosIn, DetId(tk.innerDetId()));
0638       } else {
0639         // use outer state
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         // Only perform hit checking and verification of number of layers if recHits are to be truncated
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       // loop on tracks
0682       hits.clear();  // extra safety
0683
0684       for (trackingRecHit_iterator ith = itt->recHitsBegin(), edh = itt->recHitsEnd(); ith != edh; ++ith) {
0685         const TrackingRecHit *hit = (*ith);  // ith is an iterator on edm::Ref to rechit
0686
0687         DetId detid = hit->geographicalId();
0688
0689         //check that the hit is a real hit and not a constraint
0690         if (hit->isValid() && hit == nullptr && detid.rawId() == 0)
0691           continue;
0692
0693         int verdict = checkHit(iSetup, detid, hit);
0694         if (verdict == 0) {
0695           // just copy the hit
0696           hits.push_back(hit->clone());
0697         } else if (verdict < -2) {  //hit rejected because did not pass the selections
0698                                     // still, if replaceWithInactiveHits is true we have to put a new hit
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());  //hit not in the tracker
0704         else if (verdict == -1) {        //hit not valid
0705           if (!stripAllInvalidHits_) {
0706             hits.push_back(hit->clone());
0707           }
0708         }
0709       }  // loop on hits
0710
0711     }  //end TrackerTrackHitFilter::produceFromTrack
0712
0713     void TrackerTrackHitFilter::produceFromTrajectory(const edm::EventSetup &iSetup,
0714                                                       const Trajectory *itt,
0715                                                       std::vector<TrackingRecHit *> &hits) {
0716       hits.clear();  // extra safety
0717       nOverlaps = 0;
0718
0719       //Retrieve tracker topology from geometry
0720       edm::ESHandle<TrackerTopology> tTopoHand = iSetup.getHandle(tokenTrackerTopo);
0721       const TrackerTopology *tTopo = tTopoHand.product();
0722
0723       std::vector<TrajectoryMeasurement> tmColl = itt->measurements();
0724
0725       //---OverlapBegin needed eventually for overlaps, but I must create them here in any case
0726       const TrajectoryMeasurement *previousTM(nullptr);
0727       DetId previousId(0);
0728       //int previousLayer(-1);
0729       ///---OverlapEnd
0730
0731       for (std::vector<TrajectoryMeasurement>::const_iterator itTrajMeas = tmColl.begin(); itTrajMeas != tmColl.end();
0732            itTrajMeas++) {
0733         TransientTrackingRecHit::ConstRecHitPointer hitpointer = itTrajMeas->recHit();
0734
0735         //check that the hit is a real hit and not a constraint
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)) {  //check angle of track on module if requested
0746             verdict = -6;                                            //override previous verdicts
0747           }
0748         }
0749
0750         /*
0751     //this has been included in checkHitAngle(*itTrajMeas)
0752     if (verdict == 0) {
0753     if( PXLcorrClusChargeCut_>0.0  && !checkPXLCorrClustCharge(*itTrajMeas) ){//check angle of track on module if requested
0754     verdict=-7;//override previous verdicts
0755     }
0756     }
0757     */
0758
0759         if (verdict == 0) {  // Hit TAKEN !!!!
0760
0761           if (tagOverlaps_) {  ///---OverlapBegin
0762             //std::cout<<"Looking for overlaps in Run="<<iRun<<" , Event ="<<iEvt<<std::flush;
0763
0764             int side(sideFromId(detid, tTopo));    //side 0=barrel, 1=minus , 2=plus
0765             int layer(layerFromId(detid, tTopo));  //layer or disk
0766             int subDet = detid.subdetId();
0767             //std::cout  << "  Check Subdet #" <<subDet << ", layer = " <<layer<<" stereo: "<< ((subDet > 2)?(SiStripDetId(detid).stereo()):2);
0768
0769             if ((previousTM != nullptr) && (layer != -1)) {
0770               //std::cout<<"A previous TM exists! "<<std::endl;
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())) {  //if either pixel or strip stereo module
0784                   //  overlapHits.push_back(std::make_pair(&(*itmCompare),&(*itm)));
0785                   //std::cout<< "Adding pair "<< ((subDet >2)?(SiStripDetId(detid).stereo()):2)
0786                   //     << " from SubDet = "<<subDet<<" , layer = " << layer<<"  Run:"<<iRun<<"\tEv: "<<iEvt<<"\tId1: "<<compareId.rawId()<<"\tId2: "<<detid.rawId()<<std::endl;
0787                   // if(abs(compareId.rawId()-detid.rawId())==1)std::cout<<"These two are from the same det! Id1= "<<detid.rawId()<<" has stereo type "<<SiStripDetId(detid).stereo() <<"\tId2: "<<compareId.rawId()<<" has stereo type "<<SiStripDetId(compareId).stereo()<<std::endl;
0788                   ///
0789                   // if(detid.rawId()<compareId.rawId()){
0790                   // std::cout<< "+++ "<< "\t"<<iRun<<"\t"<<iEvt<<"\t"<<detid.rawId()<<"\t"<<compareId.rawId()<<std::endl;
0791                   // }
0792                   //else  std::cout<< "+++ "<< "\t"<<iRun<<"\t"<<iEvt<<"\t"<<compareId.rawId()<<"\t"<<detid.rawId()<<std::endl;
0793
0794                   nOverlaps++;
0795                   break;
0796                 }
0797               }  //end second loop on TM for overlap tagging
0798
0799             }  //end   if ( (layer!=-1 )&&(acceptLayer[subDet]))
0800
0801             previousTM = &(*itTrajMeas);
0802             previousId = detid;
0803             //previousLayer = layer;
0804           }  //end if look for overlaps
0805           ///---OverlapEnd
0806
0807           hits.push_back(hit->clone());  //just copy it
0808         }  //end if HIT TAKEN
0809         else if (verdict < -2) {  //hit rejected because did not pass the selections
0810           // still, if replaceWithInactiveHits is true we have to put a new hit
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());  //hit not in the tracker
0816         else if (verdict == -1) {        //hit not valid
0817           if (!stripAllInvalidHits_) {
0818             hits.push_back(hit->clone());
0819           }
0820         }
0821       }  // loop on hits
0822
0823       std::reverse(hits.begin(), hits.end());
0824     }  //end TrackerTrackHitFilter::produceFromTrajectories
0825
0826     int TrackerTrackHitFilter::checkHit(const edm::EventSetup &iSetup, const DetId &detid, const TrackingRecHit *hit) {
0827       //Retrieve tracker topology from geometry
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) {  // check for tracker hits
0834           bool verdict = true;
0835           // first check at structure level
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           // if the hit is good, check again at module level
0841           if (verdict) {
0842             if (std::binary_search(detsToIgnore_.begin(), detsToIgnore_.end(), detid.rawId())) {
0843               hitresult = -4;
0844             }
0845           } else
0846             hitresult = -3;
0847           //if the hit is in the desired part of the det, check other things
0848           if (hitresult == 0 && rejectBadStoNHits_) {
0849             if (!checkStoN(detid, hit))
0850               hitresult = -5;
0851           }  //end if S/N is ok
0852         }  //end hit in tracker
0853         else
0854           hitresult = -2;
0855       }  //end hit is valid
0856       else
0857         hitresult = -1;  //invalid hit
0858       return hitresult;
0859     }  //end  TrackerTrackHitFilter::checkHit()
0860
0861     bool TrackerTrackHitFilter::checkStoN(const DetId &id, const TrackingRecHit *therechit) {
0862       bool keepthishit = true;
0863       // const uint32_t& recHitDetId = id.rawId();
0864
0865       //check StoN only if subdet was set by the user
0866       //  int subdet_cnt=0;
0867       int subdet_cnt = id.subdetId();
0868
0869       //  for(subdet_cnt=0;subdet_cnt<6; ++subdet_cnt){
0870
0871       //  if( subdetStoN_[subdet_cnt-1]&& (id.subdetId()==subdet_cnt)  ){//check that hit is in a det belonging to a subdet where we decided to apply a S/N cut
0872
0873       // for phase-2 OT placehold, do nothing
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           //check that hit is in a det belonging to a subdet where we decided to apply a S/N cut
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               // if it's matched rechit, use the average of stereo and mono clusters
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             //if(!keepthishit)std::cout<<"Hit rejected because of bad S/N: "<<siStripClusterInfo_->signalOverNoise()<<std::endl;
0962           }
0963         }  //end if  subdetStoN_[subdet_cnt]&&...
0964
0965       }  //end if subdet_cnt >2
0966       else if (GeomDetEnumerators::isTrackerPixel(theGeometry->geomDetSubDetector(id.subdetId()))) {  //pixel
0967         //pixels have naturally a very low noise (because of their low capacitance). So the S/N cut is
0968         //irrelevant in this case. Leave it dummy
0969         keepthishit = true;
0970
0971         /**************
0972        * Cut on cluster charge corr by angle embedded in the checkHitAngle() function
0973        *
0974        *************/
0975
0976         if (checkPXLQuality_) {
0977           const SiPixelRecHit *pixelhit = dynamic_cast<const SiPixelRecHit *>(therechit);
0978           if (pixelhit != nullptr) {
0979             //std::cout << "ClusterCharge=" <<std::flush<<pixelhit->cluster()->charge() << std::flush;
0980             float xyprob = pixelhit->clusterProbability(0);        //x-y combined log_e probability of the pixel cluster
0981                                                                    //singl x- and y-prob not stored sicne CMSSW 3_9_0
0982             float xychargeprob = pixelhit->clusterProbability(1);  //xy-prob * charge prob
0983             //  float chargeprob   =pixelhit->clusterProbability(2);//charge prob
0984             bool haspassed_tplreco = pixelhit->hasFilledProb();  //the cluster was associted to a template
0985             int qbin = pixelhit->qBin();  //R==meas_charge/pred_charge:  Qbin=0 ->R>1.5 , =1->1<R<1.5 ,=2->0.85<R<1 ,
0986                                           // Qbin=3->0.95*Qminpred<R<0.85 ,=4->, =5->meas_charge<0.95*Qminpred
0987
0988             //  if(haspassed_tplreco)   std::cout<<"  CLUSTPROB=\t"<<xprob<<"\t"<<yprob<<"\t"<<combprob<<"\t"<<qbin<<std::endl;
0989             //  else std::cout<<"CLUSTPROBNOTDEF=\t"<<xprob<<"\t"<<yprob<<"\t"<<combprob<<"\t"<<qbin<<std::endl;
0990
0991             keepthishit = false;
0992             //  std::cout<<"yyyyy "<<qbin<<" "<<xprob<<"  "<<yprob<<std::endl;
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         }  //end if check pixel quality flag
1001       }
1002       //    else  throw cms::Exception("TrackerTrackHitFilter") <<"Loop over subdetector out of range when applying the S/N cut: "<<subdet_cnt;
1003
1004       //  }//end loop on subdets
1005
1006       return keepthishit;
1007     }  //end CheckStoN
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   edm::LogDebug("TrackerTrackHitFilter")<<"TSOS parameters: ";
1015   edm::LogDebug("TrackerTrackHitFilter") <<"Global momentum: "<<tsos.globalMomentum().x()<<"  "<<tsos.globalMomentum().y()<<"  "<<tsos.globalMomentum().z();
1016   edm::LogDebug("TrackerTrackHitFilter") <<"Local momentum: "<<tsos.localMomentum().x()<<"  "<<tsos.localMomentum().y()<<"  "<<tsos.localMomentum().z();
1017   edm::LogDebug("TrackerTrackHitFilter") <<"Track charge: "  <<tsos.charge();
1018   edm::LogDebug("TrackerTrackHitFilter")<<"Local position: "  <<tsos.localPosition().x()<<"  "<<tsos.localPosition().y()<<"  "<<tsos.localPosition().z();
1019   */
1020       if (tsos.isValid()) {
1021         //check the angle of this tsos
1022         float mom_x = tsos.localDirection().x();
1023         float mom_y = tsos.localDirection().y();
1024         float mom_z = tsos.localDirection().z();
1025         //we took LOCAL momentum, i.e. respect to surface. Thus the plane is z=0
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;  // keep this hit
1029         // else  std::cout<<"Hit rejected because angle is "<< angle<<" ( <"<<TrackAngleCut_<<" )"<<std::endl;
1030
1031         if (angle_ok && PXLcorrClusChargeCut_ > 0.0) {
1032           //
1033           //get the hit from the TM and check that it is in the pixel
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()))) {  //do it only for pixel hits
1039               corrcharge_ok = false;
1040               float clust_alpha = atan2(mom_z, mom_x);
1041               float clust_beta = atan2(mom_z, mom_y);
1042
1043               //Now get the cluster charge
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               //std::cout<<"xxxxx "<<clust_charge<<" "<<corr_clust_charge<<"  " <<pixelhit->qBin()<<"  "<<pixelhit->clusterProbability(1)<<"  "<<pixelhit->clusterProbability(2)<< std::endl;
1050               if (corr_clust_charge > PXLcorrClusChargeCut_) {
1051                 corrcharge_ok = true;
1052               }
1053             }  //end if hit is in pixel
1054           }  //end if hit is valid
1055
1056         }  //check corr cluster charge for pixel hits
1057
1058       }  //end if TSOS is valid
1059       else {
1060         edm::LogWarning("TrackerTrackHitFilter") << "TSOS not valid ! Impossible to calculate track angle.";
1061       }
1062
1063       return angle_ok && corrcharge_ok;
1064     }  //end TrackerTrackHitFilter::checkHitAngle
1065
1066     bool TrackerTrackHitFilter::checkPXLCorrClustCharge(const TrajectoryMeasurement &meas) {
1067       /*
1068     Code taken from DPGAnalysis/SiPixelTools/plugins/PixelNtuplizer_RealData.cc
1069   */
1070
1071       bool corrcharge_ok = false;
1072       //get the hit from the TM and check that it is in the pixel
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()))) {  //SiStrip hit, skip
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         //Now get the cluster charge
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       }  //end if TSOS is valid
1100       return corrcharge_ok;
1101
1102     }  //end TrackerTrackHitFilter::checkPXLCorrClustCharge
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     // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
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   }  // namespace modules
1152 }  // namespace reco
1153
1154 // ========= MODULE DEF ==============
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);