PixelVertexFinding/src/DivisiveVertexFinder.cc

0001 #include "RecoVertex/PixelVertexFinding/interface/DivisiveVertexFinder.h"
0002 #include "RecoVertex/PixelVertexFinding/interface/PVPositionBuilder.h"
0003 #include "RecoVertex/PixelVertexFinding/interface/PVCluster.h"
0004 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "FWCore/Utilities/interface/isFinite.h"
0007 #include <utility>
0008 #include <vector>
0009 #include <map>
0010 #include <algorithm>
0011 #include <cmath>
0012
0013 DivisiveVertexFinder::DivisiveVertexFinder(double track_pt_min,
0014                                            double track_pt_max,
0015                                            double track_chi2_max,
0016                                            double track_prob_min,
0017                                            double zOffset,
0018                                            int ntrkMin,
0019                                            bool useError,
0020                                            double zSeparation,
0021                                            bool wtAverage,
0022                                            int verbosity)
0023     : zOffset_(zOffset),
0024       zSeparation_(zSeparation),
0025       ntrkMin_(ntrkMin),
0026       useError_(useError),
0027       wtAverage_(wtAverage),
0028       divmeth_(zOffset, ntrkMin, useError, zSeparation, wtAverage),
0029       verbose_(verbosity) {
0030   pvComparer_ = new PVClusterComparer(track_pt_min, track_pt_max, track_chi2_max, track_prob_min);
0031 }
0032
0033 DivisiveVertexFinder::~DivisiveVertexFinder() {}
0034
0035 bool DivisiveVertexFinder::findVertexes(const reco::TrackRefVector &trks,    // input
0036                                         reco::VertexCollection &vertexes) {  // output
0037   PVPositionBuilder pos;
0038   Measurement1D vz;
0039   if (wtAverage_) {
0040     vz = pos.wtAverage(trks);
0041   } else {
0042     vz = pos.average(trks);
0043   }
0044   reco::Vertex::Error err;
0045   err(2, 2) = vz.error() * vz.error();
0046
0047   reco::Vertex v(reco::Vertex::Point(0, 0, vz.value()), err, 0, 1, trks.size());
0048   for (unsigned int i = 0; i < trks.size(); i++) {
0049     double vz = trks[i]->vz();
0050     if (edm::isNotFinite(vz))
0051       continue;
0052     v.add(reco::TrackBaseRef(trks[i]));
0053   }
0054
0055   vertexes.push_back(v);
0056
0057   return true;
0058 }
0059
0060 bool DivisiveVertexFinder::findVertexesAlt(const reco::TrackRefVector &trks,  // input
0061                                            reco::VertexCollection &vertexes,
0062                                            const math::XYZPoint &bs) {  // output
0063   std::vector<PVCluster> in;
0064   std::pair<std::vector<PVCluster>, std::vector<const reco::Track *> > out;
0065
0066   // Convert input track collection into container needed by Wolfgang's templated code
0067   // Need to save a map to reconvert from bare pointers, oy vey
0068   std::map<const reco::Track *, reco::TrackRef> mapa;
0069   //  std::vector< std::vector< const reco::Track* > > trkps;
0070   for (unsigned int i = 0; i < trks.size(); ++i) {
0071     double vz = trks[i]->vz();
0072     if (edm::isNotFinite(vz))
0073       continue;
0074     std::vector<const reco::Track *> temp;
0075     temp.clear();
0076     temp.push_back(&(*trks[i]));
0077
0078     in.push_back(PVCluster(Measurement1D(trks[i]->dz(bs), trks[i]->dzError()), temp));
0079     mapa[temp[0]] = trks[i];
0080   }
0081
0082   if (verbose_ > 0) {
0083     edm::LogInfo("DivisiveVertexFinder") << "size of input vector of clusters " << in.size();
0084     for (unsigned int i = 0; i < in.size(); ++i) {
0085       edm::LogInfo("DivisiveVertexFinder")
0086           << "Track " << i << " addr " << in[i].tracks()[0] << " dz " << in[i].tracks()[0]->dz(bs) << " +- "
0087           << in[i].tracks()[0]->dzError() << " prodID " << mapa[in[i].tracks()[0]].id() << " dz from RefTrack "
0088           << mapa[in[i].tracks()[0]]->dz(bs) << " +- " << mapa[in[i].tracks()[0]]->dzError();
0089     }
0090   }
0091
0092   // Run the darn thing
0093   divmeth_.setBeamSpot(bs);
0094   out = divmeth_(in);
0095
0096   if (verbose_ > 0)
0097     edm::LogInfo("DivisiveVertexFinder") << " DivisiveClusterizer1D found " << out.first.size() << " vertexes";
0098
0099   // Now convert the output yet again into something we can safely store in the event
0100   for (unsigned int iv = 0; iv < out.first.size(); ++iv) {  // loop over output vertexes
0101     reco::Vertex::Error err;
0102     err(2, 2) = out.first[iv].position().error() * out.first[iv].position().error();
0103
0104     reco::Vertex v(reco::Vertex::Point(0, 0, out.first[iv].position().value()), err, 0, 1, out.second.size());
0105     if (verbose_ > 0)
0106       edm::LogInfo("DivisiveVertexFinder")
0107           << " DivisiveClusterizer1D vertex " << iv << " has " << out.first[iv].tracks().size()
0108           << " tracks and a position of " << v.z() << " +- " << std::sqrt(v.covariance(2, 2));
0109     for (unsigned int itrk = 0; itrk < out.first[iv].tracks().size(); ++itrk) {
0110       v.add(reco::TrackBaseRef(mapa[out.first[iv].tracks()[itrk]]));
0111     }
0112     vertexes.push_back(v);  // Done with horrible conversion, save it
0113   }
0114
0115   // Finally, sort the vertexes in decreasing sumPtSquared
0116   //  std::sort(vertexes.begin(), vertexes.end(), PVClusterComparer());
0117   std::sort(vertexes.begin(), vertexes.end(), *pvComparer_);
0118
0119   return true;
0120 }