TrimmedKalmanVertexFinder/src/ConfigurableTrimmedVertexFinder.cc

0001 #include "RecoVertex/TrimmedKalmanVertexFinder/interface/ConfigurableTrimmedVertexFinder.h"
0002 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0003
0004 using namespace reco;
0005
0006 ConfigurableTrimmedVertexFinder::ConfigurableTrimmedVertexFinder(const VertexFitter<5>* vf,
0007                                                                  const VertexUpdator<5>* vu,
0008                                                                  const VertexTrackCompatibilityEstimator<5>* ve)
0009     : theClusterFinder(vf, vu, ve),
0010       theVtxFitProbCut(0.01),
0011       theTrackCompatibilityToPV(0.05),
0012       theTrackCompatibilityToSV(0.01),
0013       theMaxNbOfVertices(0) {
0014   // default pt cut is 1.5 GeV
0015   theFilter.setPtCut(1.5);
0016 }
0017
0018 void ConfigurableTrimmedVertexFinder::setParameters(const edm::ParameterSet& s) {
0019   theFilter.setPtCut(s.getParameter<double>("ptCut"));
0020   theTrackCompatibilityToPV = s.getParameter<double>("trackCompatibilityToPVcut");
0021   theTrackCompatibilityToSV = s.getParameter<double>("trackCompatibilityToSVcut");
0022   theVtxFitProbCut = s.getParameter<double>("vtxFitProbCut");
0023   theMaxNbOfVertices = s.getParameter<int>("maxNbOfVertices");
0024 }
0025
0026 std::vector<TransientVertex> ConfigurableTrimmedVertexFinder::vertices(const std::vector<TransientTrack>& tracks) const {
0027   std::vector<TransientTrack> remaining;
0028
0029   return vertices(tracks, remaining, reco::BeamSpot(), false);
0030 }
0031
0032 std::vector<TransientVertex> ConfigurableTrimmedVertexFinder::vertices(const std::vector<TransientTrack>& tracks,
0033                                                                        const reco::BeamSpot& spot) const {
0034   std::vector<TransientTrack> remaining;
0035   return vertices(tracks, remaining, spot, true);
0036 }
0037
0038 std::vector<TransientVertex> ConfigurableTrimmedVertexFinder::vertices(const std::vector<TransientTrack>& tracks,
0039                                                                        std::vector<TransientTrack>& unused,
0040                                                                        const reco::BeamSpot& spot,
0041                                                                        bool use_spot) const {
0042   resetEvent(tracks);
0043   analyseInputTracks(tracks);
0044
0045   std::vector<TransientTrack> filtered;
0046   for (std::vector<TransientTrack>::const_iterator it = tracks.begin(); it != tracks.end(); it++) {
0047     if (theFilter(*it)) {
0048       filtered.push_back(*it);
0049     } else {
0050       unused.push_back(*it);
0051     }
0052   }
0053
0054   std::vector<TransientVertex> all = vertexCandidates(filtered, unused, spot, use_spot);
0055
0056   analyseVertexCandidates(all);
0057
0058   std::vector<TransientVertex> sel = clean(all);
0059
0060   analyseFoundVertices(sel);
0061
0062   return sel;
0063 }
0064
0065 std::vector<TransientVertex> ConfigurableTrimmedVertexFinder::vertexCandidates(const std::vector<TransientTrack>& tracks,
0066                                                                                std::vector<TransientTrack>& unused,
0067                                                                                const reco::BeamSpot& spot,
0068                                                                                bool use_spot) const {
0069   std::vector<TransientVertex> cand;
0070
0071   std::vector<TransientTrack> remain = tracks;
0072
0073   while (true) {
0074     float tkCompCut = (cand.empty() ? theTrackCompatibilityToPV : theTrackCompatibilityToSV);
0075
0076     //    std::cout << "PVR:compat cut " << tkCompCut << std::endl;
0077     theClusterFinder.setTrackCompatibilityCut(tkCompCut);
0078     //    std::cout << "PVCF:compat cut after setting "
0079     //   << theClusterFinder.trackCompatibilityCut() << std::endl;
0080
0081     std::vector<TransientVertex> newVertices;
0082     if (cand.empty() && use_spot) {
0083       newVertices = theClusterFinder.vertices(remain, spot);
0084     } else {
0085       newVertices = theClusterFinder.vertices(remain);
0086     }
0087     if (newVertices.empty())
0088       break;
0089
0090     analyseClusterFinder(newVertices, remain);
0091
0092     for (std::vector<TransientVertex>::const_iterator iv = newVertices.begin(); iv != newVertices.end(); iv++) {
0093       if (iv->originalTracks().size() > 1) {
0094         cand.push_back(*iv);
0095       } else {
0096         // candidate has too few tracks - get them back into the vector
0097         for (std::vector<TransientTrack>::const_iterator trk = iv->originalTracks().begin();
0098              trk != iv->originalTracks().end();
0099              ++trk) {
0100           unused.push_back(*trk);
0101         }
0102       }
0103     }
0104
0105     // when max number of vertices reached, stop
0106     if (theMaxNbOfVertices != 0) {
0107       if (cand.size() >= (unsigned int)theMaxNbOfVertices)
0108         break;
0109     }
0110   }
0111
0112   for (std::vector<TransientTrack>::const_iterator it = remain.begin(); it != remain.end(); it++) {
0113     unused.push_back(*it);
0114   }
0115
0116   return cand;
0117 }
0118
0119 std::vector<TransientVertex> ConfigurableTrimmedVertexFinder::clean(
0120     const std::vector<TransientVertex>& candidates) const {
0121   std::vector<TransientVertex> sel;
0122   for (std::vector<TransientVertex>::const_iterator i = candidates.begin(); i != candidates.end(); i++) {
0123     if (ChiSquaredProbability((*i).totalChiSquared(), (*i).degreesOfFreedom()) > theVtxFitProbCut) {
0124       sel.push_back(*i);
0125     }
0126   }
0127
0128   return sel;
0129 }