MultiVertexFit/src/MultiVertexReconstructor.cc

0001 #include "RecoVertex/MultiVertexFit/interface/MultiVertexReconstructor.h"
0002
0003 using namespace std;
0004
0005 namespace {
0006   typedef MultiVertexFitter::TrackAndWeight TrackAndWeight;
0007
0008   int verbose() { return 0; }
0009
0010 #ifndef __clang__
0011   inline void remove(vector<TransientVertex>& vtces, const vector<reco::TransientTrack>& trks) {
0012     cout << "[MultiVertexReconstructor] fixme remove not yet implemented" << endl;
0013     // remove trks from vtces
0014   }
0015 #endif
0016
0017   vector<vector<TrackAndWeight> > recover(const vector<TransientVertex>& vtces,
0018                                           const vector<reco::TransientTrack>& trks) {
0019     set<reco::TransientTrack> st;
0020     for (vector<reco::TransientTrack>::const_iterator i = trks.begin(); i != trks.end(); ++i) {
0021       st.insert(*i);
0022     }
0023
0024     vector<vector<TrackAndWeight> > bundles;
0025     for (vector<TransientVertex>::const_iterator vtx = vtces.begin(); vtx != vtces.end(); ++vtx) {
0026       vector<reco::TransientTrack> trks = vtx->originalTracks();
0027       vector<TrackAndWeight> tnws;
0028       for (vector<reco::TransientTrack>::const_iterator trk = trks.begin(); trk != trks.end(); ++trk) {
0029         float w = vtx->trackWeight(*trk);
0030         if (w > 1e-5) {
0031           TrackAndWeight tmp(*trk, w);
0032           set<reco::TransientTrack>::iterator pos = st.find(*trk);
0033           if (pos != st.end()) {
0034             st.erase(pos);
0035           }
0036           tnws.push_back(tmp);
0037         };
0038       };
0039       bundles.push_back(tnws);
0040     };
0041
0042     if (bundles.empty())
0043       return bundles;
0044
0045     // now add not-yet assigned tracks
0046     for (set<reco::TransientTrack>::const_iterator i = st.begin(); i != st.end(); ++i) {
0047       // cout << "[MultiVertexReconstructor] recovering " << i->id() << endl;
0048       TrackAndWeight tmp(*i, 0.);
0049       bundles[0].push_back(tmp);
0050     }
0051     return bundles;
0052   }
0053 }  // namespace
0054
0055 MultiVertexReconstructor::MultiVertexReconstructor(const VertexReconstructor& o,
0056                                                    const AnnealingSchedule& s,
0057                                                    float revive)
0058     : theOldReconstructor(o.clone()), theFitter(MultiVertexFitter(s, DefaultLinearizationPointFinder(), revive)) {}
0059
0060 MultiVertexReconstructor::~MultiVertexReconstructor() { delete theOldReconstructor; }
0061
0062 MultiVertexReconstructor* MultiVertexReconstructor::clone() const { return new MultiVertexReconstructor(*this); }
0063
0064 MultiVertexReconstructor::MultiVertexReconstructor(const MultiVertexReconstructor& o)
0065     : theOldReconstructor(o.theOldReconstructor->clone()), theFitter(o.theFitter) {}
0066
0067 vector<TransientVertex> MultiVertexReconstructor::vertices(const vector<reco::TransientTrack>& trks,
0068                                                            const reco::BeamSpot& s) const {
0069   return vertices(trks);
0070 }
0071
0072 vector<TransientVertex> MultiVertexReconstructor::vertices(const vector<reco::TransientTrack>& trks,
0073                                                            const vector<reco::TransientTrack>& primaries,
0074                                                            const reco::BeamSpot& s) const {
0075   return vertices(trks, primaries);
0076 }
0077
0078 vector<TransientVertex> MultiVertexReconstructor::vertices(const vector<reco::TransientTrack>& trks) const {
0079   /*
0080   cout << "[MultiVertexReconstructor] input trks: ";
0081   for ( vector< reco::TransientTrack >::const_iterator i=trks.begin();
0082         i!=trks.end() ; ++i )
0083   {
0084     cout << i->id() << "  ";
0085   }
0086   cout << endl;*/
0087   vector<TransientVertex> tmp = theOldReconstructor->vertices(trks);
0088   if (verbose()) {
0089     cout << "[MultiVertexReconstructor] non-freezing seeder found " << tmp.size() << " vertices from " << trks.size()
0090          << " tracks." << endl;
0091   }
0092   vector<vector<TrackAndWeight> > rc = recover(tmp, trks);
0093   vector<CachingVertex<5> > cvts = theFitter.vertices(rc);
0094   vector<TransientVertex> ret;
0095   for (vector<CachingVertex<5> >::const_iterator i = cvts.begin(); i != cvts.end(); ++i) {
0096     ret.push_back(*i);
0097   };
0098
0099   if (verbose()) {
0100     cout << "[MultiVertexReconstructor] input " << tmp.size() << " vertices, output " << ret.size() << " vertices."
0101          << endl;
0102   };
0103   return ret;
0104 }
0105
0106 vector<TransientVertex> MultiVertexReconstructor::vertices(const vector<reco::TransientTrack>& trks,
0107                                                            const vector<reco::TransientTrack>& primaries) const {
0108   /*
0109   cout << "[MultiVertexReconstructor] with " << primaries.size()
0110        << " primaries!" << endl;
0111        */
0112
0113   map<reco::TransientTrack, bool> st;
0114
0115   vector<reco::TransientTrack> total = trks;
0116   for (vector<reco::TransientTrack>::const_iterator i = trks.begin(); i != trks.end(); ++i) {
0117     st[(*i)] = true;
0118   }
0119
0120   // cout << "[MultiVertexReconstructor] FIXME dont just add up tracks. superpose" << endl;
0121   for (vector<reco::TransientTrack>::const_iterator i = primaries.begin(); i != primaries.end(); ++i) {
0122     if (!(st[(*i)])) {
0123       total.push_back(*i);
0124     }
0125   }
0126
0127   vector<TransientVertex> tmp = theOldReconstructor->vertices(total);
0128
0129   if (verbose()) {
0130     cout << "[MultiVertexReconstructor] freezing seeder found " << tmp.size() << " vertices from " << total.size()
0131          << " tracks." << endl;
0132   }
0133   vector<vector<TrackAndWeight> > rc = recover(tmp, trks);
0134   vector<CachingVertex<5> > cvts = theFitter.vertices(rc, primaries);
0135
0136   vector<TransientVertex> ret;
0137   for (vector<CachingVertex<5> >::const_iterator i = cvts.begin(); i != cvts.end(); ++i) {
0138     ret.push_back(*i);
0139   };
0140   if (verbose()) {
0141     cout << "[MultiVertexReconstructor] input " << tmp.size() << " vertices, output " << ret.size() << " vertices."
0142          << endl;
0143   };
0144   return ret;
0145 }
0146
0147 VertexReconstructor* MultiVertexReconstructor::reconstructor() const { return theOldReconstructor; }