Reference Manual

Go to the documentation of this file.
00001 #include "Calibration/TkAlCaRecoProducers/interface/CalibrationTrackSelector.h"
00002 
00003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00005 
00006 #include "FWCore/Framework/interface/Event.h"
00007 
00008 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
00009 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
00010 #include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
00011 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
00012 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
00013 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2DCollection.h"
00014 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2DCollection.h"
00015 
00016 #include "DataFormats/DetId/interface/DetId.h"
00017 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
00018 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
00019 
00020 const int kBPIX = PixelSubdetector::PixelBarrel;
00021 const int kFPIX = PixelSubdetector::PixelEndcap;
00022 
00023 // constructor ----------------------------------------------------------------
00024 
00025 CalibrationTrackSelector::CalibrationTrackSelector(const edm::ParameterSet & cfg) :
00026   applyBasicCuts_( cfg.getParameter<bool>( "applyBasicCuts" ) ),
00027   applyNHighestPt_( cfg.getParameter<bool>( "applyNHighestPt" ) ),
00028   applyMultiplicityFilter_( cfg.getParameter<bool>( "applyMultiplicityFilter" ) ),
00029   seedOnlyFromAbove_( cfg.getParameter<int>( "seedOnlyFrom" ) ),
00030   applyIsolation_( cfg.getParameter<bool>( "applyIsolationCut" ) ),
00031   chargeCheck_( cfg.getParameter<bool>( "applyChargeCheck" ) ),
00032   nHighestPt_( cfg.getParameter<int>( "nHighestPt" ) ),
00033   minMultiplicity_ ( cfg.getParameter<int>( "minMultiplicity" ) ),
00034   maxMultiplicity_ ( cfg.getParameter<int>( "maxMultiplicity" ) ),
00035   multiplicityOnInput_ ( cfg.getParameter<bool>( "multiplicityOnInput" ) ),
00036   ptMin_( cfg.getParameter<double>( "ptMin" ) ),
00037   ptMax_( cfg.getParameter<double>( "ptMax" ) ),
00038   etaMin_( cfg.getParameter<double>( "etaMin" ) ),
00039   etaMax_( cfg.getParameter<double>( "etaMax" ) ),
00040   phiMin_( cfg.getParameter<double>( "phiMin" ) ),
00041   phiMax_( cfg.getParameter<double>( "phiMax" ) ),
00042   nHitMin_( cfg.getParameter<double>( "nHitMin" ) ),
00043   nHitMax_( cfg.getParameter<double>( "nHitMax" ) ),
00044   chi2nMax_( cfg.getParameter<double>( "chi2nMax" ) ),
00045   minHitChargeStrip_( cfg.getParameter<double>( "minHitChargeStrip" ) ),
00046   minHitIsolation_( cfg.getParameter<double>( "minHitIsolation" ) ),
00047   rphirecHitsTag_( cfg.getParameter<edm::InputTag>("rphirecHits") ),
00048   matchedrecHitsTag_( cfg.getParameter<edm::InputTag>("matchedrecHits") ),
00049   nHitMin2D_( cfg.getParameter<unsigned int>( "nHitMin2D" ) ),
00050   // Ugly to use the same getParameter 6 times, but this allows const cut variables...
00051   minHitsinTIB_(cfg.getParameter<edm::ParameterSet>( "minHitsPerSubDet" ).getParameter<int>( "inTIB" ) ),
00052   minHitsinTOB_ (cfg.getParameter<edm::ParameterSet>( "minHitsPerSubDet" ).getParameter<int>( "inTOB" ) ),
00053   minHitsinTID_ (cfg.getParameter<edm::ParameterSet>( "minHitsPerSubDet" ).getParameter<int>( "inTID" ) ),
00054   minHitsinTEC_ (cfg.getParameter<edm::ParameterSet>( "minHitsPerSubDet" ).getParameter<int>( "inTEC" ) ),
00055   minHitsinBPIX_ (cfg.getParameter<edm::ParameterSet>( "minHitsPerSubDet" ).getParameter<int>( "inBPIX" ) ),
00056   minHitsinFPIX_ (cfg.getParameter<edm::ParameterSet>( "minHitsPerSubDet" ).getParameter<int>( "inFPIX" ) )
00057 {
00058 
00059   if (applyBasicCuts_)
00060         edm::LogInfo("CalibrationTrackSelector") 
00061           << "applying basic track cuts ..."
00062           << "\nptmin,ptmax:     " << ptMin_   << "," << ptMax_ 
00063           << "\netamin,etamax:   " << etaMin_  << "," << etaMax_
00064           << "\nphimin,phimax:   " << phiMin_  << "," << phiMax_
00065           << "\nnhitmin,nhitmax: " << nHitMin_ << "," << nHitMax_
00066           << "\nnhitmin2D:       " << nHitMin2D_
00067           << "\nchi2nmax:        " << chi2nMax_;
00068 
00069   if (applyNHighestPt_)
00070         edm::LogInfo("CalibrationTrackSelector") 
00071           << "filter N tracks with highest Pt N=" << nHighestPt_;
00072 
00073   if (applyMultiplicityFilter_)
00074         edm::LogInfo("CalibrationTrackSelector") 
00075           << "apply multiplicity filter N>= " << minMultiplicity_ << "and N<= " << maxMultiplicity_
00076           << " on " << (multiplicityOnInput_ ? "input" : "output");
00077 
00078   if (applyIsolation_)
00079     edm::LogInfo("CalibrationTrackSelector") 
00080       << "only retain tracks isolated at least by " << minHitIsolation_ << 
00081       " cm from other rechits"; 
00082 
00083   if (chargeCheck_)
00084     edm::LogInfo("CalibrationTrackSelector") 
00085       << "only retain hits with at least " << minHitChargeStrip_ <<  
00086       " ADC counts of total cluster charge"; 
00087 
00088   edm::LogInfo("CalibrationTrackSelector") 
00089     << "Minimum number of hits in TIB/TID/TOB/TEC/BPIX/FPIX = " 
00090     << minHitsinTIB_ << "/" << minHitsinTID_ << "/" << minHitsinTOB_
00091     << "/" << minHitsinTEC_ << "/" << minHitsinBPIX_ << "/" << minHitsinFPIX_;
00092 
00093 }
00094 
00095 // destructor -----------------------------------------------------------------
00096 
00097 CalibrationTrackSelector::~CalibrationTrackSelector()
00098 {}
00099 
00100 
00101 // do selection ---------------------------------------------------------------
00102 
00103 CalibrationTrackSelector::Tracks 
00104 CalibrationTrackSelector::select(const Tracks& tracks, const edm::Event& evt) const 
00105 {
00106   
00107   if (applyMultiplicityFilter_ && multiplicityOnInput_ && 
00108       (tracks.size() < static_cast<unsigned int>(minMultiplicity_) 
00109        || tracks.size() > static_cast<unsigned int>(maxMultiplicity_))) {
00110     return Tracks(); // empty collection
00111   }
00112   
00113   Tracks result = tracks;
00114   // apply basic track cuts (if selected)
00115   if (applyBasicCuts_) result= this->basicCuts(result, evt);
00116   
00117   // filter N tracks with highest Pt (if selected)
00118   if (applyNHighestPt_) result = this->theNHighestPtTracks(result);
00119   
00120   // apply minimum multiplicity requirement (if selected)
00121   if (applyMultiplicityFilter_ && !multiplicityOnInput_) {
00122     if (result.size() < static_cast<unsigned int>(minMultiplicity_) 
00123         || result.size() > static_cast<unsigned int>(maxMultiplicity_) ) {
00124 
00125       result.clear();
00126     }
00127   }
00128   
00129   // edm::LogDebug("CalibrationTrackSelector") << "tracks all,kept: " << tracks.size() << "," << result.size();
00130   
00131   return result;
00132 }
00133 
00134 // make basic cuts ------------------------------------------------------------
00135 
00136 CalibrationTrackSelector::Tracks 
00137 CalibrationTrackSelector::basicCuts(const Tracks& tracks, const edm::Event& evt) const 
00138 {
00139   Tracks result;
00140 
00141   for (Tracks::const_iterator it=tracks.begin(); it != tracks.end(); ++it) {
00142     const reco::Track* trackp=*it;
00143     float pt=trackp->pt();
00144     float eta=trackp->eta();
00145     float phi=trackp->phi();
00146     int nhit = trackp->numberOfValidHits(); 
00147     float chi2n = trackp->normalizedChi2();
00148 
00149     // edm::LogDebug("CalibrationTrackSelector") << " pt,eta,phi,nhit: "
00150     //  <<pt<<","<<eta<<","<<phi<<","<<nhit;
00151 
00152     if (pt>ptMin_ && pt<ptMax_ 
00153        && eta>etaMin_ && eta<etaMax_ 
00154        && phi>phiMin_ && phi<phiMax_ 
00155        && nhit>=nHitMin_ && nhit<=nHitMax_
00156        && chi2n<chi2nMax_) {
00157       if (this->detailedHitsCheck(trackp, evt)) result.push_back(trackp);
00158     }
00159   }
00160 
00161   return result;
00162 }
00163 
00164 //-----------------------------------------------------------------------------
00165 
00166 bool CalibrationTrackSelector::detailedHitsCheck(const reco::Track *trackp, const edm::Event& evt) const
00167 {
00168   // checking hit requirements beyond simple number of valid hits
00169 
00170   if (minHitsinTIB_ || minHitsinTOB_ || minHitsinTID_ || minHitsinTEC_
00171       || minHitsinFPIX_ || minHitsinBPIX_ || nHitMin2D_ || chargeCheck_
00172       || applyIsolation_) { // any detailed hit cut is active, so have to check
00173     
00174     int nhitinTIB = 0, nhitinTOB = 0, nhitinTID = 0;
00175     int nhitinTEC = 0, nhitinBPIX = 0, nhitinFPIX = 0;
00176     unsigned int nHit2D = 0;
00177     unsigned int thishit = 0;
00178 
00179     for (trackingRecHit_iterator iHit = trackp->recHitsBegin(); iHit != trackp->recHitsEnd(); ++iHit) {
00180 
00181       thishit++;
00182       int type = (*iHit)->geographicalId().subdetId(); 
00183 
00184       // *** thishit == 1 means last hit in CTF *** 
00185       // (FIXME: assumption might change or not be valid for all tracking algorthms)
00186       // ==> for cosmics
00187       // seedOnlyFrom = 1 is TIB-TOB-TEC tracks only
00188       // seedOnlyFrom = 2 is TOB-TEC tracks only 
00189       if (seedOnlyFromAbove_ == 1 && thishit == 1 && (type == int(StripSubdetector::TOB) || type == int(StripSubdetector::TEC))) return false; 
00190 
00191       if (seedOnlyFromAbove_ == 2 && thishit == 1 && type == int(StripSubdetector::TIB)) return false;  
00192 
00193       if (!(*iHit)->isValid()) continue; // only real hits count as in trackp->numberOfValidHits()
00194       const DetId detId((*iHit)->geographicalId());
00195       if (detId.det() != DetId::Tracker) {
00196         edm::LogError("DetectorMismatch") << "@SUB=CalibrationTrackSelector::detailedHitsCheck"
00197                                           << "DetId.det() != DetId::Tracker (=" << DetId::Tracker
00198                                           << "), but " << detId.det() << ".";
00199       }
00200       const TrackingRecHit* therechit = (*iHit).get();
00201       if (chargeCheck_ && !(this->isOkCharge(therechit))) return false;
00202       if (applyIsolation_ && (!this->isIsolated(therechit, evt))) return false;
00203       if      (StripSubdetector::TIB == detId.subdetId()) ++nhitinTIB;
00204       else if (StripSubdetector::TOB == detId.subdetId()) ++nhitinTOB;
00205       else if (StripSubdetector::TID == detId.subdetId()) ++nhitinTID;
00206       else if (StripSubdetector::TEC == detId.subdetId()) ++nhitinTEC;
00207       else if (                kBPIX == detId.subdetId()) ++nhitinBPIX;
00208       else if (                kFPIX == detId.subdetId()) ++nhitinFPIX;
00209       // Do not call isHit2D(..) if already enough 2D hits for performance reason:
00210       if (nHit2D < nHitMin2D_ && this->isHit2D(**iHit)) ++nHit2D;
00211     } // end loop on hits
00212     return (nhitinTIB >= minHitsinTIB_ && nhitinTOB >= minHitsinTOB_ 
00213             && nhitinTID >= minHitsinTID_ && nhitinTEC >= minHitsinTEC_ 
00214             && nhitinBPIX >= minHitsinBPIX_ && nhitinFPIX >= minHitsinFPIX_ 
00215             && nHit2D >= nHitMin2D_);
00216   } else { // no cuts set, so we are just fine and can avoid loop on hits
00217     return true;
00218   }
00219 }
00220 
00221 //-----------------------------------------------------------------------------
00222 
00223 bool CalibrationTrackSelector::isHit2D(const TrackingRecHit &hit) const
00224 {
00225   if (hit.dimension() < 2) {
00226     return false; // some (muon...) stuff really has RecHit1D
00227   } else {
00228     const DetId detId(hit.geographicalId());
00229     if (detId.det() == DetId::Tracker) {
00230       if (detId.subdetId() == kBPIX || detId.subdetId() == kFPIX) {
00231         return true; // pixel is always 2D
00232       } else { // should be SiStrip now
00233         if (dynamic_cast<const SiStripRecHit2D*>(&hit)) return false; // normal hit
00234         else if (dynamic_cast<const SiStripMatchedRecHit2D*>(&hit)) return true; // matched is 2D
00235         else if (dynamic_cast<const ProjectedSiStripRecHit2D*>(&hit)) return false; // crazy hit...
00236         else {
00237           edm::LogError("UnkownType") << "@SUB=CalibrationTrackSelector::isHit2D"
00238                                       << "Tracker hit not in pixel and neither SiStripRecHit2D nor "
00239                                       << "SiStripMatchedRecHit2D nor ProjectedSiStripRecHit2D.";
00240           return false;
00241         }
00242       }
00243     } else { // not tracker??
00244       edm::LogWarning("DetectorMismatch") << "@SUB=CalibrationTrackSelector::isHit2D"
00245                                           << "Hit not in tracker with 'official' dimension >=2.";
00246       return true; // dimension() >= 2 so accept that...
00247     }
00248   }
00249   // never reached...
00250 }
00251 
00252 //-----------------------------------------------------------------------------
00253 
00254 bool CalibrationTrackSelector::isOkCharge(const TrackingRecHit* therechit) const
00255 {
00256   // const TrackingRecHit* therechit = hit.get();
00257   float charge1 = 0;
00258   float charge2 = 0;
00259   const SiStripMatchedRecHit2D* matchedhit = dynamic_cast<const SiStripMatchedRecHit2D*>(therechit);
00260   const SiStripRecHit2D* hit = dynamic_cast<const SiStripRecHit2D*>(therechit);
00261   const ProjectedSiStripRecHit2D* unmatchedhit = dynamic_cast<const ProjectedSiStripRecHit2D*>(therechit);
00262   
00263   if (matchedhit) {  
00264     const SiStripCluster & monocluster = matchedhit->monoCluster();
00265     const std::vector<uint16_t> amplitudesmono( monocluster.amplitudes().begin(),
00266                                                 monocluster.amplitudes().end());
00267     for(size_t ia=0; ia<amplitudesmono.size();++ia)
00268       { charge1+=amplitudesmono[ia];} 
00269     
00270     const SiStripCluster & stereocluster = matchedhit->stereoCluster();
00271     const std::vector<uint16_t> amplitudesstereo( stereocluster.amplitudes().begin(),
00272                                                   stereocluster.amplitudes().end());
00273     for(size_t ia=0; ia<amplitudesstereo.size();++ia)
00274       {charge2+=amplitudesstereo[ia];}
00275     // std::cout << "charge1 = " << charge1 << "\n";
00276     // std::cout << "charge2 = " << charge2 << "\n";
00277     if (charge1 < minHitChargeStrip_ || charge2 < minHitChargeStrip_) return false;
00278   }
00279   else if (hit) {
00280     
00281     const SiStripCluster* cluster = &*(hit->cluster());
00282     const std::vector<uint16_t> amplitudes( cluster->amplitudes().begin(),
00283                                             cluster->amplitudes().end());
00284     for(size_t ia=0; ia<amplitudes.size();++ia)
00285       {charge1+=amplitudes[ia];}
00286     // std::cout << "charge1 = " << charge1 << "\n";
00287     if (charge1 < minHitChargeStrip_) return false;
00288   }
00289   else if (unmatchedhit) {
00290     
00291     const SiStripRecHit2D &orighit = unmatchedhit->originalHit(); 
00292     const SiStripCluster* origcluster = &*(orighit.cluster());
00293     const std::vector<uint16_t> amplitudes( origcluster->amplitudes().begin(),
00294                                             origcluster->amplitudes().end());
00295     for(size_t ia=0; ia<amplitudes.size();++ia)
00296       {charge1+=amplitudes[ia];}
00297     // std::cout << "charge1 = " << charge1 << "\n";
00298     if (charge1 < minHitChargeStrip_) return false;
00299   } 
00300   return true;
00301 } 
00302 
00303 //-----------------------------------------------------------------------------
00304 
00305 bool CalibrationTrackSelector::isIsolated(const TrackingRecHit* therechit, const edm::Event& evt) const
00306 {
00307            
00308   // edm::ESHandle<TrackerGeometry> tracker;
00309   edm::Handle<SiStripRecHit2DCollection> rphirecHits;
00310   edm::Handle<SiStripMatchedRecHit2DCollection> matchedrecHits;
00311   // es.get<TrackerDigiGeometryRecord>().get(tracker);
00312   evt.getByLabel( rphirecHitsTag_, rphirecHits );
00313   evt.getByLabel( matchedrecHitsTag_, matchedrecHits ); 
00314 
00315   SiStripRecHit2DCollection::DataContainer::const_iterator istripSt; 
00316   SiStripMatchedRecHit2DCollection::DataContainer::const_iterator istripStm; 
00317   const SiStripRecHit2DCollection& stripcollSt = *rphirecHits;
00318   const SiStripMatchedRecHit2DCollection& stripcollStm = *matchedrecHits;
00319   
00320   DetId idet = therechit->geographicalId(); 
00321   
00322   // FIXME: instead of looping the full hit collection, we should explore the features of 
00323   // SiStripRecHit2DCollection::rangeRphi = rphirecHits.get(idet) and loop
00324   // only from rangeRphi.first until rangeRphi.second
00325   for( istripSt=stripcollSt.data().begin(); istripSt!=stripcollSt.data().end(); ++istripSt ) {
00326     const SiStripRecHit2D *aHit = &*(istripSt);
00327     DetId mydet1 = aHit->geographicalId(); 
00328     if (idet.rawId() != mydet1.rawId()) continue; 
00329     float theDistance = ( therechit->localPosition() - aHit->localPosition() ).mag();
00330     // std::cout << "theDistance1 = " << theDistance << "\n";
00331     if (theDistance > 0.001 && theDistance < minHitIsolation_) return false;
00332   }
00333   
00334   // FIXME: see above
00335   for( istripStm=stripcollStm.data().begin(); istripStm!=stripcollStm.data().end(); ++istripStm ) {
00336     const SiStripMatchedRecHit2D *aHit = &*(istripStm);
00337     DetId mydet2 = aHit->geographicalId(); 
00338     if (idet.rawId() != mydet2.rawId()) continue;
00339     float theDistance = (therechit->localPosition() - aHit->localPosition()).mag();
00340     // std::cout << "theDistance1 = " << theDistance << "\n";
00341     if (theDistance > 0.001 && theDistance < minHitIsolation_) return false;
00342   }
00343   return true;
00344 }
00345 //-----------------------------------------------------------------------------
00346 
00347 CalibrationTrackSelector::Tracks 
00348 CalibrationTrackSelector::theNHighestPtTracks(const Tracks& tracks) const
00349 {
00350   Tracks sortedTracks=tracks;
00351   Tracks result;
00352 
00353   // sort in pt
00354   std::sort(sortedTracks.begin(),sortedTracks.end(),ptComparator);
00355 
00356   // copy theTrackMult highest pt tracks to result vector
00357   int n=0;
00358   for (Tracks::const_iterator it=sortedTracks.begin();
00359            it!=sortedTracks.end(); ++it) {
00360         if (n<nHighestPt_) { result.push_back(*it); n++; }
00361   }
00362 
00363   return result;
00364 }
00365
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/Calibration/TkAlCaRecoProducers/src/CalibrationTrackSelector.cc
