MultiVertexFit/src/MultiVertexBSeeder.cc

0001 #include "RecoVertex/MultiVertexFit/interface/MultiVertexBSeeder.h"
0002 // #include "RecoVertex/AdaptiveVertexFit/interface/AdaptiveVertexFitter.h"
0003 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
0004 // #include "RecoVertex/VertexTools/interface/BeamSpot.h"
0005 #include "TrackingTools/PatternTools/interface/TwoTrackMinimumDistance.h"
0006 #include "RecoVertex/VertexTools/interface/FsmwModeFinder3d.h"
0007 #include "CommonTools/Clustering1D/interface/FsmwClusterizer1D.h"
0008 // #include "CommonTools/Clustering1D/interface/MultiClusterizer1D.h"
0009 #include "CommonTools/Clustering1D/interface/OutermostClusterizer1D.h"
0010 #include "DataFormats/GeometryCommonDetAlgo/interface/GlobalError.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012
0013 // #define MVBS_DEBUG
0014 #ifdef MVBS_DEBUG
0015 #include <map>
0016 #include "RecoVertex/MultiVertexFit/interface/DebuggingHarvester.h"
0017 #endif
0018
0019 using namespace std;
0020
0021 namespace {
0022   vector<reco::TransientTrack> convert(const vector<const reco::TransientTrack*>& ptrs) {
0023     vector<reco::TransientTrack> ret;
0024     for (vector<const reco::TransientTrack*>::const_iterator i = ptrs.begin(); i != ptrs.end(); ++i) {
0025       ret.push_back(**i);
0026     }
0027     return ret;
0028   }
0029
0030 #ifndef __clang__
0031   inline GlobalPoint& operator+=(GlobalPoint& a, const GlobalPoint& b) {
0032     a = GlobalPoint(a.x() + b.x(), a.y() + b.y(), a.z() + b.z());
0033     return a;
0034   }
0035 #endif
0036
0037   inline GlobalPoint& operator/=(GlobalPoint& a, float b) {
0038     a = GlobalPoint(a.x() / b, a.y() / b, a.z() / b);
0039     return a;
0040   }
0041
0042 #ifndef __clang__
0043   inline bool element(const reco::TransientTrack& rt, const TransientVertex& rv) {
0044     const vector<reco::TransientTrack> trks = rv.originalTracks();
0045     for (vector<reco::TransientTrack>::const_iterator i = trks.begin(); i != trks.end(); ++i) {
0046       if (rt == (*i))
0047         return true;
0048     };
0049     return false;
0050   }
0051 #endif
0052
0053   GlobalPoint toPoint(const GlobalVector& v) { return GlobalPoint(v.x(), v.y(), v.z()); }
0054
0055   GlobalVector toVector(const GlobalPoint& v) { return GlobalVector(v.x(), v.y(), v.z()); }
0056
0057   GlobalPoint computeJetOrigin(const vector<reco::TransientTrack>& trks) {
0058     FsmwModeFinder3d f;
0059     vector<ModeFinder3d::PointAndDistance> input;
0060     for (vector<reco::TransientTrack>::const_iterator i = trks.begin(); i != trks.end(); ++i) {
0061       input.push_back(ModeFinder3d::PointAndDistance(i->impactPointState().globalPosition(), 1.));
0062     }
0063     return f(input);
0064   }
0065
0066   GlobalVector computeJetDirection(const vector<reco::TransientTrack>& trks) {
0067     FsmwModeFinder3d f;
0068     vector<ModeFinder3d::PointAndDistance> input;
0069     for (vector<reco::TransientTrack>::const_iterator i = trks.begin(); i != trks.end(); ++i) {
0070       input.push_back(ModeFinder3d::PointAndDistance(toPoint(i->impactPointState().globalMomentum()), 1.));
0071     }
0072     GlobalPoint pt(f(input));
0073     pt /= pt.mag();
0074     return toVector(pt);
0075   }
0076
0077   GlobalTrajectoryParameters computeJetTrajectory(const vector<reco::TransientTrack>& trks) {
0078     /**
0079      *  construct a trajectory at the mean of
0080      *  the impact points of all tracks,
0081      *  momentum = total momentum of all tracks
0082      *
0083      */
0084     if (trks.empty())
0085       return GlobalTrajectoryParameters();
0086
0087     GlobalVector mom = computeJetDirection(trks);
0088     GlobalPoint pos = computeJetOrigin(trks);
0089
0090     GlobalTrajectoryParameters ret(pos, mom, 0, &(trks[0].impactPointState().globalParameters().magneticField()));
0091 #ifdef MVBS_DEBUG
0092     DebuggingHarvester("out.txt").save(ret, "p<sub>tot</sub>");
0093 #endif
0094     return ret;
0095   }
0096
0097   vector<Cluster1D<reco::TransientTrack> > computeIPs(const vector<reco::TransientTrack>& trks) {
0098     GlobalTrajectoryParameters jet = computeJetTrajectory(trks);
0099     FreeTrajectoryState axis(jet);
0100     TwoTrackMinimumDistance ttmd;
0101     vector<Cluster1D<reco::TransientTrack> > pts;
0102     for (vector<reco::TransientTrack>::const_iterator i = trks.begin(); i != trks.end(); ++i) {
0103       bool status = ttmd.calculate(axis, *(i->impactPointState().freeState()));
0104       if (status) {
0105         pair<GlobalPoint, GlobalPoint> pt = ttmd.points();
0106         double d = (pt.first - pt.second).mag();
0107         double w = 1. / (0.002 + d);  // hard coded weights
0108         double s = (pt.first - axis.position()).mag();
0109         Measurement1D ms(s, 1.0);
0110         vector<const reco::TransientTrack*> trk;
0111         trk.push_back(&(*i));
0112         pts.push_back(Cluster1D<reco::TransientTrack>(ms, trk, w));
0113       }
0114     }
0115     /*
0116     #ifdef MVBS_DEBUG
0117     map < string, harvest::MultiType > attrs;
0118     attrs["point:mag"]=0.5;
0119     attrs["point:color"]="khaki";
0120     DebuggingHarvester("out.txt").save ( pts, jet, attrs, "ip" );
0121     #endif
0122     */
0123     return pts;
0124   }
0125
0126 #ifndef __clang__
0127   inline GlobalPoint computePos(const GlobalTrajectoryParameters& jet, double s) {
0128     GlobalPoint ret = jet.position();
0129     ret += s * jet.momentum();
0130     return ret;
0131   }
0132 #endif
0133
0134   TransientVertex pseudoVertexFit(const Cluster1D<reco::TransientTrack>& src,
0135                                   bool ascending = false,
0136                                   bool kalmanfit = false) {
0137     // cout << "[MultiVertexBSeeder] debug: pseudoVertexFit with " << flush;
0138     vector<const reco::TransientTrack*> trkptrs = src.tracks();
0139     vector<reco::TransientTrack> trks = convert(trkptrs);
0140     // cout << trks.size() << " tracks.";
0141     GlobalPoint gp;
0142     GlobalError ge;
0143     if (kalmanfit) {
0144       TransientVertex v = KalmanVertexFitter().vertex(trks);
0145       gp = v.position();
0146     }
0147     TransientVertex ret = TransientVertex(gp, ge, trks, -1.);
0148     TransientVertex::TransientTrackToFloatMap mp;
0149     float w = 1.0;
0150     float r = 0.5;
0151     if (ascending) {
0152       w = pow((float)0.5, (int)(trks.size() - 1));
0153       r = 2.0;
0154     };
0155     for (vector<reco::TransientTrack>::const_iterator i = trks.begin(); i != trks.end(); ++i) {
0156       mp[*i] = w;
0157       w *= r;
0158     }
0159     ret.weightMap(mp);
0160     // cout << "[MultiVertexBSeeder] debug: return pseudoVertexFit with " << endl;
0161     return ret;
0162   }
0163
0164   /* TransientVertex kalmanVertexFit ( const Cluster1D < reco::TransientTrack > & src )
0165   {
0166     KalmanVertexFitter fitter;
0167       vector < const reco::TransientTrack * > trkptrs=src.tracks();
0168       vector < reco::TransientTrack > trks = convert ( trkptrs );
0169       return fitter.vertex ( trks );
0170     return TransientVertex();
0171   }*/
0172 }  // namespace
0173
0174 MultiVertexBSeeder* MultiVertexBSeeder::clone() const { return new MultiVertexBSeeder(*this); }
0175
0176 MultiVertexBSeeder::MultiVertexBSeeder(double nsigma) : theNSigma(nsigma) {}
0177
0178 vector<TransientVertex> MultiVertexBSeeder::vertices(const vector<reco::TransientTrack>& trks,
0179                                                      const reco::BeamSpot& s) const {
0180   return vertices(trks);
0181 }
0182
0183 vector<TransientVertex> MultiVertexBSeeder::vertices(const vector<reco::TransientTrack>& trks) const {
0184   vector<Cluster1D<reco::TransientTrack> > ips = computeIPs(trks);
0185   /*
0186   FsmwClusterizer1D < reco::TransientTrack > fsmw ( .4, theNSigma );
0187   MultiClusterizer1D<reco::TransientTrack> finder ( fsmw );
0188   */
0189   OutermostClusterizer1D<reco::TransientTrack> finder;
0190
0191   pair<vector<Cluster1D<reco::TransientTrack> >, vector<const reco::TransientTrack*> > res;
0192   res = finder(ips);
0193 #ifdef MVBS_DEBUG
0194   // need to compute jet trajectory again :-(
0195   GlobalTrajectoryParameters jet = computeJetTrajectory(trks);
0196   map<string, harvest::MultiType> attrs;
0197   attrs["point:mag"] = 0.75;
0198   attrs["point:color"] = "blue";
0199   attrs["point:name"] = "mode";
0200   DebuggingHarvester("out.txt").save(res.first, jet, attrs, "mode");
0201 #endif
0202
0203   vector<TransientVertex> ret;
0204   for (vector<Cluster1D<reco::TransientTrack> >::const_iterator i = res.first.begin(); i != res.first.end(); ++i) {
0205     ret.push_back(pseudoVertexFit(*i,
0206                                   i != res.first.begin(),
0207                                   /* kalman fit*/ true));
0208   }
0209   return ret;
0210 }
0211
0212 #ifdef MVBS_DEBUG
0213 #undef MVBS_DEBUG
0214 #endif