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  *     stripFrontInvalidHits   = strip invalid hits at the beginning of the track
0066  *     stripBackInvalidHits    = strip invalid hits at the end of the track
0067  *     stripAllInvalidHits     = remove ALL invald hits (might be a problem for multiple scattering, use with care!)
0068  *
0069  *   Per structure:
0070  *      commands = list of commands, to be applied in order as they are written
0071  *      commands can be:
0072  *          "keep XYZ"  , "drop XYZ"    (XYZ = PXB, PXE, TIB, TID, TOB, TEC)
0073  *          "keep XYZ l", "drop XYZ n"  (XYZ as above, n = layer, wheel or disk = 1 .. 6 ; positive and negative are the same )
0074  *
0075  *   Individual modules:
0076  *     detsToIgnore        = individual list of detids on which hits must be discarded
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         // parse a rule from a string
0097         Rule(const std::string &str);
0098         // check this DetId, update the verdict if the rule has anything to say about it
0099         void apply(DetId detid, const TrackerTopology *tTopo, bool &verdict) const {
0100           // check detector
0101           if (detid.subdetId() == subdet_) {
0102             // check layer
0103             if ((layer_ == 0) || (layer_ == layer(detid, tTopo))) {
0104               // override verdict
0105               verdict = keep_;
0106               //  std::cout<<"Verdict is "<<verdict<<" for subdet "<<subdet_<<", layer "<<layer_<<"! "<<std::endl;
0107             }
0108             // else std::cout<<"No, sorry, wrong layer.Retry..."<<std::endl;
0109           }
0110           //    else{ std::cout<<"No, sorry, wrong subdet.Retry..."<<std::endl;}
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_;           //(6); //,std::bool(false));
0137       std::vector<double> subdetStoNlowcut_;   //(6,-1.0);
0138       std::vector<double> subdetStoNhighcut_;  //(6,-1.0);
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       // bool checkOverlapHit();
0173
0174       TrackCandidate makeCandidate(const reco::Track &tk,
0175                                    std::vector<TrackingRecHit *>::iterator hitsBegin,
0176                                    std::vector<TrackingRecHit *>::iterator hitsEnd);
0177       //const TransientTrackingRecHitBuilder *RHBuilder;
0178     };  // class
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       // match and check it works
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       // Set up fields:
0193       //  rule type
0194       keep_ = (strncmp(match[1].first, "keep", 4) == 0);
0195
0196       //  subdet
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       //   layer (if present)
0215       if (match[4].first != match[4].second) {
0216         layer_ = atoi(match[4].first);
0217       } else {
0218         layer_ = 0;
0219       }
0220     }  //end Rule::Rule
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       // match a set of capital case chars (preceded by an arbitrary number of leading blanks),
0228       //followed b an arbitrary number of blanks, one or more digits (not necessary, they cannot also be,
0229       // another set of blank spaces and, again another *eventual* digit
0230       // static boost::regex rule("\\s+([A-Z]+)(\\s+(\\d+)(\\.)?(\\d+))?(\\s+(\\d+)(\\.)?(\\d+))?");
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       // match and check it works
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++) {  //loop on subdets
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++) {  //check that the regex was correct
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       //std::cout<<"Reached end of parseStoN"<<std::endl;
0299     }  //end parseStoN
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       // construct the SiStripClusterInfo object only for Phase-0 / Phase-1
0322       // no Strip modules in Phase-2
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       // sanity check
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       // read and parse commands
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_) {  //commands for S/N cut
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         ////edm::LogDebug("TrackerTrackHitFilter")
0380         edm::LogInfo("TrackerTrackHitFilter") << "Finished parsing S/N. Applying following cuts to subdets:";
0381         for (cnt = 0; cnt < 6; cnt++) {
0382           ////edm::LogDebug("TrackerTrackHitFilter")
0383           edm::LogVerbatim("TrackerTrackHitFilter")
0384               << "Subdet #" << cnt + 1 << " -> " << subdetStoNlowcut_[cnt] << " , " << subdetStoNhighcut_[cnt];
0385         }
0386       }  //end if rejectBadStoNHits_
0387
0388       if (rejectLowAngleHits_)
0389         edm::LogInfo("TrackerTrackHitFilter") << "\nApplying cut on angle track = " << TrackAngleCut_;
0390
0391       // sort detids to ignore
0392       std::sort(detsToIgnore_.begin(), detsToIgnore_.end());
0393
0394       // issue the produce<>
0395       produces<TrackCandidateCollection>();
0396     }
0397
0398     void TrackerTrackHitFilter::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0399       //Dump Run and Event
0400       iRun = iEvent.id().run();
0401       iEvt = iEvent.id().event();
0402
0403       // read with View, so we can read also a TrackRefVector
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       // read from EventSetup
0413       theGeometry = iSetup.getHandle(tokenGeometry);
0414       theMagField = iSetup.getHandle(tokenMagField);
0415       if (!isPhase2_)
0416         siStripClusterInfo_->initEvent(iSetup);
0417
0418       // prepare output collection
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       // working area and tools
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;  //trajectory in the collection
0435           const reco::TrackRef tkref = itass->val;                     //associated track track in the collection
0436           //std::cout<<"The hit collection has size "<<hits.size()<<" (should be 0) while the track contains initially "<< tkref->recHitsEnd() - tkref->recHitsBegin()<<std::endl;
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           // strip invalid hits at the beginning
0445           if (stripFrontInvalidHits_) {
0446             while ((begin != end) && ((*begin)->isValid() == false))
0447               ++begin;
0448           }
0449           // strip invalid hits at the end
0450           if (stripBackInvalidHits_ && (begin != end)) {
0451             --end;
0452             while ((begin != end) && ((*end)->isValid() == false))
0453               --end;
0454             ++end;
0455           }
0456
0457           // if we still have some hits build the track candidate
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 {  //all invalid hits have been already kicked out
0468             if ((end - begin) >= int(minimumHits_)) {
0469               output->push_back(makeCandidate(*trk, begin, end));
0470             }
0471           }
0472
0473           // if we still have some hits build the track candidate
0474           //if ((end - begin) >= int(minimumHits_)) {
0475           //    output->push_back( makeCandidate ( *trk, begin, end ) );
0476           //}
0477
0478           // now delete the hits not used by the candidate
0479           for (begin = hits.begin(), end = hits.end(); begin != end; ++begin) {
0480             if (*begin)
0481               delete *begin;
0482           }
0483           hits.clear();
0484         }  // loop on trajectories
0485
0486       } else {  //use plain tracks
0487
0488         // loop on tracks
0489         for (std::vector<reco::Track>::const_iterator ittrk = tracks->begin(), edtrk = tracks->end(); ittrk != edtrk;
0490              ++ittrk) {
0491           //    std::cout<<"The hit collection has size "<<hits.size()<<" (should be 0) while the track contains initially "<< ittrk->recHitsEnd() - ittrk->recHitsBegin()<<std::endl;
0492
0493           const Track *trk = &(*ittrk);
0494
0495           produceFromTrack(iSetup, trk, hits);
0496           //-----------------------
0497           /*
0498       std::cout<<"Hit collection in output has size "<<hits.size()<<". Dumping hit positions..."<<std::endl;
0499         for (std::vector<TrackingRecHit *>::iterator ith = hits.begin(), edh = hits.end(); ith != edh; ++ith) {
0500       const TrackingRecHit *myhit = *(ith);
0501         TransientTrackingRecHit::RecHitPointer ttrh;
0502         float radius=0.0;
0503         float xx=-999.0,yy=-999.0,zz=-999.0;
0504         unsigned int myid=0;
0505         if(myhit->isValid()){
0506           ttrh = RHBuilder->build(myhit);
0507           xx=ttrh->globalPosition().x();
0508           yy=ttrh->globalPosition().y();
0509           zz=ttrh->globalPosition().z();
0510           radius = sqrt( pow(xx,2)+pow(yy,2) );
0511           myid=myhit->geographicalId().rawId();
0512         }
0513         std::cout<<"-$-$ OUTPUT Hit position: ( "<<xx<<" , " <<yy<<" , " <<zz<<" ) , RADIUS = "  <<radius<<"  on DetID= "<< myid<<std::endl;
0514     }//end loop on hits
0515       */
0516           //-----------------------
0517
0518           std::vector<TrackingRecHit *>::iterator begin = hits.begin(), end = hits.end();
0519           // std::cout << "Back in the main producer (TRK), the final hit collection has size " << hits.size() << std::endl;
0520           // strip invalid hits at the beginning
0521           if (stripFrontInvalidHits_) {
0522             while ((begin != end) && ((*begin)->isValid() == false))
0523               ++begin;
0524           }
0525           // strip invalid hits at the end
0526           if (stripBackInvalidHits_ && (begin != end)) {
0527             --end;
0528             while ((begin != end) && ((*end)->isValid() == false))
0529               --end;
0530             ++end;
0531           }
0532
0533           // if we still have some hits build the track candidate
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 {  //all invalid hits have been already kicked out
0545             if ((end - begin) >= int(minimumHits_)) {
0546               output->push_back(makeCandidate(*ittrk, begin, end));
0547             }
0548           }
0549
0550           // now delete the hits not used by the candidate
0551           for (begin = hits.begin(), end = hits.end(); begin != end; ++begin) {
0552             if (*begin)
0553               delete *begin;
0554           }
0555           hits.clear();
0556         }  // loop on tracks
0557       }    //end else useTracks
0558
0559       // std::cout<<"OUTPUT SIZE: "<<output->size()<<std::endl;
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       //  double innerP=sqrt( pow(tk.innerMomentum().X(),2)+pow(tk.innerMomentum().Y(),2)+pow(tk.innerMomentum().Z(),2) );
0573       //  if ( (pdir == alongMomentum) == ( innerP >= tk.outerP() ) ) {
0574
0575       if ((pdir == alongMomentum) == ((tk.outerPosition() - tk.innerPosition()).Dot(tk.momentum()) >= 0)) {
0576         // use inner state
0577         TrajectoryStateOnSurface originalTsosIn(
0578             trajectoryStateTransform::innerStateOnSurface(tk, *theGeometry, &*theMagField));
0579         state = trajectoryStateTransform::persistentState(originalTsosIn, DetId(tk.innerDetId()));
0580       } else {
0581         // use outer state
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         //if(! (*hitsBegin)->isValid() ) std::cout<<"Putting in the trackcandidate an INVALID HIT !"<<std::endl;
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       // loop on tracks
0603       hits.clear();  // extra safety
0604
0605       for (trackingRecHit_iterator ith = itt->recHitsBegin(), edh = itt->recHitsEnd(); ith != edh; ++ith) {
0606         const TrackingRecHit *hit = (*ith);  // ith is an iterator on edm::Ref to rechit
0607
0608         DetId detid = hit->geographicalId();
0609
0610         //check that the hit is a real hit and not a constraint
0611         if (hit->isValid() && hit == nullptr && detid.rawId() == 0)
0612           continue;
0613
0614         int verdict = checkHit(iSetup, detid, hit);
0615         if (verdict == 0) {
0616           // just copy the hit
0617           hits.push_back(hit->clone());
0618         } else if (verdict < -2) {  //hit rejected because did not pass the selections
0619                                     // still, if replaceWithInactiveHits is true we have to put a new hit
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());  //hit not in the tracker
0625         else if (verdict == -1) {        //hit not valid
0626           if (!stripAllInvalidHits_) {
0627             hits.push_back(hit->clone());
0628           }
0629         }
0630       }  // loop on hits
0631
0632     }  //end TrackerTrackHitFilter::produceFromTrack
0633
0634     void TrackerTrackHitFilter::produceFromTrajectory(const edm::EventSetup &iSetup,
0635                                                       const Trajectory *itt,
0636                                                       std::vector<TrackingRecHit *> &hits) {
0637       hits.clear();  // extra safety
0638       nOverlaps = 0;
0639
0640       //Retrieve tracker topology from geometry
0641       edm::ESHandle<TrackerTopology> tTopoHand = iSetup.getHandle(tokenTrackerTopo);
0642       const TrackerTopology *tTopo = tTopoHand.product();
0643
0644       std::vector<TrajectoryMeasurement> tmColl = itt->measurements();
0645
0646       //---OverlapBegin needed eventually for overlaps, but I must create them here in any case
0647       const TrajectoryMeasurement *previousTM(nullptr);
0648       DetId previousId(0);
0649       //int previousLayer(-1);
0650       ///---OverlapEnd
0651
0652       for (std::vector<TrajectoryMeasurement>::const_iterator itTrajMeas = tmColl.begin(); itTrajMeas != tmColl.end();
0653            itTrajMeas++) {
0654         TransientTrackingRecHit::ConstRecHitPointer hitpointer = itTrajMeas->recHit();
0655
0656         //check that the hit is a real hit and not a constraint
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)) {  //check angle of track on module if requested
0667             verdict = -6;                                            //override previous verdicts
0668           }
0669         }
0670
0671         /*
0672     //this has been included in checkHitAngle(*itTrajMeas)
0673     if (verdict == 0) {
0674     if( PXLcorrClusChargeCut_>0.0  && !checkPXLCorrClustCharge(*itTrajMeas) ){//check angle of track on module if requested
0675     verdict=-7;//override previous verdicts
0676     }
0677     }
0678     */
0679
0680         if (verdict == 0) {  // Hit TAKEN !!!!
0681
0682           if (tagOverlaps_) {  ///---OverlapBegin
0683             //std::cout<<"Looking for overlaps in Run="<<iRun<<" , Event ="<<iEvt<<std::flush;
0684
0685             int side(sideFromId(detid, tTopo));    //side 0=barrel, 1=minus , 2=plus
0686             int layer(layerFromId(detid, tTopo));  //layer or disk
0687             int subDet = detid.subdetId();
0688             //std::cout  << "  Check Subdet #" <<subDet << ", layer = " <<layer<<" stereo: "<< ((subDet > 2)?(SiStripDetId(detid).stereo()):2);
0689
0690             if ((previousTM != nullptr) && (layer != -1)) {
0691               //std::cout<<"A previous TM exists! "<<std::endl;
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())) {  //if either pixel or strip stereo module
0705                   //  overlapHits.push_back(std::make_pair(&(*itmCompare),&(*itm)));
0706                   //std::cout<< "Adding pair "<< ((subDet >2)?(SiStripDetId(detid).stereo()):2)
0707                   //     << " from SubDet = "<<subDet<<" , layer = " << layer<<"  Run:"<<iRun<<"\tEv: "<<iEvt<<"\tId1: "<<compareId.rawId()<<"\tId2: "<<detid.rawId()<<std::endl;
0708                   // 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;
0709                   ///
0710                   // if(detid.rawId()<compareId.rawId()){
0711                   // std::cout<< "+++ "<< "\t"<<iRun<<"\t"<<iEvt<<"\t"<<detid.rawId()<<"\t"<<compareId.rawId()<<std::endl;
0712                   // }
0713                   //else  std::cout<< "+++ "<< "\t"<<iRun<<"\t"<<iEvt<<"\t"<<compareId.rawId()<<"\t"<<detid.rawId()<<std::endl;
0714
0715                   nOverlaps++;
0716                   break;
0717                 }
0718               }  //end second loop on TM for overlap tagging
0719
0720             }  //end   if ( (layer!=-1 )&&(acceptLayer[subDet]))
0721
0722             previousTM = &(*itTrajMeas);
0723             previousId = detid;
0724             //previousLayer = layer;
0725           }  //end if look for overlaps
0726           ///---OverlapEnd
0727
0728           hits.push_back(hit->clone());  //just copy it
0729         }                                //end if HIT TAKEN
0730         else if (verdict < -2) {         //hit rejected because did not pass the selections
0731           // still, if replaceWithInactiveHits is true we have to put a new hit
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());  //hit not in the tracker
0737         else if (verdict == -1) {        //hit not valid
0738           if (!stripAllInvalidHits_) {
0739             hits.push_back(hit->clone());
0740           }
0741         }
0742       }  // loop on hits
0743
0744       std::reverse(hits.begin(), hits.end());
0745     }  //end TrackerTrackHitFilter::produceFromTrajectories
0746
0747     int TrackerTrackHitFilter::checkHit(const edm::EventSetup &iSetup, const DetId &detid, const TrackingRecHit *hit) {
0748       //Retrieve tracker topology from geometry
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) {  // check for tracker hits
0755           bool verdict = true;
0756           // first check at structure level
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           // if the hit is good, check again at module level
0762           if (verdict) {
0763             if (std::binary_search(detsToIgnore_.begin(), detsToIgnore_.end(), detid.rawId())) {
0764               hitresult = -4;
0765             }
0766           } else
0767             hitresult = -3;
0768           //if the hit is in the desired part of the det, check other things
0769           if (hitresult == 0 && rejectBadStoNHits_) {
0770             if (!checkStoN(detid, hit))
0771               hitresult = -5;
0772           }  //end if S/N is ok
0773         }    //end hit in tracker
0774         else
0775           hitresult = -2;
0776       }  //end hit is valid
0777       else
0778         hitresult = -1;  //invalid hit
0779       return hitresult;
0780     }  //end  TrackerTrackHitFilter::checkHit()
0781
0782     bool TrackerTrackHitFilter::checkStoN(const DetId &id, const TrackingRecHit *therechit) {
0783       bool keepthishit = true;
0784       // const uint32_t& recHitDetId = id.rawId();
0785
0786       //check StoN only if subdet was set by the user
0787       //  int subdet_cnt=0;
0788       int subdet_cnt = id.subdetId();
0789
0790       //  for(subdet_cnt=0;subdet_cnt<6; ++subdet_cnt){
0791
0792       //  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
0793
0794       // for phase-2 OT placehold, do nothing
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           //check that hit is in a det belonging to a subdet where we decided to apply a S/N cut
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               // if it's matched rechit, use the average of stereo and mono clusters
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             //if(!keepthishit)std::cout<<"Hit rejected because of bad S/N: "<<siStripClusterInfo_->signalOverNoise()<<std::endl;
0883           }
0884         }  //end if  subdetStoN_[subdet_cnt]&&...
0885
0886       }  //end if subdet_cnt >2
0887       else if (GeomDetEnumerators::isTrackerPixel(theGeometry->geomDetSubDetector(id.subdetId()))) {  //pixel
0888         //pixels have naturally a very low noise (because of their low capacitance). So the S/N cut is
0889         //irrelevant in this case. Leave it dummy
0890         keepthishit = true;
0891
0892         /**************
0893        * Cut on cluster charge corr by angle embedded in the checkHitAngle() function
0894        *
0895        *************/
0896
0897         if (checkPXLQuality_) {
0898           const SiPixelRecHit *pixelhit = dynamic_cast<const SiPixelRecHit *>(therechit);
0899           if (pixelhit != nullptr) {
0900             //std::cout << "ClusterCharge=" <<std::flush<<pixelhit->cluster()->charge() << std::flush;
0901             float xyprob = pixelhit->clusterProbability(0);        //x-y combined log_e probability of the pixel cluster
0902                                                                    //singl x- and y-prob not stored sicne CMSSW 3_9_0
0903             float xychargeprob = pixelhit->clusterProbability(1);  //xy-prob * charge prob
0904             //  float chargeprob   =pixelhit->clusterProbability(2);//charge prob
0905             bool haspassed_tplreco = pixelhit->hasFilledProb();  //the cluster was associted to a template
0906             int qbin = pixelhit->qBin();  //R==meas_charge/pred_charge:  Qbin=0 ->R>1.5 , =1->1<R<1.5 ,=2->0.85<R<1 ,
0907                                           // Qbin=3->0.95*Qminpred<R<0.85 ,=4->, =5->meas_charge<0.95*Qminpred
0908
0909             //  if(haspassed_tplreco)   std::cout<<"  CLUSTPROB=\t"<<xprob<<"\t"<<yprob<<"\t"<<combprob<<"\t"<<qbin<<std::endl;
0910             //  else std::cout<<"CLUSTPROBNOTDEF=\t"<<xprob<<"\t"<<yprob<<"\t"<<combprob<<"\t"<<qbin<<std::endl;
0911
0912             keepthishit = false;
0913             //  std::cout<<"yyyyy "<<qbin<<" "<<xprob<<"  "<<yprob<<std::endl;
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         }  //end if check pixel quality flag
0922       }
0923       //    else  throw cms::Exception("TrackerTrackHitFilter") <<"Loop over subdetector out of range when applying the S/N cut: "<<subdet_cnt;
0924
0925       //  }//end loop on subdets
0926
0927       return keepthishit;
0928     }  //end CheckStoN
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   edm::LogDebug("TrackerTrackHitFilter")<<"TSOS parameters: ";
0936   edm::LogDebug("TrackerTrackHitFilter") <<"Global momentum: "<<tsos.globalMomentum().x()<<"  "<<tsos.globalMomentum().y()<<"  "<<tsos.globalMomentum().z();
0937   edm::LogDebug("TrackerTrackHitFilter") <<"Local momentum: "<<tsos.localMomentum().x()<<"  "<<tsos.localMomentum().y()<<"  "<<tsos.localMomentum().z();
0938   edm::LogDebug("TrackerTrackHitFilter") <<"Track charge: "  <<tsos.charge();
0939   edm::LogDebug("TrackerTrackHitFilter")<<"Local position: "  <<tsos.localPosition().x()<<"  "<<tsos.localPosition().y()<<"  "<<tsos.localPosition().z();
0940   */
0941       if (tsos.isValid()) {
0942         //check the angle of this tsos
0943         float mom_x = tsos.localDirection().x();
0944         float mom_y = tsos.localDirection().y();
0945         float mom_z = tsos.localDirection().z();
0946         //we took LOCAL momentum, i.e. respect to surface. Thus the plane is z=0
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;  // keep this hit
0950         // else  std::cout<<"Hit rejected because angle is "<< angle<<" ( <"<<TrackAngleCut_<<" )"<<std::endl;
0951
0952         if (angle_ok && PXLcorrClusChargeCut_ > 0.0) {
0953           //
0954           //get the hit from the TM and check that it is in the pixel
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()))) {  //do it only for pixel hits
0960               corrcharge_ok = false;
0961               float clust_alpha = atan2(mom_z, mom_x);
0962               float clust_beta = atan2(mom_z, mom_y);
0963
0964               //Now get the cluster charge
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               //std::cout<<"xxxxx "<<clust_charge<<" "<<corr_clust_charge<<"  " <<pixelhit->qBin()<<"  "<<pixelhit->clusterProbability(1)<<"  "<<pixelhit->clusterProbability(2)<< std::endl;
0971               if (corr_clust_charge > PXLcorrClusChargeCut_) {
0972                 corrcharge_ok = true;
0973               }
0974             }  //end if hit is in pixel
0975           }    //end if hit is valid
0976
0977         }  //check corr cluster charge for pixel hits
0978
0979       }  //end if TSOS is valid
0980       else {
0981         edm::LogWarning("TrackerTrackHitFilter") << "TSOS not valid ! Impossible to calculate track angle.";
0982       }
0983
0984       return angle_ok && corrcharge_ok;
0985     }  //end TrackerTrackHitFilter::checkHitAngle
0986
0987     bool TrackerTrackHitFilter::checkPXLCorrClustCharge(const TrajectoryMeasurement &meas) {
0988       /*
0989     Code taken from DPGAnalysis/SiPixelTools/plugins/PixelNtuplizer_RealData.cc
0990   */
0991
0992       bool corrcharge_ok = false;
0993       //get the hit from the TM and check that it is in the pixel
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()))) {  //SiStrip hit, skip
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         //Now get the cluster charge
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       }  //end if TSOS is valid
1021       return corrcharge_ok;
1022
1023     }  //end TrackerTrackHitFilter::checkPXLCorrClustCharge
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     // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
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   }  // namespace modules
1067 }  // namespace reco
1068
1069 // ========= MODULE DEF ==============
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);