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
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
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
00144 typedef map<unsigned int,int,less<unsigned int> > TrackMap;
00145
00146 for(unsigned int i = 0; i < tracks.size(); i++)
00147 {
00148
00149 if(!keep[i]) continue;
00150
00151 TrackMap trackMap;
00152 vector<DetId> detIds;
00153 vector<int> detLayers;
00154
00155
00156 for(vector<const TrackingRecHit *>::const_iterator
00157 recHit = tracks[i].second.begin();
00158 recHit!= tracks[i].second.end(); recHit++)
00159 {
00160
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
00168 detIds.push_back((*recHit)->geographicalId());
00169 detLayers.push_back(getLayer((*recHit)->geographicalId()));
00170 }
00171
00172
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 {
00181 if((*sharing).second > min(tracks[i].second.size(),
00182 tracks[j].second.size())/2)
00183 {
00184
00185 if(!hasCommonLayer(tracks[i].second,tracks[j].second,detLayers))
00186 {
00187
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
00198 keep[j] = false;
00199
00200 changes++;
00201 }
00202 else
00203 {
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 {
00237 if((*sharing).second > 0)
00238 {
00239
00240 keep[j] = false;
00241
00242 changes++;
00243 }
00244 }
00245 }
00246 }
00247 }
00248 while(changes > 0);
00249
00250
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