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
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
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
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
00087 Rule(const std::string &str) ;
00088
00089 void apply(DetId detid, bool & verdict) const {
00090
00091 if (detid.subdetId() == subdet_) {
00092
00093 if ( (layer_ == 0) || (layer_ == layer(detid)) ) {
00094
00095 verdict = keep_;
00096
00097 }
00098
00099 }
00100
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_;
00126 std::vector<double> subdetStoNlowcut_;
00127 std::vector<double> subdetStoNhighcut_;
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
00157
00158 TrackCandidate makeCandidate(const reco::Track &tk, std::vector<TrackingRecHit *>::iterator hitsBegin, std::vector<TrackingRecHit *>::iterator hitsEnd) ;
00159
00160 };
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
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
00179
00180 keep_ = (strncmp(match[1].first,"keep",4 )== 0);
00181
00182
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
00194 if (match[4].first != match[4].second) {
00195 layer_ = atoi(match[4].first);
00196 } else {
00197 layer_ = 0;
00198 }
00199 }
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;
00211 }
00212
00213 void TrackerTrackHitFilter::parseStoN(const std::string &str) {
00214
00215
00216
00217
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
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++ ){
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++ ){
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
00279 }
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
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
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_){
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 }
00351
00352
00353 if(rejectLowAngleHits_ ) std::cout<<"\nApplying cut on angle track = "<<TrackAngleCut_<<std::endl;
00354
00355
00356
00357 std::sort(detsToIgnore_.begin(), detsToIgnore_.end());
00358
00359
00360 produces<TrackCandidateCollection>();
00361 }
00362
00363 void
00364 TrackerTrackHitFilter::produce(edm::Event &iEvent, const edm::EventSetup &iSetup)
00365 {
00366
00367 iRun=iEvent.id().run();
00368 iEvt=iEvent.id().event();
00369
00370
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
00378 iSetup.get<TrackerDigiGeometryRecord>().get(theGeometry);
00379 iSetup.get<IdealMagneticFieldRecord>().get(theMagField);
00380
00381
00382
00383
00384
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
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;
00399 const reco::TrackRef tkref = itass->val;
00400
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
00409 if (stripFrontInvalidHits_) {
00410 while ( (begin != end) && ( (*begin)->isValid() == false ) ) ++begin;
00411 }
00412
00413 if (stripBackInvalidHits_ && (begin != end)) {
00414 --end;
00415 while ( (begin != end) && ( (*end)->isValid() == false ) ) --end;
00416 ++end;
00417 }
00418
00419
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{
00431 if ((end - begin) >= int(minimumHits_)) {
00432 output->push_back( makeCandidate ( *trk, begin, end ) );
00433 }
00434 }
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444 for (begin = hits.begin(), end = hits.end(); begin != end; ++begin) {
00445 if (*begin) delete *begin;
00446 }
00447 hits.clear();
00448 }
00449
00450
00451 }
00452 else{
00453
00454
00455 for (std::vector<reco::Track>::const_iterator ittrk = tracks->begin(), edtrk = tracks->end(); ittrk != edtrk; ++ittrk) {
00456
00457
00458
00459 const Track *trk = &(*ittrk);
00460
00461 produceFromTrack(iSetup,trk,hits);
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485 std::vector<TrackingRecHit *>::iterator begin = hits.begin(), end = hits.end();
00486
00487
00488 if (stripFrontInvalidHits_) {
00489 while ( (begin != end) && ( (*begin)->isValid() == false ) ) ++begin;
00490 }
00491
00492 if (stripBackInvalidHits_ && (begin != end)) {
00493 --end;
00494 while ( (begin != end) && ( (*end)->isValid() == false ) ) --end;
00495 ++end;
00496 }
00497
00498
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{
00510 if ((end - begin) >= int(minimumHits_)) {
00511 output->push_back( makeCandidate ( *ittrk, begin, end ) );
00512 }
00513 }
00514
00515
00516 for (begin = hits.begin(), end = hits.end(); begin != end; ++begin) {
00517 if (*begin) delete *begin;
00518 }
00519 hits.clear();
00520 }
00521 }
00522
00523
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
00537
00538
00539 if ( (pdir == alongMomentum) == ( (tk.outerPosition()-tk.innerPosition()).Dot(tk.momentum()) >= 0 ) ) {
00540
00541 TrajectoryStateOnSurface originalTsosIn(transform.innerStateOnSurface(tk, *theGeometry, &*theMagField));
00542 state = transform.persistentState( originalTsosIn, DetId(tk.innerDetId()) );
00543 } else {
00544
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
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
00565 hits.clear();
00566
00567 for (trackingRecHit_iterator ith = itt->recHitsBegin(), edh = itt->recHitsEnd(); ith != edh; ++ith) {
00568 const TrackingRecHit * hit = ith->get();
00569
00570 DetId detid = hit->geographicalId();
00571
00572
00573 if(hit->isValid() && hit==0 && detid.rawId()==0) continue;
00574
00575 int verdict=checkHit(iSetup,detid,hit);
00576 if (verdict == 0) {
00577
00578 hits.push_back(hit->clone());
00579 }
00580 else if(verdict<-2){
00581
00582 if (replaceWithInactiveHits_) {
00583 hits.push_back(new InvalidTrackingRecHit(detid, TrackingRecHit::inactive));
00584 }
00585 }
00586 else if(verdict==-2) hits.push_back(hit->clone());
00587 else if(verdict==-1){
00588 if (!stripAllInvalidHits_) {
00589 hits.push_back(hit->clone());
00590 }
00591 }
00592 }
00593
00594 }
00595
00596
00597 void TrackerTrackHitFilter::produceFromTrajectory(const edm::EventSetup &iSetup, const Trajectory *itt, std::vector<TrackingRecHit *>&hits){
00598 hits.clear();
00599 nOverlaps=0;
00600
00601
00602 std::vector<TrajectoryMeasurement> tmColl =itt->measurements();
00603
00604
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
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) ){
00624 verdict=-6;
00625 }
00626 }
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638 if(verdict==0){
00639
00640 if(tagOverlaps_){
00641
00642
00643 int layer(layerFromId(detid));
00644 int subDet = detid.subdetId();
00645
00646
00647 if ( ( previousTM!=0 )&& (layer!=-1 )) {
00648
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 {
00658
00659
00660
00661
00663
00664
00665
00666
00667
00668 nOverlaps++;
00669 break;
00670 }
00671 }
00672
00673 }
00674
00675 previousTM = &(* itTrajMeas);
00676 previousId = detid;
00677 previousLayer = layer;
00678 }
00680
00681 hits.push_back(hit->clone());
00682 }
00683 else if(verdict<-2){
00684
00685 if (replaceWithInactiveHits_) {
00686 hits.push_back(new InvalidTrackingRecHit(detid, TrackingRecHit::inactive));
00687 }
00688 }
00689 else if(verdict==-2) hits.push_back(hit->clone());
00690 else if(verdict==-1){
00691 if (!stripAllInvalidHits_) {
00692 hits.push_back(hit->clone());
00693 }
00694 }
00695 }
00696
00697 std::reverse(hits.begin(),hits.end());
00698
00699 }
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) {
00708 bool verdict = true;
00709
00710 for (std::vector<Rule>::const_iterator itr = rules_.begin(), edr = rules_.end(); itr != edr; ++itr) {
00711 itr->apply(detid, verdict);
00712 }
00713
00714
00715 if ( verdict){
00716 if(std::binary_search(detsToIgnore_.begin(), detsToIgnore_.end(), detid.rawId())) {
00717 hitresult=-4;
00718 }
00719 }
00720 else hitresult=-3;
00721
00722 if( hitresult==0 && rejectBadStoNHits_){
00723 if(!checkStoN(iSetup, detid, hit))hitresult=-5;
00724 }
00725 }
00726 else hitresult =-2;
00727 }
00728 else hitresult = -1;
00729 return hitresult;
00730 }
00731
00732
00733
00734 bool TrackerTrackHitFilter::checkStoN(const edm::EventSetup &iSetup, const DetId &id, const TrackingRecHit* therechit){
00735
00736 bool keepthishit=true;
00737
00738
00739
00740
00741 int subdet_cnt=id.subdetId();
00742
00743
00744
00745
00746
00747
00748 if(subdet_cnt>2){
00749 if( subdetStoN_[subdet_cnt-1]){
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
00769
00770
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
00781 }
00782
00783 }
00784
00785 }
00786 else if (subdet_cnt<=2){
00787
00788
00789 keepthishit = true;
00790
00791
00792
00793
00794
00795
00796 if(checkPXLQuality_){
00797 const SiPixelRecHit* pixelhit = dynamic_cast<const SiPixelRecHit*>(therechit);
00798 if(pixelhit!=0){
00799
00800 float xyprob=pixelhit->clusterProbability(0);
00801
00802 float xychargeprob =pixelhit->clusterProbability(1);
00803
00804 bool haspassed_tplreco= pixelhit->hasFilledProb();
00805 int qbin =pixelhit->qBin();
00806
00807
00808
00809
00810
00811 keepthishit = false;
00812
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 }
00818 }
00819
00820
00821
00822
00823
00824 return keepthishit;
00825 }
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
00836
00837
00838
00839
00840
00841 if(tsos.isValid()){
00842
00843 float mom_x=tsos.localDirection().x();
00844 float mom_y=tsos.localDirection().y();
00845 float mom_z=tsos.localDirection().z();
00846
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;
00849
00850
00851 if(angle_ok && PXLcorrClusChargeCut_>0.0){
00852
00853
00854 TransientTrackingRecHit::ConstRecHitPointer hitpointer = meas.recHit();
00855 if(hitpointer->isValid()){
00856 const TrackingRecHit *hit=(*hitpointer).hit();
00857 if( (hit->geographicalId()).subdetId()<=2 ){
00858 corrcharge_ok=false;
00859 float clust_alpha= atan2( mom_z, mom_x );
00860 float clust_beta= atan2( mom_z, mom_y );
00861
00862
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
00871 if(corr_clust_charge>PXLcorrClusChargeCut_){
00872 corrcharge_ok=true;
00873 }
00874 }
00875 }
00876
00877 }
00878
00879 }
00880 else{
00881 edm::LogWarning("TrackerTrackHitFilter") <<"TSOS not valid ! Impossible to calculate track angle.";
00882 }
00883
00884 return angle_ok&&corrcharge_ok;
00885 }
00886
00887
00888 bool TrackerTrackHitFilter::checkPXLCorrClustCharge(const TrajectoryMeasurement &meas){
00889
00890
00891
00892
00893 bool corrcharge_ok=false;
00894
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 ){
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
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 }
00921 return corrcharge_ok;
00922
00923 }
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 }}
00958
00959
00960
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);