PrimaryVertexProducer/src/GapClusterizerInZ.cc

0001 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
0002 #include "RecoVertex/PrimaryVertexProducer/interface/GapClusterizerInZ.h"
0003 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
0004
0005 using namespace std;
0006
0007 namespace {
0008
0009   bool recTrackLessZ(const reco::TransientTrack& tk1, const reco::TransientTrack& tk2) {
0010     return tk1.stateAtBeamLine().trackStateAtPCA().position().z() <
0011            tk2.stateAtBeamLine().trackStateAtPCA().position().z();
0012   }
0013
0014 }  // namespace
0015
0016 GapClusterizerInZ::GapClusterizerInZ(const edm::ParameterSet& conf) {
0017   // some defaults to avoid uninitialized variables
0018   verbose_ = conf.getUntrackedParameter<bool>("verbose", false);
0019   zSep = conf.getParameter<double>("zSeparation");
0020   if (verbose_) {
0021     std::cout << "TrackClusterizerInZ:  algorithm=gap, zSeparation=" << zSep << std::endl;
0022   }
0023 }
0024
0025 float GapClusterizerInZ::zSeparation() const { return zSep; }
0026
0027 vector<vector<reco::TransientTrack> > GapClusterizerInZ::clusterize(const vector<reco::TransientTrack>& tracks) const {
0028   vector<reco::TransientTrack> tks = tracks;  // copy to be sorted
0029
0030   vector<vector<reco::TransientTrack> > clusters;
0031   if (tks.empty())
0032     return clusters;
0033
0034   // sort in increasing order of z
0035   stable_sort(tks.begin(), tks.end(), recTrackLessZ);
0036
0037   // init first cluster
0038   vector<reco::TransientTrack>::const_iterator it = tks.begin();
0039   vector<reco::TransientTrack> currentCluster;
0040   currentCluster.push_back(*it);
0041
0042   it++;
0043   for (; it != tks.end(); it++) {
0044     double zPrev = currentCluster.back().stateAtBeamLine().trackStateAtPCA().position().z();
0045     double zCurr = (*it).stateAtBeamLine().trackStateAtPCA().position().z();
0046
0047     if (abs(zCurr - zPrev) < zSeparation()) {
0048       // close enough ? cluster together
0049       currentCluster.push_back(*it);
0050     } else {
0051       // store current cluster, start new one
0052       clusters.push_back(currentCluster);
0053       currentCluster.clear();
0054       currentCluster.push_back(*it);
0055     }
0056   }
0057
0058   // store last cluster
0059   clusters.push_back(currentCluster);
0060
0061   return clusters;
0062 }
0063
0064 vector<TransientVertex> GapClusterizerInZ::vertices(const vector<reco::TransientTrack>& tracks) const {
0065   /* repackage track clusters, compatibility with newer clusterizers */
0066   std::vector<TransientVertex> primary_vertices;
0067   auto trackClusters = clusterize(tracks);
0068
0069   GlobalError dummyError(0.01, 0, 0.01, 0., 0., 0.01);
0070   for (auto& vertexTracks : trackClusters) {
0071     GlobalPoint position(0, 0, 0);  // dummy
0072     primary_vertices.push_back(TransientVertex(position, dummyError, vertexTracks, 0));
0073   }
0074
0075   return primary_vertices;
0076 }
0077
0078 void GapClusterizerInZ::fillPSetDescription(edm::ParameterSetDescription& desc) {
0079   desc.add<double>("zSeparation", 1.0);
0080   desc.addUntracked<bool>("verbose", false);
0081 }