TrimmedKalmanVertexFinder/src/TrimmedVertexFinder.cc

0001 #include "RecoVertex/TrimmedKalmanVertexFinder/interface/TrimmedVertexFinder.h"
0002 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0003 #include "RecoVertex/VertexTools/interface/PerigeeLinearizedTrackState.h"
0004 #include "FWCore/Utilities/interface/Exception.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006
0007 using namespace reco;
0008
0009 TrimmedVertexFinder::TrimmedVertexFinder(const VertexFitter<5>* vf,
0010                                          const VertexUpdator<5>* vu,
0011                                          const VertexTrackCompatibilityEstimator<5>* ve)
0012     : theFitter(vf->clone()), theUpdator(vu->clone()), theEstimator(ve->clone()), theMinProb(0.05) {}
0013
0014 TrimmedVertexFinder::TrimmedVertexFinder(const TrimmedVertexFinder& other)
0015     : theFitter(other.theFitter->clone()),
0016       theUpdator(other.theUpdator->clone()),
0017       theEstimator(other.theEstimator->clone()),
0018       theMinProb(other.theMinProb) {}
0019
0020 TrimmedVertexFinder::~TrimmedVertexFinder() {
0021   delete theFitter;
0022   delete theUpdator;
0023   delete theEstimator;
0024 }
0025
0026 std::vector<TransientVertex> TrimmedVertexFinder::vertices(std::vector<TransientTrack>& tks) const {
0027   // FIXME write this!!!
0028   return vertices(tks, reco::BeamSpot(), false);
0029 }
0030
0031 std::vector<TransientVertex> TrimmedVertexFinder::vertices(std::vector<TransientTrack>& tks,
0032                                                            const reco::BeamSpot& spot,
0033                                                            bool use_spot) const {
0034   std::vector<TransientVertex> all;
0035   if (tks.size() < 2)
0036     return all;
0037
0038   // prepare vertex tracks and initial vertex
0039   CachingVertex<5> vtx;
0040   if (use_spot) {
0041     vtx = theFitter->vertex(tks, spot);
0042   } else {
0043     vtx = theFitter->vertex(tks);
0044   }
0045   if (!vtx.isValid()) {
0046     edm::LogWarning("TrimmedVertexFinder") << "initial vertex invalid."
0047                                            << " vertex finding stops here.";
0048     return all;
0049   }
0050
0051   std::vector<RefCountedVertexTrack> selected = vtx.tracks();
0052
0053   // reject incompatible tracks starting from the worst
0054   std::vector<RefCountedVertexTrack> remain;
0055   bool found = false;
0056   while (!found && selected.size() >= 2) {
0057     // find track with worst compatibility
0058     std::vector<RefCountedVertexTrack>::iterator iWorst = theWorst(vtx, selected, theMinProb);
0059
0060     if (iWorst != selected.end()) {
0061       // reject track
0062       remain.push_back(*iWorst);
0063       selected.erase(iWorst);
0064
0065       if (selected.size() == 1) {
0066         // not enough tracks to build new vertices
0067         remain.push_back(selected.front());
0068       } else {
0069         // removing track from vertex
0070         // need to redo the vertex fit instead of removing the track;
0071         // successive removals lead to numerical problems
0072         // this is however quick since intermediate quantities are cached
0073         // in the RefCountedVertexTracks
0074         if (use_spot)  // && all.size()==0 )
0075         {
0076           vtx = theFitter->vertex(selected, spot);
0077         } else {
0078           vtx = theFitter->vertex(selected);
0079         }
0080         if (!vtx.isValid()) {
0081           edm::LogWarning("TrimmedVertexFinder") << "current vertex invalid"
0082                                                  << "vertex finding stops here.";
0083           return all;
0084         }
0085
0086         // ref-counted tracks may have changed during the fit
0087         // if the linearization point has moved too much
0088         selected = vtx.tracks();
0089       }
0090
0091     } else {
0092       // no incompatible track remains, make vertex
0093       found = true;
0094       int n_tracks_in_vertex = selected.size();
0095
0096       // now return all tracks with weight < 0.5 to 'remain'.
0097       for (std::vector<RefCountedVertexTrack>::const_iterator t = selected.begin(); t != selected.end(); ++t) {
0098         if ((**t).weight() < 0.5) {
0099           /*
0100           cout << "[TrimmedVertexFinder] recycling track with weight "
0101                << (**t).weight() << endl;*/
0102           remain.push_back(*t);
0103           n_tracks_in_vertex--;  // one 'good' track less in the vertex
0104         };
0105       };
0106
0107       if (n_tracks_in_vertex > 1) {
0108         all.push_back(vtx);
0109       } else {
0110         edm::LogError("TrimmedVertexFinder") << "found vertex has less than 2 tracks";
0111       }
0112     }
0113   }
0114
0115   // modify list of incompatible tracks
0116   tks.clear();
0117   for (std::vector<RefCountedVertexTrack>::const_iterator i = remain.begin(); i != remain.end(); i++) {
0118     const PerigeeLinearizedTrackState* plts =
0119         dynamic_cast<const PerigeeLinearizedTrackState*>((**i).linearizedTrack().get());
0120     if (plts == nullptr) {
0121       throw cms::Exception("TrimmedVertexFinder: can't take track from non-perigee track state");
0122     }
0123
0124     tks.push_back(plts->track());
0125   }
0126
0127   // return 0 or 1 vertex
0128   return all;
0129 }
0130
0131 std::vector<TrimmedVertexFinder::RefCountedVertexTrack>::iterator TrimmedVertexFinder::theWorst(
0132     const CachingVertex<5>& vtx, std::vector<RefCountedVertexTrack>& vtxTracks, float cut) const {
0133   //  cout << "Cut is now " << cut << endl;
0134
0135   // find track with worst compatibility
0136   std::vector<RefCountedVertexTrack>::iterator iWorst = vtxTracks.end();
0137   float worseChi2 = 0.;
0138   for (std::vector<RefCountedVertexTrack>::iterator itr = vtxTracks.begin(); itr != vtxTracks.end(); itr++) {
0139     CachingVertex<5> newV = theUpdator->remove(vtx, *itr);
0140     if (!newV.isValid())
0141       return itr;
0142     std::pair<bool, double> result = theEstimator->estimate(newV, *itr);
0143     if (!result.first)
0144       return itr;
0145     float chi2 = result.second;
0146
0147     // compute number of degrees of freedom
0148     if (chi2 > 0.) {
0149       // we want to keep negative chi squares, and avoid calling
0150       // ChiSquaredProbability with a negative chi2.
0151       int ndf = 2;
0152       if (vtx.tracks().size() == 2)
0153         ndf = 1;
0154
0155       float prob = ChiSquaredProbability(chi2, ndf);
0156       if (prob < cut && chi2 >= worseChi2) {
0157         // small inconsistency: sorting on chi2, not on chi2 probability
0158         // because large chi2 are not rounded off while small probabilities
0159         // are rounded off to 0....
0160         // should not matter since at a given iteration
0161         // all track chi2 increments have the same ndf
0162         iWorst = itr;
0163         worseChi2 = chi2;
0164       }
0165     }
0166   }
0167
0168   return iWorst;
0169 }