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