00001 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h" 00002 #include "RecoVertex/PrimaryVertexProducer/interface/GapClusterizerInZ.h" 00003 #include "RecoVertex/VertexPrimitives/interface/VertexException.h" 00004 00005 00006 using namespace std; 00007 00008 00009 namespace { 00010 00011 bool recTrackLessZ(const reco::TransientTrack & tk1, 00012 const reco::TransientTrack & tk2) 00013 { 00014 return tk1.stateAtBeamLine().trackStateAtPCA().position().z() < tk2.stateAtBeamLine().trackStateAtPCA().position().z(); 00015 } 00016 00017 } 00018 00019 00020 00021 GapClusterizerInZ::GapClusterizerInZ(const edm::ParameterSet& conf) 00022 { 00023 // some defaults to avoid uninitialized variables 00024 verbose_= conf.getUntrackedParameter<bool>("verbose", false); 00025 zSep = conf.getParameter<double>("zSeparation"); 00026 if(verbose_) {std::cout << "TrackClusterizerInZ: algorithm=gap, zSeparation="<< zSep << std::endl;} 00027 } 00028 00029 00030 00031 float GapClusterizerInZ::zSeparation() const 00032 { 00033 return zSep; 00034 } 00035 00036 00037 00038 00039 vector< vector<reco::TransientTrack> > 00040 GapClusterizerInZ::clusterize(const vector<reco::TransientTrack> & tracks) 00041 const 00042 { 00043 00044 vector<reco::TransientTrack> tks = tracks; // copy to be sorted 00045 00046 vector< vector<reco::TransientTrack> > clusters; 00047 if (tks.empty()) return clusters; 00048 00049 // sort in increasing order of z 00050 stable_sort(tks.begin(), tks.end(), recTrackLessZ); 00051 00052 // init first cluster 00053 vector<reco::TransientTrack>::const_iterator it = tks.begin(); 00054 vector <reco::TransientTrack> currentCluster; currentCluster.push_back(*it); 00055 00056 it++; 00057 for ( ; it != tks.end(); it++) { 00058 00059 double zPrev = currentCluster.back().stateAtBeamLine().trackStateAtPCA().position().z(); 00060 double zCurr = (*it).stateAtBeamLine().trackStateAtPCA().position().z(); 00061 00062 if ( abs(zCurr - zPrev) < zSeparation() ) { 00063 // close enough ? cluster together 00064 currentCluster.push_back(*it); 00065 } 00066 else { 00067 // store current cluster, start new one 00068 clusters.push_back(currentCluster); 00069 currentCluster.clear(); 00070 currentCluster.push_back(*it); 00071 } 00072 } 00073 00074 // store last cluster 00075 clusters.push_back(currentCluster); 00076 00077 return clusters; 00078 00079 } 00080 00081