Reference Manual

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    TrackHitFilter
00004 // Class:      TrackHitFilter
00005 // 
00013 //
00014 // Original Author:  Roberto Covarelli
00015 //         Created:  Mon Jan 15 10:39:42 CET 2007
00016 // $Id: TrackHitFilter.cc,v 1.15 2010/03/31 20:40:55 mussgill Exp $
00017 //
00018 //
00019 
00020 #include "Alignment/TrackHitFilter/interface/TrackHitFilter.h"
00021 #include "FWCore/Framework/interface/Event.h" 
00022 #include "FWCore/Framework/interface/MakerMacros.h" 
00023 #include "FWCore/ParameterSet/interface/ParameterSet.h" 
00024 #include "FWCore/MessageLogger/interface/MessageLogger.h" 
00025 #include "DataFormats/TrackReco/interface/Track.h" 
00026 #include "DataFormats/TrackReco/interface/TrackExtra.h" 
00027 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h" 
00028 #include "Alignment/TrackerAlignment/interface/TrackerAlignableId.h" 
00029 #include "DataFormats/TrajectorySeed/interface/PropagationDirection.h" 
00030 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00031 
00032 #include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
00033 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
00034 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
00035 #include "AnalysisDataFormats/SiStripClusterInfo/interface/SiStripClusterInfo.h"
00036 
00037 #include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
00038 
00039 // RC
00040 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
00041 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
00042 
00043 //
00044 // constructors and destructor
00045 //
00046 TrackHitFilter::TrackHitFilter(const edm::ParameterSet& iConfig):
00047   theSrc( iConfig.getParameter<edm::InputTag>( "src" ) ),
00048   theHitSel( iConfig.getParameter<std::string>( "hitSelection" ) ),
00049   theMinHits( iConfig.getParameter<unsigned int>( "minHitsForRefit" ) ),
00050   rejectBadMods( iConfig.getParameter<bool>( "rejectBadMods" ) ),
00051   theBadMods( iConfig.getParameter<std::vector<unsigned int> >( "theBadModules" ) ),
00052   rejectBadStoNHits( iConfig.getParameter<bool>( "rejectBadStoNHits" ) ),
00053   theCMNSubtractionMode(iConfig.getUntrackedParameter<std::string>( "CMNSubtractionMode" ,"Median") ) ,
00054   theStoNthreshold( iConfig.getParameter<double>( "theStoNthreshold" ) ),  
00055   rejectBadClusterPixelHits( iConfig.getParameter<bool>( "rejectBadClusterPixelHits" ) ),
00056   thePixelClusterthreshold( iConfig.getParameter<double>( "thePixelClusterthreshold" ) )
00057 {
00058 
00059    //register your products, and/or set an "alias" label
00060   produces<reco::TrackCollection>();
00061   produces<reco::TrackExtraCollection>();
00062   produces<TrackingRecHitCollection>();
00063  
00064 }
00065 
00066 
00067 TrackHitFilter::~TrackHitFilter()
00068 {
00069 }
00070 
00071 void TrackHitFilter::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
00072 {
00073 
00074    //Read track collection from the Event
00075    edm::Handle<reco::TrackCollection> trackAllHits;
00076    iEvent.getByLabel(theSrc,trackAllHits);
00077 
00078    // Create empty Track, TrackExtra and TrackingRecHits collections
00079    std::auto_ptr<reco::TrackCollection> trackSelectedHits( new reco::TrackCollection() );
00080    std::auto_ptr<reco::TrackExtraCollection> txSelectedHits( new reco::TrackExtraCollection() );
00081    std::auto_ptr<TrackingRecHitCollection> trhSelectedHits( new TrackingRecHitCollection() );
00082 
00083    //get references
00084    TrackingRecHitRefProd rHits = iEvent.getRefBeforePut<TrackingRecHitCollection>();
00085    reco::TrackExtraRefProd rTrackExtras = iEvent.getRefBeforePut<reco::TrackExtraCollection>();
00086    reco::TrackRefProd rTracks = iEvent.getRefBeforePut<reco::TrackCollection>();
00087    edm::Ref<reco::TrackExtraCollection>::key_type idx = 0;
00088    edm::Ref<reco::TrackExtraCollection>::key_type hidx = 0;
00089 
00090    TrackerAlignableId* TkMap = new TrackerAlignableId();
00091    
00092    LogDebug("HitFilter")<< trackAllHits->size() << " track(s) found in the event with label " << theSrc<<std::endl;
00093 
00094    unsigned int nTr = 0;
00095    std::vector<unsigned int> accHits;
00096    std::vector<unsigned int> allHits;
00097 
00098    // first iteration : count hits
00099    for( reco::TrackCollection::const_iterator iTrack = trackAllHits->begin();
00100         iTrack != trackAllHits->end();
00101         ++iTrack ) {
00102        
00103      const reco::Track * trk = &(*iTrack);
00104      
00105      unsigned int allhits = 0;
00106      unsigned int acchits = 0;
00107      for (trackingRecHit_iterator iHit = trk->recHitsBegin(); iHit != trk->recHitsEnd(); iHit++) {
00108        
00109        allhits++;
00110        TrackingRecHit* hit = (*iHit)->clone();
00111        //std::cout<<"See if it crashes. Geographical Id of the hit: "<<hit->geographicalId()<<std::endl;
00112        std::pair<int,int> typeAndLay = TkMap->typeAndLayerFromDetId( hit->geographicalId() );
00113        int type = typeAndLay.first;   
00114        int layer = typeAndLay.second;
00115        DetId detid= hit->geographicalId();
00116        if (keepThisHit( detid, type, layer, hit, iEvent,iSetup )) acchits++; 
00117        
00118      }//end loop on hits
00119      
00120      if (!nTr) {
00121        LogDebug("HitFilter") << "TrackHitFilter **** In first track " << acchits << " RecHits retained out of " << allhits;
00122        nTr++;
00123      }
00124      allHits.push_back(allhits);
00125      accHits.push_back(acchits);
00126      
00127    }//end loop on tracks
00128 
00129    nTr = 0;
00130    // second iteration : store tracks
00131    for( reco::TrackCollection::const_iterator iTrack2 = trackAllHits->begin();
00132         iTrack2 != trackAllHits->end();
00133         ++iTrack2 ) {
00134        
00135      if (accHits.at(nTr) >= theMinHits) {
00136        const reco::Track * trk = &(*iTrack2);
00137        reco::Track * myTrk = new reco::Track(*trk);
00138        PropagationDirection seedDir = trk->seedDirection();
00139        myTrk->setExtra( reco::TrackExtraRef( rTrackExtras, idx ++ ) );
00140        reco::TrackExtra * tx = new reco::TrackExtra( trk->outerPosition(), trk->outerMomentum(), 
00141                                                      trk->outerOk(), trk->innerPosition(), 
00142                                                      trk->innerMomentum(), trk->innerOk(),
00143                                                      trk->outerStateCovariance(), trk->outerDetId(),
00144                                                      trk->innerStateCovariance(), trk->innerDetId() ,
00145                                                      seedDir ) ;        
00146 
00147        unsigned int i = 0; 
00148        for (trackingRecHit_iterator iHit = trk->recHitsBegin(); iHit != trk->recHitsEnd(); iHit++) {
00149 
00150          TrackingRecHit* hit = (*iHit)->clone();
00151          std::pair<int,int> typeAndLay = TkMap->typeAndLayerFromDetId( hit->geographicalId() );
00152          int type = typeAndLay.first;   
00153          int layer = typeAndLay.second;
00154          
00155          if (keepThisHit( hit->geographicalId(), type, layer, hit, iEvent, iSetup )) {
00156    
00157            myTrk->setHitPattern( * hit, i ++ );
00158            trhSelectedHits->push_back( hit );
00159            tx->add( TrackingRecHitRef( rHits, hidx ++ ) );
00160 
00161          }
00162        } 
00163        
00164        trackSelectedHits->push_back( *myTrk );
00165        txSelectedHits->push_back( *tx );
00166      }
00167        
00168      nTr++;
00169    }
00170     
00171    iEvent.put( trackSelectedHits );     
00172    iEvent.put( txSelectedHits );
00173    iEvent.put( trhSelectedHits );
00174    //   std::cout<<"Finished this call to <TrackHitFilter::produce>"<<std::endl;
00175 }
00176 
00177 
00178 void 
00179 TrackHitFilter::beginJob()
00180 {
00181 }
00182 
00183 void 
00184 TrackHitFilter::endJob() {
00185 }
00186 
00187 bool TrackHitFilter::keepThisHit(DetId id, int type, int layer, const TrackingRecHit* therechit, const edm::Event &iEvent, const edm::EventSetup& iSetup) 
00188 {
00189   bool keepthishit = true;
00190   
00191   if ( theHitSel == "PixelOnly" ) {
00192     if (abs(type)<1 || abs(type)>2) keepthishit = false; 
00193   }
00194   else if ( theHitSel == "PixelBarrelOnly" ) {
00195     if (abs(type)!=1) keepthishit = false; 
00196   }
00197   else if ( theHitSel == "PixelAndDSStripBarrelOnly" ) {
00198     if (!((abs(type)==1)
00199           || ((abs(type)==3 || abs(type)==5) && layer<2))) keepthishit = false;
00200   }        
00201   else if ( theHitSel == "SiStripOnly" ) {
00202     if (abs(type)>=1 && abs(type)<=2) keepthishit = false; 
00203   } 
00204   else if ( theHitSel == "TOBOnly" ) {
00205      if (abs(type)!=5) keepthishit = false;
00206   }
00207   else if ( theHitSel == "TOBandTIBOnly" ) {
00208      if ( !(abs(type)==3 || abs(type)==5) )  keepthishit = false;
00209   }
00210   else if ( theHitSel == "TOBandTIBl4Only" ) {
00211     if (!(abs(type)==5 || (abs(type)==3 && layer>=4))) keepthishit = false;
00212   }
00213   else if ( theHitSel == "TOBandTIBl34Only" ) {
00214     if (!(abs(type)==5 || (abs(type)==3 && layer>=3))) keepthishit = false;
00215   }
00216   else if ( theHitSel == "TOBandTIBl234Only" ) {
00217     if (!(abs(type)==5 || (abs(type)==3 && layer>=2))) keepthishit = false;
00218   }
00219   else if ( theHitSel == "noPXE" ) {
00220     if (abs(type)==2) keepthishit = false;
00221   }
00222   
00223   if (rejectBadMods) {
00224     for( std::vector<unsigned int>::const_iterator iMod = theBadMods.begin(); iMod != theBadMods.end(); ++iMod ) {  
00225       if (id.rawId() == *iMod) {
00226         keepthishit = false;
00227         LogDebug("HitFilter") << "TrackHitFilter **** Rejected a hit on module " << id.rawId();
00228       }
00229     }
00230   }
00231 
00232 
00233   //********* Rejects pixel hits with 'bad' cluster 
00234 
00235   if(rejectBadClusterPixelHits && (abs(type)==1 || abs(type)==2)){
00236 
00237     const SiPixelRecHit* pixelhit = dynamic_cast<const SiPixelRecHit*>(therechit);
00238     const SiPixelCluster* pixelcluster = &*(pixelhit->cluster());
00239     if(pixelcluster->size()==1 && pixelcluster->charge()<thePixelClusterthreshold){keepthishit = false;}
00240   }
00241 
00242 
00243   // ****************        
00244   // Reject hits with bad S/N
00245   if (rejectBadStoNHits && (abs(type)>2) ) { //apply it only to Strip hits
00246     /*** RC ****/ 
00247     //   const uint32_t& recHitDetId = id.rawId();
00248     const SiStripMatchedRecHit2D* matchedhit = dynamic_cast<const SiStripMatchedRecHit2D*>(therechit);
00249     const SiStripRecHit2D* hit = dynamic_cast<const SiStripRecHit2D*>(therechit);
00250     const ProjectedSiStripRecHit2D* unmatchedhit = dynamic_cast<const ProjectedSiStripRecHit2D*>(therechit);
00251 
00252     if (matchedhit) {  
00253       bool keepmonohit=true;
00254       bool keepstereohit=true;
00255 
00256       const SiStripRecHit2D* monohit=matchedhit->monoHit();    
00257       const SiStripCluster* monocluster = &*(monohit->cluster());
00258       SiStripClusterInfo monoclusterInfo = SiStripClusterInfo( *monocluster, iSetup); 
00259       if (monoclusterInfo.signalOverNoise() < theStoNthreshold ) keepmonohit = false;
00260         
00261 
00262       const SiStripRecHit2D* stereohit=matchedhit->stereoHit();   
00263       const SiStripCluster* stereocluster = &*(stereohit->cluster());
00264       SiStripClusterInfo stereoclusterInfo = SiStripClusterInfo(*stereocluster, iSetup);   
00265       if (stereoclusterInfo.signalOverNoise() < theStoNthreshold )keepstereohit = false;
00266            
00267       if (!keepmonohit || !keepstereohit) keepthishit = false;    
00268     }
00269     else if (hit) {
00270       const SiStripCluster* cluster = &*(hit->cluster());
00271       SiStripClusterInfo clusterInfo = SiStripClusterInfo(*cluster, iSetup);     
00272        if (clusterInfo.signalOverNoise() < theStoNthreshold )keepthishit = false;
00273      
00274     }
00275     else if (unmatchedhit) {
00276       const SiStripRecHit2D &orighit = unmatchedhit->originalHit(); 
00277       const SiStripCluster* origcluster = &*(orighit.cluster());
00278       SiStripClusterInfo clusterInfo = SiStripClusterInfo(*origcluster, iSetup);     
00279 
00280       if (clusterInfo.signalOverNoise() < theStoNthreshold ) keepthishit = false;   
00281          
00282     }
00283   } // end reject bad S/N
00284 
00285   
00286   return keepthishit;
00287 }
00288 
00289 //define this as a plug-in
00290 DEFINE_FWK_MODULE(TrackHitFilter);
/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/Alignment/TrackHitFilter/src/TrackHitFilter.cc
