CMS 3D CMS Logo

/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/RecoTracker/FinalTrackSelectors/src/TrackerTrackHitFilter.cc

Go to the documentation of this file.
00001 #include "FWCore/Framework/interface/EDProducer.h"
00002 #include "FWCore/Framework/interface/Event.h"
00003 #include "FWCore/Framework/interface/EventSetup.h"
00004 #include "FWCore/Framework/interface/ESHandle.h"
00005 
00006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00008 #include "FWCore/Utilities/interface/InputTag.h"
00009 
00010 #include "DataFormats/Common/interface/View.h"
00011 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00012 #include "DataFormats/TrackReco/interface/Track.h"
00013 #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
00014 #include "DataFormats/TrackingRecHit/interface/InvalidTrackingRecHit.h"
00015 
00016 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00017 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00018 #include "MagneticField/Engine/interface/MagneticField.h"
00019 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00020 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00021 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00022 
00023 #include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
00024 #include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
00025 #include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
00026 #include "DataFormats/SiStripDetId/interface/TIBDetId.h"
00027 #include "DataFormats/SiStripDetId/interface/TIDDetId.h"
00028 #include "DataFormats/SiStripDetId/interface/TOBDetId.h"
00029 #include "DataFormats/SiStripDetId/interface/TECDetId.h"
00030 
00031 //for S/N cut
00032 #include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
00033 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
00034 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
00035 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h"
00036 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
00037 #include "AnalysisDataFormats/SiStripClusterInfo/interface/SiStripClusterInfo.h"
00038 //for angle cut
00039 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00040 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
00041 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00042 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
00043 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
00044 #include "TMath.h"
00045 
00046 #include "TrackingTools/Records/interface/TransientRecHitRecord.h" 
00047 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
00048 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
00049 
00050 #include <boost/regex.hpp>
00051 #include <map>
00052 //#include <math>
00053 
00075 namespace reco { namespace modules {
00076   class TrackerTrackHitFilter : public edm::EDProducer {
00077   public:
00078     TrackerTrackHitFilter(const edm::ParameterSet &iConfig) ; 
00079     virtual void produce(edm::Event &iEvent, const edm::EventSetup &iSetup) ;
00080     int checkHit(const edm::EventSetup &iSetup,const  DetId &detid,  const TrackingRecHit * hit);
00081     void produceFromTrajectory( const edm::EventSetup &iSetup, const Trajectory *itt, std::vector<TrackingRecHit *>&hits);
00082     void produceFromTrack( const edm::EventSetup &iSetup, const Track *itt, std::vector<TrackingRecHit *>&hits);
00083   private:
00084     class Rule {
00085     public:
00086       // parse a rule from a string
00087       Rule(const std::string &str) ;
00088       // check this DetId, update the verdict if the rule has anything to say about it
00089       void apply(DetId detid, bool & verdict) const {
00090         // check detector
00091         if (detid.subdetId() == subdet_) {
00092           // check layer
00093           if ( (layer_ == 0) || (layer_ == layer(detid)) ) {
00094             // override verdict
00095             verdict = keep_;
00096             //  std::cout<<"Verdict is "<<verdict<<" for subdet "<<subdet_<<", layer "<<layer_<<"! "<<std::endl;
00097           }
00098           // else std::cout<<"No, sorry, wrong layer.Retry..."<<std::endl;
00099         }
00100         //      else{ std::cout<<"No, sorry, wrong subdet.Retry..."<<std::endl;}
00101       }
00102     private:
00103       int  subdet_;
00104 
00105       int  layer_;
00106       bool keep_;
00107       int layer(DetId detid) const ;
00108 
00109     };
00110     
00111     edm::InputTag src_;
00112 
00113     int iRun;
00114     int iEvt;
00115 
00116     size_t minimumHits_;
00117     
00118     bool replaceWithInactiveHits_;
00119     bool stripFrontInvalidHits_;
00120     bool stripBackInvalidHits_;
00121     bool stripAllInvalidHits_;
00122     
00123     bool rejectBadStoNHits_;
00124     std::string CMNSubtractionMode_;
00125     std::vector<bool>   subdetStoN_;//(6); //,std::bool(false));
00126     std::vector<double> subdetStoNlowcut_;//(6,-1.0);
00127     std::vector<double> subdetStoNhighcut_;//(6,-1.0);
00128     bool checkStoN( const edm::EventSetup &iSetup,const DetId &id, const TrackingRecHit* therechit);
00129     void parseStoN(const std::string &str);
00130     
00131     std::vector<uint32_t> detsToIgnore_;
00132     std::vector<Rule> rules_;
00133 
00134     bool useTrajectories_;
00135     bool rejectLowAngleHits_;
00136     double TrackAngleCut_;
00137     bool checkHitAngle(const TrajectoryMeasurement &meas);
00138     bool checkPXLQuality_;
00139     double pxlTPLProbXY_;
00140     double pxlTPLProbXYQ_;
00141     std::vector<int32_t> pxlTPLqBin_;
00142   
00143     bool checkPXLCorrClustCharge(const TrajectoryMeasurement &meas);
00144     double PXLcorrClusChargeCut_;
00145 
00146     edm::ESHandle<TrackerGeometry> theGeometry;
00147     edm::ESHandle<MagneticField>   theMagField;
00148 
00149     edm::ESHandle<TransientTrackingRecHitBuilder> theBuilder;
00150 
00151 
00152 
00153     bool tagOverlaps_;
00154     int nOverlaps;
00155     int layerFromId (const DetId& id) const;
00156     // bool checkOverlapHit();
00157 
00158     TrackCandidate makeCandidate(const reco::Track &tk, std::vector<TrackingRecHit *>::iterator hitsBegin, std::vector<TrackingRecHit *>::iterator hitsEnd) ;
00159     //const TransientTrackingRecHitBuilder *RHBuilder;
00160   }; // class
00161   
00162   
00163 
00164 TrackerTrackHitFilter::Rule::Rule(const std::string &str) {
00165     static boost::regex rule("(keep|drop)\\s+([A-Z]+)(\\s+(\\d+))?");
00166     boost::cmatch match;
00167     std::string match_1;
00168     std::string match_2;
00169     std::string match_3;
00170     // match and check it works
00171     if (!regex_match(str.c_str(), match, rule)) {
00172         throw cms::Exception("Configuration") << "Rule '" << str << "' not understood.\n";
00173     }
00174     else{
00175       std::cout<<"*** Rule Command given to TrackerTrackHitFilter:\t"<<str<<std::endl;
00176 
00177     }
00178     // Set up fields:
00179     //  rule type
00180     keep_  = (strncmp(match[1].first,"keep",4 )== 0);
00181    
00182     //  subdet
00183     subdet_ = -1;
00184     if      (strncmp(match[2].first, "PXB", 3) == 0) subdet_ = PixelSubdetector::PixelBarrel;
00185     else if (strncmp(match[2].first, "PXE", 3) == 0) subdet_ = PixelSubdetector::PixelEndcap;
00186     else if (strncmp(match[2].first, "TIB", 3) == 0) subdet_ = StripSubdetector::TIB;
00187     else if (strncmp(match[2].first, "TID", 3) == 0) subdet_ = StripSubdetector::TID;
00188     else if (strncmp(match[2].first, "TOB", 3) == 0) subdet_ = StripSubdetector::TOB;
00189     else if (strncmp(match[2].first, "TEC", 3) == 0) subdet_ = StripSubdetector::TEC;
00190     if (subdet_ == -1) {
00191         throw cms::Exception("Configuration") << "Detector '" << match[2].first << "' not understood. Should be PXB, PXE, TIB, TID, TOB, TEC.\n";
00192     }
00193     //   layer (if present)
00194     if (match[4].first != match[4].second) {
00195         layer_ = atoi(match[4].first);        
00196     } else {
00197         layer_ = 0;
00198     }
00199 }//end Rule::Rule
00200 
00201 int TrackerTrackHitFilter::Rule::layer(DetId detid) const {
00202     switch (detid.subdetId()) {
00203         case PixelSubdetector::PixelBarrel: return PXBDetId(detid).layer();
00204         case PixelSubdetector::PixelEndcap: return PXFDetId(detid).disk();
00205         case StripSubdetector::TIB:         return TIBDetId(detid).layer();
00206         case StripSubdetector::TID:         return TIDDetId(detid).wheel();
00207         case StripSubdetector::TOB:         return TOBDetId(detid).layer();
00208         case StripSubdetector::TEC:         return TECDetId(detid).wheel();
00209     }
00210     return -1; // never match
00211 }
00212 
00213 void TrackerTrackHitFilter::parseStoN(const std::string &str) {
00214   // match a set of capital case chars (preceded by an arbitrary number of leading blanks),
00215   //followed b an arbitrary number of blanks, one or more digits (not necessary, they cannot also be,
00216   // another set of blank spaces and, again another *eventual* digit  
00217   // static boost::regex rule("\\s+([A-Z]+)(\\s+(\\d+)(\\.)?(\\d+))?(\\s+(\\d+)(\\.)?(\\d+))?");
00218  static boost::regex rule("([A-Z]+)"
00219                           "\\s*(\\d+\\.*\\d*)?"
00220                           "\\s*(\\d+\\.*\\d*)?");
00221  
00222 
00223   boost::cmatch match;
00224   std::string match_1;
00225   std::string match_2;
00226   std::string match_3;
00227   // match and check it works
00228   if (!regex_match(str.c_str(), match, rule)) {
00229     throw cms::Exception("Configuration") << "Rule for S to N cut '" << str << "' not understood.\n";
00230   }
00231   else{
00232     std::string match_0=(match[0].first,match[0].second);
00233     match_1=(match[1].first,match[1].second);
00234     match_2=(match[2].first,match[2].second);
00235     match_3=(match[3].first,match[3].second);
00236 
00237   }
00238 
00239   int cnt=0;
00240   float subdet_ind[6];
00241   for(cnt=0;cnt<6;cnt++ ){
00242     subdet_ind[cnt]=-1.0;
00243   }
00244   
00245 
00246   bool doALL=false;
00247   std::string match_1a(match[1].first,match[1].second);
00248  if (strncmp(match[1].first, "ALL", 3) == 0) doALL=true;  
00249   if (doALL || strncmp(match[1].first, "PXB", 3) == 0)  subdet_ind[0] = +1.0;
00250   if (doALL || strncmp(match[1].first, "PXE", 3) == 0)  subdet_ind[1] = +1.0;
00251   if (doALL || strncmp(match[1].first, "TIB", 3) == 0)  subdet_ind[2] = +1.0;
00252   if (doALL || strncmp(match[1].first, "TID", 3) == 0)  subdet_ind[3] = +1.0;
00253   if (doALL || strncmp(match[1].first, "TOB", 3) == 0)  subdet_ind[4] = +1.0;
00254   if (doALL || strncmp(match[1].first, "TEC", 3) == 0)  subdet_ind[5] = +1.0;
00255 
00256   for(cnt=0;cnt<6;cnt++ ){//loop on subdets
00257     if(subdet_ind[cnt]>0.0){
00258       subdetStoN_[cnt]=true;
00259       if (match[2].first != match[2].second) {
00260         subdetStoNlowcut_[cnt] = atof(match[2].first);        
00261       }
00262       if (match[3].first != match[3].second ) {
00263         subdetStoNhighcut_[cnt] = atof(match[3].first);        
00264       }
00265       std::cout<<"Setting thresholds*&^ for subdet #"<<cnt+1<<" = "<<subdetStoNlowcut_[cnt]<<" - "<<subdetStoNhighcut_[cnt]<<std::endl;
00266     }
00267   }
00268 
00269   bool correct_regex=false;
00270   for(cnt=0;cnt<6;cnt++ ){//check that the regex was correct
00271     if(subdetStoN_[cnt])correct_regex=true;
00272   }
00273 
00274   if (!correct_regex) {
00275     throw cms::Exception("Configuration") << "Detector '" << match_1a << "' not understood in parseStoN. Should be PXB, PXE, TIB, TID, TOB, TEC.\n";
00276     }
00277 
00278   //std::cout<<"Reached end of parseStoN"<<std::endl;
00279 }//end parseStoN
00280 
00281 
00282 TrackerTrackHitFilter::TrackerTrackHitFilter(const edm::ParameterSet &iConfig) :
00283     src_(iConfig.getParameter<edm::InputTag>("src")),
00284     minimumHits_(iConfig.getParameter<uint32_t>("minimumHits")),
00285     replaceWithInactiveHits_(iConfig.getParameter<bool>("replaceWithInactiveHits")),
00286     stripFrontInvalidHits_(iConfig.getParameter<bool>("stripFrontInvalidHits")),
00287     stripBackInvalidHits_( iConfig.getParameter<bool>("stripBackInvalidHits") ),
00288     stripAllInvalidHits_(  iConfig.getParameter<bool>("stripAllInvalidHits")  ),
00289     rejectBadStoNHits_(  iConfig.getParameter<bool>("rejectBadStoNHits")  ),
00290     CMNSubtractionMode_(  iConfig.getParameter<std::string>("CMNSubtractionMode")  ),
00291     detsToIgnore_( iConfig.getParameter<std::vector<uint32_t> >("detsToIgnore") ),
00292     useTrajectories_( iConfig.getParameter<bool>("useTrajectories") ),
00293     rejectLowAngleHits_( iConfig.getParameter<bool>("rejectLowAngleHits") ),
00294     TrackAngleCut_( iConfig.getParameter<double>("TrackAngleCut") ),
00295     checkPXLQuality_(iConfig.getParameter<bool>("usePixelQualityFlag") ),
00296     pxlTPLProbXY_(iConfig.getParameter<double>("PxlTemplateProbXYCut")),
00297     pxlTPLProbXYQ_(iConfig.getParameter<double>("PxlTemplateProbXYChargeCut")),
00298     pxlTPLqBin_(iConfig.getParameter<std::vector<int32_t> >("PxlTemplateqBinCut")),
00299     PXLcorrClusChargeCut_(iConfig.getParameter<double>("PxlCorrClusterChargeCut")),
00300     tagOverlaps_( iConfig.getParameter<bool>("tagOverlaps") )
00301 {
00302 
00303 
00304     // sanity check 
00305     if (stripAllInvalidHits_ && replaceWithInactiveHits_) {
00306         throw cms::Exception("Configuration") << "Inconsistent Configuration: you can't set both 'stripAllInvalidHits' and 'replaceWithInactiveHits' to true\n";
00307     }
00308     if(rejectLowAngleHits_ && !useTrajectories_){
00309       throw cms::Exception("Configuration") << "Wrong configuration of TrackerTrackHitFilter. You cannot apply the cut on the track angle without using Trajectories!\n";
00310     }
00311     if (!useTrajectories_ && PXLcorrClusChargeCut_>0){
00312       throw cms::Exception("Configuration") << "Wrong configuration of TrackerTrackHitFilter. You cannot apply the cut on the corrected pixel cluster charge without using Trajectories!\n";
00313     }
00314 
00315     if(pxlTPLqBin_.size()>2){
00316       std::cout<<"Warning from TrackerTrackHitFilter: vector with qBin cuts has size > 2. Additional items will be ignored."<<std::endl;
00317     }
00318 
00319 
00320 
00321     // read and parse commands
00322     std::vector<std::string> str_rules = iConfig.getParameter<std::vector<std::string> >("commands");
00323     rules_.reserve(str_rules.size());
00324     for (std::vector<std::string>::const_iterator it = str_rules.begin(), ed = str_rules.end(); it != ed; ++it) {
00325         rules_.push_back(Rule(*it));
00326     }
00327 
00328     if(rejectBadStoNHits_){//commands for S/N cut
00329 
00330       subdetStoN_.reserve(6);
00331       subdetStoNlowcut_.reserve(6);
00332       subdetStoNhighcut_.reserve(6);
00333       int cnt=0;
00334       for(cnt=0;cnt<6;cnt++ ){
00335         subdetStoN_[cnt]=false;
00336         subdetStoNlowcut_[cnt]=-1.0;
00337         subdetStoNhighcut_[cnt]=-1.0;
00338       }
00339 
00340       std::vector<std::string> str_StoNrules = iConfig.getParameter<std::vector<std::string> >("StoNcommands");
00341        for (std::vector<std::string>::const_iterator str_StoN = str_StoNrules.begin(); str_StoN != str_StoNrules.end();  ++str_StoN) {
00342          parseStoN(*str_StoN);
00343        }
00345        std::cout<<"Finished parsing S/N. Applying following cuts to subdets:";
00346       for(cnt=0;cnt<6;cnt++ ){
00348         std::cout<<"Subdet #"<<cnt+1<<" -> "<<subdetStoNlowcut_[cnt]<<" , "<<subdetStoNhighcut_[cnt];
00349       }
00350     }//end if rejectBadStoNHits_
00351 
00352 
00353     if(rejectLowAngleHits_ )    std::cout<<"\nApplying cut on angle track = "<<TrackAngleCut_<<std::endl;
00354 
00355 
00356     // sort detids to ignore
00357     std::sort(detsToIgnore_.begin(), detsToIgnore_.end());
00358     
00359     // issue the produce<>
00360     produces<TrackCandidateCollection>();
00361 }
00362 
00363 void
00364 TrackerTrackHitFilter::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) 
00365 {
00366   //Dump Run and Event
00367   iRun=iEvent.id().run();
00368   iEvt=iEvent.id().event();
00369  
00370   // read with View, so we can read also a TrackRefVector
00371   edm::Handle<std::vector<reco::Track> > tracks;
00372   edm::Handle<TrajTrackAssociationCollection> assoMap;
00373 
00374   if(useTrajectories_)iEvent.getByLabel(src_,  assoMap);
00375   else iEvent.getByLabel(src_, tracks);
00376   
00377   // read from EventSetup
00378   iSetup.get<TrackerDigiGeometryRecord>().get(theGeometry);
00379   iSetup.get<IdealMagneticFieldRecord>().get(theMagField);
00380   //iSetup.get<TransientRecHitRecord>().get("WithTrackAngle",theBuilder);
00381   // RHBuilder=   theBuilder.product();
00382 
00383 
00384   // prepare output collection
00385   size_t  candcollsize;
00386   if(useTrajectories_)candcollsize= assoMap->size();
00387   else candcollsize=tracks->size();
00388   std::auto_ptr<TrackCandidateCollection> output(new TrackCandidateCollection());
00389 
00390   output->reserve(candcollsize);
00391   
00392   // working area and tools
00393   std::vector<TrackingRecHit *> hits;
00394 
00395   if(useTrajectories_){
00396     for (TrajTrackAssociationCollection::const_iterator itass = assoMap->begin();  itass != assoMap->end(); ++itass){
00397       
00398       const edm::Ref<std::vector<Trajectory> >traj = itass->key;//trajectory in the collection
00399       const reco::TrackRef tkref = itass->val;//associated track track in the collection
00400       //std::cout<<"The hit collection has size "<<hits.size()<<" (should be 0) while the track contains initially "<< tkref->recHitsEnd() - tkref->recHitsBegin()<<std::endl;
00401       
00402       const Track *trk = &(*tkref);
00403       const Trajectory * myTrajectory= &(*traj);      
00404       produceFromTrajectory(iSetup,myTrajectory,hits);
00405 
00406       std::vector<TrackingRecHit *>::iterator begin = hits.begin(), end = hits.end();
00407       
00408       // strip invalid hits at the beginning
00409       if (stripFrontInvalidHits_) {
00410         while ( (begin != end) && ( (*begin)->isValid() == false ) ) ++begin;
00411       }
00412       // strip invalid hits at the end
00413       if (stripBackInvalidHits_ && (begin != end)) {
00414         --end;
00415         while ( (begin != end) && ( (*end)->isValid() == false ) ) --end;
00416         ++end;
00417       }
00418  
00419          // if we still have some hits build the track candidate
00420       if(replaceWithInactiveHits_){
00421         int nvalidhits=0;
00422         for( std::vector<TrackingRecHit *>::iterator ithit=begin;ithit!=end;++ithit){
00423           if( (*ithit)->isValid())nvalidhits++;
00424         }
00425         if(nvalidhits >= int(minimumHits_)){
00426           output->push_back( makeCandidate ( *trk, begin, end ) );
00427         }
00428         nvalidhits=0;
00429       }
00430       else{//all invalid hits have been already kicked out
00431         if ((end - begin) >= int(minimumHits_)) {
00432           output->push_back( makeCandidate ( *trk, begin, end ) );
00433         } 
00434       }
00435 
00436 
00437       // if we still have some hits build the track candidate
00438       //if ((end - begin) >= int(minimumHits_)) {
00439       //        output->push_back( makeCandidate ( *trk, begin, end ) );
00440       //}
00441 
00442  
00443       // now delete the hits not used by the candidate
00444       for (begin = hits.begin(), end = hits.end(); begin != end; ++begin) {
00445         if (*begin) delete *begin;
00446       } 
00447       hits.clear();
00448     } // loop on trajectories
00449     
00450     
00451   }
00452   else{ //use plain tracks
00453   
00454     // loop on tracks
00455     for (std::vector<reco::Track>::const_iterator ittrk = tracks->begin(), edtrk = tracks->end(); ittrk != edtrk; ++ittrk) {
00456 
00457       //    std::cout<<"The hit collection has size "<<hits.size()<<" (should be 0) while the track contains initially "<< ittrk->recHitsEnd() - ittrk->recHitsBegin()<<std::endl;
00458 
00459       const Track *trk = &(*ittrk);
00460 
00461       produceFromTrack(iSetup,trk,hits);
00462       //-----------------------
00463       /*
00464       std::cout<<"Hit collection in output has size "<<hits.size()<<". Dumping hit positions..."<<std::endl;
00465         for (std::vector<TrackingRecHit *>::iterator ith = hits.begin(), edh = hits.end(); ith != edh; ++ith) {
00466           const TrackingRecHit *myhit = *(ith);
00467             TransientTrackingRecHit::RecHitPointer ttrh;
00468             float radius=0.0;
00469             float xx=-999.0,yy=-999.0,zz=-999.0;
00470             unsigned int myid=0;
00471             if(myhit->isValid()){
00472               ttrh = RHBuilder->build(myhit);
00473               xx=ttrh->globalPosition().x();
00474               yy=ttrh->globalPosition().y();
00475               zz=ttrh->globalPosition().z();
00476               radius = sqrt( pow(xx,2)+pow(yy,2) );
00477               myid=myhit->geographicalId().rawId();
00478             }
00479             std::cout<<"-$-$ OUTPUT Hit position: ( "<<xx<<" , " <<yy<<" , " <<zz<<" ) , RADIUS = "  <<radius<<"  on DetID= "<< myid<<std::endl;
00480         }//end loop on hits
00481       */ 
00482       //-----------------------
00483 
00484 
00485       std::vector<TrackingRecHit *>::iterator begin = hits.begin(), end = hits.end();
00486       // std::cout << "Back in the main producer (TRK), the final hit collection has size " << hits.size() << std::endl;
00487       // strip invalid hits at the beginning
00488       if (stripFrontInvalidHits_) {
00489         while ( (begin != end) && ( (*begin)->isValid() == false ) ) ++begin;
00490       }
00491       // strip invalid hits at the end
00492       if (stripBackInvalidHits_ && (begin != end)) {
00493         --end;
00494         while ( (begin != end) && ( (*end)->isValid() == false ) ) --end;
00495         ++end;
00496       }
00497       
00498       // if we still have some hits build the track candidate
00499       if(replaceWithInactiveHits_){
00500         int nvalidhits=0;
00501         for( std::vector<TrackingRecHit *>::iterator ithit=begin;ithit!=end;++ithit){
00502           if( (*ithit)->isValid())nvalidhits++;
00503         }
00504         if(nvalidhits >= int(minimumHits_)){
00505           output->push_back( makeCandidate ( *ittrk, begin, end ) );
00506         }
00507 
00508       }
00509       else{//all invalid hits have been already kicked out
00510         if ((end - begin) >= int(minimumHits_)) {
00511           output->push_back( makeCandidate ( *ittrk, begin, end ) );
00512         } 
00513       }
00514 
00515       // now delete the hits not used by the candidate
00516       for (begin = hits.begin(), end = hits.end(); begin != end; ++begin) {
00517         if (*begin) delete *begin;
00518       } 
00519       hits.clear();
00520     } // loop on tracks
00521   }//end else useTracks
00522 
00523   // std::cout<<"OUTPUT SIZE: "<<output->size()<<std::endl;
00524   
00525   iEvent.put(output);
00526 }
00527 
00528 TrackCandidate
00529 TrackerTrackHitFilter::makeCandidate(const reco::Track &tk, std::vector<TrackingRecHit *>::iterator hitsBegin, std::vector<TrackingRecHit *>::iterator hitsEnd) {
00530 
00531     TrajectoryStateTransform transform;
00532     PropagationDirection   pdir = tk.seedDirection();
00533     PTrajectoryStateOnDet *state;
00534     if ( pdir == anyDirection ) throw cms::Exception("UnimplementedFeature") << "Cannot work with tracks that have 'anyDirecton' \n";
00535 
00536     //  double innerP=sqrt( pow(tk.innerMomentum().X(),2)+pow(tk.innerMomentum().Y(),2)+pow(tk.innerMomentum().Z(),2) );
00537     //  if ( (pdir == alongMomentum) == ( innerP >= tk.outerP() ) ) {
00538 
00539     if ( (pdir == alongMomentum) == (  (tk.outerPosition()-tk.innerPosition()).Dot(tk.momentum()) >= 0    ) ) {
00540         // use inner state
00541         TrajectoryStateOnSurface originalTsosIn(transform.innerStateOnSurface(tk, *theGeometry, &*theMagField));
00542         state = transform.persistentState( originalTsosIn, DetId(tk.innerDetId()) );
00543     } else { 
00544         // use outer state
00545         TrajectoryStateOnSurface originalTsosOut(transform.outerStateOnSurface(tk, *theGeometry, &*theMagField));
00546         state = transform.persistentState( originalTsosOut, DetId(tk.outerDetId()) );
00547     }
00548     TrajectorySeed seed(*state, TrackCandidate::RecHitContainer(), pdir);
00549     TrackCandidate::RecHitContainer ownHits;
00550     ownHits.reserve(hitsEnd - hitsBegin);
00551     for ( ; hitsBegin != hitsEnd; ++hitsBegin) { 
00552       //if(! (*hitsBegin)->isValid() ) std::cout<<"Putting in the trackcandidate an INVALID HIT !"<<std::endl;
00553       ownHits.push_back( *hitsBegin ); 
00554     }
00555         
00556     TrackCandidate cand(ownHits, seed, *state, tk.seedRef());
00557     delete state;
00558 
00559     return cand;
00560 }
00561 
00562 void TrackerTrackHitFilter::produceFromTrack(const edm::EventSetup &iSetup, const Track *itt , std::vector<TrackingRecHit *>&hits){
00563 
00564     // loop on tracks
00565         hits.clear(); // extra safety
00566 
00567         for (trackingRecHit_iterator ith = itt->recHitsBegin(), edh = itt->recHitsEnd(); ith != edh; ++ith) {
00568           const TrackingRecHit * hit = ith->get(); // ith is an iterator on edm::Ref to rechit
00569 
00570             DetId detid = hit->geographicalId();
00571           
00572             //check that the hit is a real hit and not a constraint
00573             if(hit->isValid() && hit==0 && detid.rawId()==0) continue;
00574 
00575             int verdict=checkHit(iSetup,detid,hit);
00576             if (verdict == 0) {
00577               // just copy the hit
00578               hits.push_back(hit->clone());
00579             }
00580             else if(verdict<-2){//hit rejected because did not pass the selections
00581                                 // still, if replaceWithInactiveHits is true we have to put a new hit
00582               if (replaceWithInactiveHits_) {
00583                 hits.push_back(new InvalidTrackingRecHit(detid, TrackingRecHit::inactive));
00584               } 
00585             }
00586             else if(verdict==-2) hits.push_back(hit->clone());//hit not in the tracker
00587             else if(verdict==-1){ //hit not valid
00588               if (!stripAllInvalidHits_) {
00589                 hits.push_back(hit->clone());
00590               } 
00591             } 
00592         } // loop on hits
00593                         
00594 }//end TrackerTrackHitFilter::produceFromTrack
00595 
00596 
00597 void TrackerTrackHitFilter::produceFromTrajectory(const edm::EventSetup &iSetup, const Trajectory *itt, std::vector<TrackingRecHit *>&hits){
00598   hits.clear(); // extra safety
00599   nOverlaps=0;
00600 
00601   
00602   std::vector<TrajectoryMeasurement> tmColl =itt->measurements();
00603 
00604   //---OverlapBegin needed eventually for overlaps, but I must create them here in any case
00605   const TrajectoryMeasurement* previousTM(0);
00606   DetId previousId(0);
00607   int previousLayer(-1);
00609 
00610   int constrhits=0;
00611 
00612   for(std::vector<TrajectoryMeasurement>::const_iterator itTrajMeas = tmColl.begin(); itTrajMeas!=tmColl.end(); itTrajMeas++){
00613      TransientTrackingRecHit::ConstRecHitPointer hitpointer = itTrajMeas->recHit();
00614 
00615      //check that the hit is a real hit and not a constraint
00616      if(hitpointer->isValid() && hitpointer->hit()==0){constrhits++; continue;}
00617     
00618     const TrackingRecHit *hit=((*hitpointer).hit());
00619     DetId detid = hit->geographicalId();
00620     int verdict=checkHit(iSetup,detid,hit);
00621   
00622     if (verdict == 0) {
00623       if( rejectLowAngleHits_ && !checkHitAngle(*itTrajMeas) ){//check angle of track on module if requested
00624         verdict=-6;//override previous verdicts
00625       }
00626     }
00627 
00628     /*
00629     //this has been included in checkHitAngle(*itTrajMeas) 
00630     if (verdict == 0) {
00631     if( PXLcorrClusChargeCut_>0.0  && !checkPXLCorrClustCharge(*itTrajMeas) ){//check angle of track on module if requested
00632     verdict=-7;//override previous verdicts
00633     }
00634     }
00635     */
00636 
00637 
00638     if(verdict==0){// Hit TAKEN !!!!
00639 
00640         if(tagOverlaps_){       
00641           //std::cout<<"Looking for overlaps in Run="<<iRun<<" , Event ="<<iEvt<<std::flush;
00642 
00643           int layer(layerFromId(detid));//layer 1-4=TIB, layer 5-10=TOB
00644           int subDet = detid.subdetId();
00645           //std::cout  << "  Check Subdet #" <<subDet << ", layer = " <<layer<<" stereo: "<< ((subDet > 2)?(SiStripDetId(detid).stereo()):2);
00646           
00647             if ( ( previousTM!=0 )&& (layer!=-1 )) {
00648               //std::cout<<"A previous TM exists! "<<std::endl;
00649               for (std::vector<TrajectoryMeasurement>::const_iterator itmCompare =itTrajMeas-1;itmCompare >= tmColl.begin() &&  itmCompare > itTrajMeas - 4;--itmCompare){
00650 
00651                 DetId compareId = itmCompare->recHit()->geographicalId();
00652                 if ( subDet != compareId.subdetId() ||
00653                      layer  != layerFromId(compareId)) break;
00654                 if (!itmCompare->recHit()->isValid()) continue;
00655                 if ( (subDet<=2) ||
00656                      (subDet > 2 && SiStripDetId(detid).stereo()==SiStripDetId(compareId).stereo()))
00657                   {//if either pixel or strip stereo module
00658                     //  overlapHits.push_back(std::make_pair(&(*itmCompare),&(*itm)));
00659                     //std::cout<< "Adding pair "<< ((subDet >2)?(SiStripDetId(detid).stereo()):2)
00660                     //     << " from SubDet = "<<subDet<<" , layer = " << layer<<"  Run:"<<iRun<<"\tEv: "<<iEvt<<"\tId1: "<<compareId.rawId()<<"\tId2: "<<detid.rawId()<<std::endl;
00661                     // 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;
00663                     // if(detid.rawId()<compareId.rawId()){
00664                       // std::cout<< "+++ "<< "\t"<<iRun<<"\t"<<iEvt<<"\t"<<detid.rawId()<<"\t"<<compareId.rawId()<<std::endl;
00665                     // }
00666                     //else  std::cout<< "+++ "<< "\t"<<iRun<<"\t"<<iEvt<<"\t"<<compareId.rawId()<<"\t"<<detid.rawId()<<std::endl;
00667                     
00668                     nOverlaps++;
00669                     break;
00670                   }
00671               }//end second loop on TM for overlap tagging
00672               
00673             }//end   if ( (layer!=-1 )&&(acceptLayer[subDet]))
00674 
00675             previousTM = &(* itTrajMeas);
00676             previousId = detid;
00677             previousLayer = layer;
00678         }//end if look for overlaps
00680 
00681         hits.push_back(hit->clone());   //just copy it 
00682     }//end if HIT TAKEN
00683     else if(verdict<-2){//hit rejected because did not pass the selections
00684       // still, if replaceWithInactiveHits is true we have to put a new hit
00685       if (replaceWithInactiveHits_) {
00686         hits.push_back(new InvalidTrackingRecHit(detid, TrackingRecHit::inactive));
00687       } 
00688     }
00689     else if(verdict==-2) hits.push_back(hit->clone());//hit not in the tracker
00690     else if(verdict==-1){ //hit not valid
00691       if (!stripAllInvalidHits_) {
00692         hits.push_back(hit->clone());
00693       } 
00694     }
00695   } // loop on hits
00696 
00697   std::reverse(hits.begin(),hits.end());
00698   //  std::cout<<"Finished producefromTrajecotries. Nhits in final coll"<<hits.size() <<"   Nconstraints="<<constrhits<<std::endl; 
00699 } //end TrackerTrackHitFilter::produceFromTrajectories
00700 
00701 int TrackerTrackHitFilter::checkHit(const edm::EventSetup &iSetup,const  DetId &detid,  const TrackingRecHit * hit){
00702 
00703 
00704   int hitresult=0;
00705   if (hit->isValid()) { 
00706  
00707     if (detid.det() == DetId::Tracker) {  // check for tracker hits
00708        bool  verdict = true;
00709       // first check at structure level
00710       for (std::vector<Rule>::const_iterator itr = rules_.begin(), edr = rules_.end(); itr != edr; ++itr) {
00711         itr->apply(detid, verdict);
00712       }
00713  
00714       // if the hit is good, check again at module level
00715       if ( verdict){
00716         if(std::binary_search(detsToIgnore_.begin(), detsToIgnore_.end(), detid.rawId())) {
00717           hitresult=-4;
00718         }
00719       }
00720       else hitresult=-3;
00721        //if the hit is in the desired part of the det, check other things
00722       if( hitresult==0 && rejectBadStoNHits_){
00723         if(!checkStoN(iSetup, detid, hit))hitresult=-5;
00724       }//end if S/N is ok
00725     }//end hit in tracker
00726     else hitresult =-2;
00727   }//end hit is valid
00728   else hitresult = -1; //invalid hit
00729   return hitresult;
00730 }//end  TrackerTrackHitFilter::checkHit()
00731 
00732 
00733 
00734 bool TrackerTrackHitFilter::checkStoN(const edm::EventSetup &iSetup, const DetId &id, const TrackingRecHit* therechit){
00735 
00736   bool  keepthishit=true;
00737   // const uint32_t& recHitDetId = id.rawId();
00738 
00739    //check StoN only if subdet was set by the user
00740   //  int subdet_cnt=0;
00741   int subdet_cnt=id.subdetId();
00742  
00743   //  for(subdet_cnt=0;subdet_cnt<6; ++subdet_cnt){
00744  
00745   //  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
00746 
00747 
00748     if(subdet_cnt>2){ //SiStrip
00749       if( subdetStoN_[subdet_cnt-1]){//check that hit is in a det belonging to a subdet where we decided to apply a S/N cut
00750         const std::type_info &type = typeid(*therechit);
00751         const SiStripCluster* cluster;
00752         if (type == typeid(SiStripRecHit2D)) {
00753           const SiStripRecHit2D* hit = dynamic_cast<const SiStripRecHit2D*>(therechit);
00754           if (hit!=0)    cluster = &*(hit->cluster());
00755           else{
00756             edm::LogError("TrackerTrackHitFilter")<< "TrackerTrackHitFilter::checkStoN : Unknown valid tracker hit in subdet " << id.subdetId()<< "(detID="<<id.rawId()<<")\n ";
00757             keepthishit = false;
00758           }
00759         }
00760         else if (type == typeid(SiStripRecHit1D)) {
00761           const SiStripRecHit1D* hit = dynamic_cast<const SiStripRecHit1D*>(therechit);
00762           if (hit!=0)    cluster = &*(hit->cluster());
00763           else{
00764             edm::LogError("TrackerTrackHitFilter")<< "TrackerTrackHitFilter::checkStoN : Unknown valid tracker hit in subdet " << id.subdetId()<< "(detID="<<id.rawId()<<")\n ";
00765             keepthishit = false;
00766           }
00767         }
00768         //the following two cases should not happen anymore since CMSSW > 2_0_X because of hit splitting in stereo modules
00769         //const SiStripMatchedRecHit2D* matchedhit = dynamic_cast<const SiStripMatchedRecHit2D*>(therechit);
00770         //const ProjectedSiStripRecHit2D* unmatchedhit = dynamic_cast<const ProjectedSiStripRecHit2D*>(therechit);
00771         else{
00772           throw cms::Exception("Unknown RecHit Type") << "RecHit of type " << type.name() << 
00773             " not supported. (use c++filt to demangle the name)";
00774         }
00775         
00776         if(keepthishit){
00777           SiStripClusterInfo clusterInfo = SiStripClusterInfo( *cluster, iSetup); 
00778           if ( (subdetStoNlowcut_[subdet_cnt-1]>0) && (clusterInfo.signalOverNoise() < subdetStoNlowcut_[subdet_cnt-1])  ) keepthishit = false; 
00779           if ( (subdetStoNhighcut_[subdet_cnt-1]>0) && (clusterInfo.signalOverNoise() > subdetStoNhighcut_[subdet_cnt-1])  ) keepthishit = false;       
00780           //if(!keepthishit)std::cout<<"Hit rejected because of bad S/N: "<<clusterInfo.signalOverNoise()<<std::endl;
00781         }
00782         
00783       }//end if  subdetStoN_[subdet_cnt]&&...
00784     
00785     }//end if subdet_cnt >2
00786     else if (subdet_cnt<=2){//pixel 
00787       //pixels have naturally a very low noise (because of their low capacitance). So the S/N cut is 
00788       //irrelevant in this case. Leave it dummy
00789       keepthishit = true;
00790 
00791       /**************
00792        * Cut on cluster charge corr by angle embedded in the checkHitAngle() function
00793        * 
00794        *************/
00795 
00796       if(checkPXLQuality_){
00797       const SiPixelRecHit* pixelhit = dynamic_cast<const SiPixelRecHit*>(therechit);
00798       if(pixelhit!=0){
00799         //std::cout << "ClusterCharge=" <<std::flush<<pixelhit->cluster()->charge() << std::flush;
00800         float xyprob=pixelhit->clusterProbability(0);//x-y combined log_e probability of the pixel cluster
00801                                                        //singl x- and y-prob not stored sicne CMSSW 3_9_0
00802         float xychargeprob   =pixelhit->clusterProbability(1);//xy-prob * charge prob
00803         //      float chargeprob   =pixelhit->clusterProbability(2);//charge prob
00804         bool haspassed_tplreco= pixelhit->hasFilledProb(); //the cluster was associted to a template
00805         int qbin      =pixelhit->qBin(); //R==meas_charge/pred_charge:  Qbin=0 ->R>1.5 , =1->1<R<1.5 ,=2->0.85<R<1 ,
00806                                          // Qbin=3->0.95*Qminpred<R<0.85 ,=4->, =5->meas_charge<0.95*Qminpred
00807 
00808         //      if(haspassed_tplreco)   std::cout<<"  CLUSTPROB=\t"<<xprob<<"\t"<<yprob<<"\t"<<combprob<<"\t"<<qbin<<std::endl;
00809         //      else std::cout<<"CLUSTPROBNOTDEF=\t"<<xprob<<"\t"<<yprob<<"\t"<<combprob<<"\t"<<qbin<<std::endl;
00810 
00811         keepthishit = false;
00812         //      std::cout<<"yyyyy "<<qbin<<" "<<xprob<<"  "<<yprob<<std::endl;
00813         if( haspassed_tplreco && xyprob>pxlTPLProbXY_ && xychargeprob>pxlTPLProbXYQ_ && qbin>pxlTPLqBin_[0] && qbin<=pxlTPLqBin_[1] )keepthishit = true;
00814 
00815       }
00816       else {std::cout<<"HIT IN PIXEL ("<<subdet_cnt <<") but PixelRecHit is EMPTY!!!"<<std::endl;}
00817       }//end if check pixel quality flag
00818     }
00819     //    else  throw cms::Exception("TrackerTrackHitFilter") <<"Loop over subdetector out of range when applying the S/N cut: "<<subdet_cnt;
00820 
00821     //  }//end loop on subdets
00822 
00823 
00824   return  keepthishit;  
00825 }//end CheckStoN
00826 
00827 
00828 
00829 bool TrackerTrackHitFilter::checkHitAngle(const TrajectoryMeasurement &meas){
00830   
00831   bool angle_ok=false;
00832   bool corrcharge_ok=true;
00833   TrajectoryStateOnSurface tsos = meas.updatedState();
00834   /*  
00835   edm::LogDebug("TrackerTrackHitFilter")<<"TSOS parameters: ";
00836   edm::LogDebug("TrackerTrackHitFilter") <<"Global momentum: "<<tsos.globalMomentum().x()<<"  "<<tsos.globalMomentum().y()<<"  "<<tsos.globalMomentum().z();
00837   edm::LogDebug("TrackerTrackHitFilter") <<"Local momentum: "<<tsos.localMomentum().x()<<"  "<<tsos.localMomentum().y()<<"  "<<tsos.localMomentum().z();
00838   edm::LogDebug("TrackerTrackHitFilter") <<"Track charge: "  <<tsos.charge();
00839   edm::LogDebug("TrackerTrackHitFilter")<<"Local position: "  <<tsos.localPosition().x()<<"  "<<tsos.localPosition().y()<<"  "<<tsos.localPosition().z();
00840   */
00841   if(tsos.isValid()){
00842     //check the angle of this tsos
00843     float mom_x=tsos.localDirection().x();
00844     float mom_y=tsos.localDirection().y();
00845     float mom_z=tsos.localDirection().z();
00846     //we took LOCAL momentum, i.e. respect to surface. Thus the plane is z=0
00847     float angle=TMath::ASin(TMath::Abs(mom_z) / sqrt(pow(mom_x,2)+pow(mom_y,2)+pow(mom_z,2) )  );
00848     if(!rejectLowAngleHits_ || angle>=TrackAngleCut_) angle_ok=true;// keep this hit
00849     // else  std::cout<<"Hit rejected because angle is "<< angle<<" ( <"<<TrackAngleCut_<<" )"<<std::endl;
00850 
00851     if(angle_ok &&  PXLcorrClusChargeCut_>0.0){
00852       //
00853       //get the hit from the TM and check that it is in the pixel
00854       TransientTrackingRecHit::ConstRecHitPointer hitpointer = meas.recHit();
00855       if(hitpointer->isValid()){
00856       const TrackingRecHit *hit=(*hitpointer).hit();
00857       if( (hit->geographicalId()).subdetId()<=2  ){//do it only for pixel hits
00858         corrcharge_ok=false;
00859         float clust_alpha= atan2( mom_z, mom_x );
00860         float clust_beta=  atan2( mom_z, mom_y );
00861         
00862         //Now get the cluster charge
00863         
00864         const SiPixelRecHit* pixelhit = dynamic_cast<const SiPixelRecHit*>(hit);
00865         float clust_charge=pixelhit->cluster()->charge();
00866         float corr_clust_charge=clust_charge *  sqrt( 1.0 / ( 1.0/pow( tan(clust_alpha), 2 ) + 
00867                                                               1.0/pow( tan(clust_beta ), 2 ) + 
00868                                                               1.0 )
00869                                                       );
00870         //std::cout<<"xxxxx "<<clust_charge<<" "<<corr_clust_charge<<"  " <<pixelhit->qBin()<<"  "<<pixelhit->clusterProbability(1)<<"  "<<pixelhit->clusterProbability(2)<< std::endl;
00871         if(corr_clust_charge>PXLcorrClusChargeCut_){
00872           corrcharge_ok=true;
00873         }
00874       }       //end if hit is in pixel
00875       }//end if hit is valid
00876       
00877     }//check corr cluster charge for pixel hits
00878 
00879   }//end if TSOS is valid
00880   else{
00881     edm::LogWarning("TrackerTrackHitFilter") <<"TSOS not valid ! Impossible to calculate track angle.";
00882   }
00883   
00884   return angle_ok&&corrcharge_ok; 
00885 }//end TrackerTrackHitFilter::checkHitAngle
00886 
00887 
00888 bool TrackerTrackHitFilter::checkPXLCorrClustCharge(const TrajectoryMeasurement &meas){
00889   /*
00890     Code taken from DPGAnalysis/SiPixelTools/plugins/PixelNtuplizer_RealData.cc
00891   */
00892 
00893   bool corrcharge_ok=false;
00894   //get the hit from the TM and check that it is in the pixel
00895   TransientTrackingRecHit::ConstRecHitPointer hitpointer = meas.recHit();
00896   if(!hitpointer->isValid()) return corrcharge_ok;
00897   const TrackingRecHit *hit=(*hitpointer).hit();
00898   if( (hit->geographicalId()).subdetId()>2  ){//SiStrip hit, skip
00899      return corrcharge_ok;
00900   }
00901 
00902   TrajectoryStateOnSurface tsos = meas.updatedState();  
00903   if(tsos.isValid()){
00904     float mom_x=tsos.localDirection().x();
00905     float mom_y=tsos.localDirection().y();
00906     float mom_z=tsos.localDirection().z();
00907     float clust_alpha= atan2( mom_z, mom_x );
00908     float clust_beta=  atan2( mom_z, mom_y );
00909 
00910     //Now get the cluster charge
00911  
00912       const SiPixelRecHit* pixelhit = dynamic_cast<const SiPixelRecHit*>(hit);
00913       float clust_charge=pixelhit->cluster()->charge();
00914       float corr_clust_charge=clust_charge *  sqrt( 1.0 / ( 1.0/pow( tan(clust_alpha), 2 ) + 
00915                                                           1.0/pow( tan(clust_beta ), 2 ) + 
00916                                                           1.0 )
00917                                                   );
00918       if(corr_clust_charge>PXLcorrClusChargeCut_)corrcharge_ok=true;
00919     
00920   }//end if TSOS is valid
00921   return corrcharge_ok;
00922 
00923 }//end TrackerTrackHitFilter::checkPXLCorrClustCharge
00924 
00925 
00926 
00927 int TrackerTrackHitFilter::layerFromId (const DetId& id) const
00928 {
00929  if ( id.subdetId()== int(PixelSubdetector::PixelBarrel) ) {
00930     PXBDetId tobId(id);
00931     return tobId.layer();
00932   }
00933   else if ( id.subdetId()== int(PixelSubdetector::PixelEndcap) ) {
00934     PXFDetId tobId(id);
00935     return tobId.disk() + (3*(tobId.side()-1));
00936   }
00937   else if ( id.subdetId()==StripSubdetector::TIB ) {
00938     TIBDetId tibId(id);
00939     return tibId.layer();
00940   }
00941   else if ( id.subdetId()==StripSubdetector::TOB ) {
00942     TOBDetId tobId(id);
00943     return tobId.layer();
00944   }
00945   else if ( id.subdetId()==StripSubdetector::TEC ) {
00946     TECDetId tobId(id);
00947     return tobId.wheel() + (9*(tobId.side()-1));
00948   }
00949   else if ( id.subdetId()==StripSubdetector::TID ) {
00950     TIDDetId tobId(id);
00951     return tobId.wheel() + (3*(tobId.side()-1));
00952   }
00953   return -1;
00954 
00955 }
00956 
00957 }} //namespaces
00958 
00959 
00960 // ========= MODULE DEF ==============
00961 #include "FWCore/PluginManager/interface/ModuleDef.h"
00962 #include "FWCore/Framework/interface/MakerMacros.h"
00963 
00964 using reco::modules::TrackerTrackHitFilter;
00965 DEFINE_FWK_MODULE(TrackerTrackHitFilter);