CMS 3D CMS Logo

TrackCleaner.cc

Go to the documentation of this file.
00001 #include "RecoPixelVertexing/PixelLowPtUtilities/interface/TrackCleaner.h"
00002 #include "RecoPixelVertexing/PixelLowPtUtilities/interface/HitInfo.h"
00003 
00004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00005 
00006 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
00007 #include "DataFormats/TrackReco/interface/Track.h"
00008 
00009 #include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
00010 #include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
00011 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
00012 
00013 using namespace std;
00014 using namespace pixeltrackfitting;
00015 
00016 /*****************************************************************************/
00017 class HitComparator
00018 {
00019 public:
00020   bool operator() (const TrackingRecHit* a, const TrackingRecHit* b) const
00021   {
00022     if(a->geographicalId() < b->geographicalId()) return true;
00023     if(b->geographicalId() < a->geographicalId()) return false;
00024 
00025     if(a->localPosition().x() < b->localPosition().x() - 1e-5) return true;
00026     if(b->localPosition().x() < a->localPosition().x() - 1e-5) return false;
00027 
00028     if(a->localPosition().y() < b->localPosition().y() - 1e-5) return true;
00029     if(b->localPosition().y() < a->localPosition().y() - 1e-5) return false;
00030 
00031     return false;
00032   }
00033 };
00034 
00035 /*****************************************************************************/
00036 TrackCleaner::TrackCleaner
00037   (const edm::ParameterSet& ps)
00038 {
00039 }
00040 
00041 /*****************************************************************************/
00042 TrackCleaner::~TrackCleaner()
00043 {
00044 }
00045 
00046 /*****************************************************************************/
00047 int TrackCleaner::getLayer(const DetId & id)
00048 {
00049   if(id.subdetId() == int(PixelSubdetector::PixelBarrel))
00050   {
00051     PXBDetId pid(id);
00052     return 0 + ((pid.layer() - 1)*2) + ((pid.ladder() - 1)%2);
00053   } 
00054   else
00055   {
00056     PXFDetId pid(id);
00057     return 6 + ((pid.disk()  - 1)*2) + ((pid.panel()  - 1)%2);
00058   } 
00059 }
00060 
00061 /*****************************************************************************/
00062 bool TrackCleaner::hasCommonDetUnit
00063   (vector<const TrackingRecHit *> recHitsA,
00064    vector<const TrackingRecHit *> recHitsB,
00065    vector<DetId> detIds)
00066 {
00067   for(vector<const TrackingRecHit *>::const_iterator
00068       recHit = recHitsB.begin(); recHit!= recHitsB.end(); recHit++)
00069     if(find(recHitsA.begin(), recHitsA.end(), *recHit) == recHitsA.end())
00070     if(find(detIds.begin(),detIds.end(),
00071             (*recHit)->geographicalId()) != detIds.end())
00072       return true;
00073 
00074   return false;
00075 }
00076 
00077 /*****************************************************************************/
00078 bool TrackCleaner::hasCommonLayer
00079   (vector<const TrackingRecHit *> recHitsA,
00080    vector<const TrackingRecHit *> recHitsB,
00081    vector<int> detLayers)
00082 {
00083   for(vector<const TrackingRecHit *>::const_iterator
00084       recHit = recHitsB.begin(); recHit!= recHitsB.end(); recHit++)
00085     if(find(recHitsA.begin(), recHitsA.end(), *recHit) == recHitsA.end())
00086     if(find(detLayers.begin(),detLayers.end(),
00087             getLayer((*recHit)->geographicalId())) != detLayers.end())
00088       return true;
00089 
00090   return false;
00091 }
00092 
00093 /*****************************************************************************/
00094 struct RadiusComparator
00095 { 
00096   bool operator() (const TrackingRecHit * h1,
00097                    const TrackingRecHit * h2)
00098   { 
00099     return (h1 < h2);
00100   };
00101 };
00102 
00103 /*****************************************************************************/
00104 TracksWithRecHits TrackCleaner::cleanTracks
00105   (const TracksWithRecHits & tracks_)
00106 {
00107   LogTrace("MinBiasTracking")
00108     << " [TrackCleaner]";
00109 
00110   // Local copy
00111   TracksWithRecHits tracks = tracks_;
00112 
00113   typedef map<const TrackingRecHit*,vector<unsigned int>,HitComparator>
00114     RecHitMap;
00115 
00116   vector<bool> keep(tracks.size(),true);
00117 
00118   int changes;
00119 
00120   LogTrace("MinBiasTracking")
00121     << " [TrackCleaner ] initial tracks : " << tracks.size();
00122 
00123   for(unsigned int i = 0; i < tracks.size(); i++)
00124   LogTrace("MinBiasTracking")
00125     << "   Track #" << i << " : " << HitInfo::getInfo(tracks[i].second);
00126 
00127   do
00128   {
00129   changes = 0;
00130 
00131   RecHitMap recHitMap;
00132 
00133   // Fill the rechit map
00134   for(unsigned int i = 0; i < tracks.size(); i++)
00135   if(keep[i])
00136   {
00137     for(vector<const TrackingRecHit *>::const_iterator
00138         recHit = tracks[i].second.begin();
00139         recHit!= tracks[i].second.end(); recHit++)
00140       recHitMap[*recHit].push_back(i);
00141   }
00142 
00143   // Look at each track
00144   typedef map<unsigned int,int,less<unsigned int> > TrackMap; 
00145 
00146   for(unsigned int i = 0; i < tracks.size(); i++)
00147   {
00148     // Skip if 'i' already removed
00149     if(!keep[i]) continue;
00150 
00151     TrackMap trackMap;
00152     vector<DetId> detIds;
00153     vector<int> detLayers;
00154 
00155     // Go trough all rechits of this track
00156     for(vector<const TrackingRecHit *>::const_iterator
00157         recHit = tracks[i].second.begin();
00158         recHit!= tracks[i].second.end(); recHit++)
00159     {
00160       // Get tracks sharing this rechit
00161       vector<unsigned int> sharing = recHitMap[*recHit];
00162 
00163       for(vector<unsigned int>::iterator j = sharing.begin();
00164                                          j!= sharing.end(); j++)
00165         if(i < *j) trackMap[*j]++;
00166 
00167       // Fill detLayers vector
00168       detIds.push_back((*recHit)->geographicalId());
00169       detLayers.push_back(getLayer((*recHit)->geographicalId()));
00170     }
00171 
00172     // Check for tracks with shared rechits
00173     for(TrackMap::iterator sharing = trackMap.begin();
00174                            sharing!= trackMap.end(); sharing++)
00175     {
00176       unsigned int j = (*sharing).first;
00177       if(!keep[i] || !keep[j]) continue;
00178 
00179       if(tracks[i].second.size() >=3) 
00180       { // triplet tracks
00181         if((*sharing).second > min(tracks[i].second.size(),
00182                                    tracks[j].second.size())/2)
00183         { // more than min(hits1,hits2)/2 rechits are shared
00184 
00185           if(!hasCommonLayer(tracks[i].second,tracks[j].second,detLayers))
00186           { 
00187             // merge tracks, add separate hits of the second to the first one
00188             for(vector<const TrackingRecHit *>::const_iterator
00189                 recHit = tracks[j].second.begin();
00190                 recHit!= tracks[j].second.end(); recHit++)
00191               if(find(tracks[i].second.begin(),
00192                       tracks[i].second.end(),*recHit) == tracks[i].second.end())
00193                 tracks[i].second.push_back(*recHit);
00194 
00195             LogTrace("MinBiasTracking") << " merge " << i << " " << j;
00196   
00197             // Remove second track
00198             keep[j] = false;
00199   
00200            changes++;
00201           }
00202           else
00203           { // remove track with higher impact / chi2
00204             if(tracks[i].first->chi2() < tracks[j].first->chi2())
00205               keep[j] = false;
00206             else
00207               keep[i] = false;
00208 
00209             LogTrace("MinBiasTracking") << " merge + keep lower chi2 " << i << " " << j;
00210   
00211             changes++;
00212           }
00213         }
00214         else
00215         {
00216           if((*sharing).second > 1)
00217           {
00218             if(tracks[i].second.size() != tracks[j].second.size())
00219             {
00220               if(tracks[i].second.size() > tracks[j].second.size()) 
00221                 keep[j] = false; else keep[i] = false; 
00222               changes++;
00223             LogTrace("MinBiasTracking") << " sharing " << (*sharing).second << " remove by size";
00224             }
00225             else
00226             { 
00227               if(tracks[i].first->chi2() < tracks[j].first->chi2())
00228                 keep[j] = false; else keep[i] = false; 
00229               changes++;
00230             LogTrace("MinBiasTracking") << " sharing " << (*sharing).second << " remove by chi2";
00231             } 
00232           }
00233         }
00234       }
00235       else
00236       { // pair tracks
00237         if((*sharing).second > 0)
00238         {
00239           // Remove second track
00240           keep[j] = false;
00241   
00242           changes++;
00243         }
00244       }
00245     }
00246   }
00247   }
00248   while(changes > 0);
00249 
00250   // Final copy
00251   TracksWithRecHits cleaned;
00252   
00253   for(unsigned int i = 0; i < tracks.size(); i++)
00254     if(keep[i]) 
00255       cleaned.push_back(tracks[i]);
00256     else delete tracks_[i].first;
00257 
00258   LogTrace("MinBiasTracking")
00259     << " [TrackCleaner ] cleaned tracks : " << cleaned.size();
00260 
00261   for(unsigned int i = 0; i < cleaned.size(); i++)
00262   LogTrace("MinBiasTracking")
00263     << "   Track #" << i << " : " << HitInfo::getInfo(cleaned[i].second);
00264 
00265   return cleaned;
00266 }
00267 

Generated on Tue Jun 9 17:44:51 2009 for CMSSW by  doxygen 1.5.4