CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoVertex/GhostTrackFitter/src/GhostTrackVertexFinder.cc

Go to the documentation of this file.
00001 #include <algorithm>
00002 #include <iterator>
00003 #include <vector>
00004 #include <map>
00005 #include <set>
00006 
00007 #include <Math/SVector.h>
00008 #include <Math/SMatrix.h>
00009 #include <Math/MatrixFunctions.h>
00010 
00011 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00012 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
00013 #include "DataFormats/GeometryCommonDetAlgo/interface/GlobalError.h"
00014 
00015 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00016 
00017 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
00018 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00019 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
00020 #include "TrackingTools/GeomPropagators/interface/AnalyticalImpactPointExtrapolator.h"
00021 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
00022 #include "TrackingTools/TransientTrack/interface/TransientTrackFromFTSFactory.h"
00023 
00024 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
00025 #include "RecoVertex/VertexPrimitives/interface/VertexFitter.h"
00026 #include "RecoVertex/VertexPrimitives/interface/VertexState.h"
00027 #include "RecoVertex/VertexPrimitives/interface/CachingVertex.h"
00028 #include "RecoVertex/VertexPrimitives/interface/LinearizedTrackState.h"
00029 #include "RecoVertex/VertexPrimitives/interface/ConvertError.h"
00030 #include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
00031 #include "RecoVertex/VertexTools/interface/LinearizedTrackStateFactory.h"
00032 #include "RecoVertex/VertexTools/interface/VertexTrackFactory.h"
00033 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
00034 #include "RecoVertex/VertexTools/interface/GeometricAnnealing.h"
00035 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
00036 #include "RecoVertex/AdaptiveVertexFit/interface/AdaptiveVertexFitter.h"
00037 
00038 #include "RecoVertex/GhostTrackFitter/interface/GhostTrack.h"
00039 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackState.h"
00040 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackFitter.h"
00041 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackPrediction.h"
00042 #include "RecoVertex/GhostTrackFitter/interface/SequentialGhostTrackFitter.h"
00043 #include "RecoVertex/GhostTrackFitter/interface/KalmanGhostTrackUpdater.h"
00044 
00045 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackVertexFinder.h"
00046 
00047 // #define DEBUG
00048 
00049 #ifdef DEBUG
00050 #       include "RecoBTag/SecondaryVertex/interface/TrackKinematics.h"
00051 #endif
00052 
00053 using namespace reco;
00054 
00055 namespace {
00056         using namespace ROOT::Math;
00057 
00058         typedef SVector<double, 3> Vector3;
00059 
00060         typedef SMatrix<double, 3, 3, MatRepSym<double, 3> > Matrix3S;
00061         typedef SMatrix<double, 2, 2, MatRepSym<double, 2> > Matrix2S;
00062         typedef SMatrix<double, 2, 3> Matrix23;
00063 
00064         typedef ReferenceCountingPointer<VertexTrack<5> >
00065                                         RefCountedVertexTrack;
00066         typedef ReferenceCountingPointer<LinearizedTrackState<5> >
00067                                         RefCountedLinearizedTrackState;
00068 
00069         struct VtxTrackIs : public std::unary_function<
00070                                         RefCountedVertexTrack, bool> {
00071                 VtxTrackIs(const TransientTrack &track) : track(track) {}
00072                 bool operator()(const RefCountedVertexTrack &vtxTrack) const
00073                 { return vtxTrack->linearizedTrack()->track() == track; }
00074 
00075                 const TransientTrack &track;
00076         };
00077 
00078         static inline Vector3 conv(const GlobalVector &vec)
00079         {
00080                 Vector3 result;
00081                 result[0] = vec.x();
00082                 result[1] = vec.y();
00083                 result[2] = vec.z();
00084                 return result;
00085         }
00086 
00087         static inline double sqr(double arg) { return arg * arg; }
00088 }
00089 
00090 struct GhostTrackVertexFinder::FinderInfo {
00091         FinderInfo(const CachingVertex<5> &primaryVertex,
00092                    const GhostTrack &ghostTrack, const BeamSpot &beamSpot,
00093                    bool hasBeamSpot, bool hasPrimaries) :
00094                 primaryVertex(primaryVertex), pred(ghostTrack.prediction()),
00095                 prior(ghostTrack.prior()), states(ghostTrack.states()),
00096                 beamSpot(beamSpot), hasBeamSpot(hasBeamSpot),
00097                 hasPrimaries(hasPrimaries), field(0) {}
00098 
00099         const CachingVertex<5>                  &primaryVertex;
00100         GhostTrackPrediction                    pred;
00101         GhostTrackPrediction                    prior;
00102         std::vector<GhostTrackState>            states;
00103         VertexState                             beamSpot;
00104         bool                                    hasBeamSpot;
00105         bool                                    hasPrimaries;
00106         const MagneticField                     *field;
00107         TransientTrack                          ghostTrack;
00108 };
00109 
00110 static TransientTrack transientGhostTrack(const GhostTrackPrediction &pred,
00111                                           const MagneticField *field)
00112 { return TransientTrackFromFTSFactory().build(pred.fts(field)); }
00113 
00114 static double vtxErrorLong(const GlobalError &error, const GlobalVector &dir)
00115 {
00116         return ROOT::Math::Similarity(conv(dir), error.matrix_new());
00117 }
00118 
00119 static GlobalPoint vtxMean(const GlobalPoint &p1, const GlobalError &e1,
00120                            const GlobalPoint &p2, const GlobalError &e2)
00121 {
00122         GlobalVector diff = p2 - p1;
00123 
00124         double err1 = vtxErrorLong(e1, diff);
00125         double err2 = vtxErrorLong(e2, diff);
00126 
00127         double weight = err1 / (err1 + err2);
00128 
00129         return p1 + weight * diff;
00130 }
00131 
00132 static VertexState stateMean(const VertexState &v1, const VertexState &v2)
00133 {
00134         return VertexState(vtxMean(v1.position(), v1.error(),
00135                                    v2.position(), v2.error()),
00136                            v1.error() + v2.error());
00137 }
00138 
00139 static bool covarianceUpdate(Matrix3S &cov, const Vector3 &residual,
00140                              const Matrix3S &error, double &chi2,
00141                              double theta, double phi)
00142 {
00143         using namespace ROOT::Math;
00144 
00145         Matrix23 jacobian;
00146         jacobian(0, 0) = std::cos(phi) * std::cos(theta);
00147         jacobian(0, 1) = std::sin(phi) * std::cos(theta);
00148         jacobian(0, 2) = -std::sin(theta);
00149         jacobian(1, 0) = -std::sin(phi);
00150         jacobian(1, 1) = std::cos(phi);
00151 
00152         Matrix2S measErr = Similarity(jacobian, error);
00153         Matrix2S combErr = Similarity(jacobian, cov) + measErr;
00154         if (!measErr.Invert() || !combErr.Invert())
00155                 return false;
00156 
00157         cov -= SimilarityT(jacobian * cov, combErr);
00158         chi2 += Similarity(jacobian * residual, measErr);
00159 
00160         return true;
00161 }
00162 
00163 static CachingVertex<5> vertexAtState(const TransientTrack &ghostTrack,
00164                                       const GhostTrackPrediction &pred,
00165                                       const GhostTrackState &state)
00166 {
00167         LinearizedTrackStateFactory linTrackFactory;
00168         VertexTrackFactory<5> vertexTrackFactory;
00169 
00170         if (!state.isValid())
00171                 return CachingVertex<5>();
00172 
00173         GlobalPoint pca1 = pred.position(state.lambda());
00174         GlobalError err1 = pred.positionError(state.lambda());
00175 
00176         GlobalPoint pca2 = state.globalPosition();
00177         GlobalError err2 = state.cartesianError();
00178 
00179         GlobalPoint point = vtxMean(pca1, err1, pca2, err2);
00180 
00181         TransientTrack recTrack = state.track();
00182 
00183         RefCountedLinearizedTrackState linState[2] = {
00184                 linTrackFactory.linearizedTrackState(point, ghostTrack),
00185                 linTrackFactory.linearizedTrackState(point, recTrack)
00186         };
00187         if ( !linState[0]->isValid() || !linState[1]->isValid() )
00188           return CachingVertex<5>();
00189 
00190         Matrix3S cov = SMatrixIdentity();
00191         cov *= 10000;
00192         double chi2 = 0.;
00193         if (!covarianceUpdate(cov, conv(pca1 - point), err1.matrix_new(), chi2,
00194                                linState[0]->predictedStateParameters()[1],
00195                                linState[0]->predictedStateParameters()[2]) ||
00196             !covarianceUpdate(cov, conv(pca2 - point), err2.matrix_new(), chi2,
00197                                linState[1]->predictedStateParameters()[1],
00198                                linState[1]->predictedStateParameters()[2]))
00199                 return CachingVertex<5>();
00200 
00201         GlobalError error(cov);
00202         VertexState vtxState(point, error);
00203 
00204         std::vector<RefCountedVertexTrack> linTracks(2);
00205         linTracks[0] = vertexTrackFactory.vertexTrack(linState[0], vtxState);
00206         linTracks[1] = vertexTrackFactory.vertexTrack(linState[1], vtxState);
00207 
00208         return CachingVertex<5>(point, error, linTracks, chi2);
00209 }
00210 
00211 static RefCountedVertexTrack relinearizeTrack(
00212                         const RefCountedVertexTrack &track,
00213                         const VertexState &state,
00214                         const VertexTrackFactory<5> factory)
00215 {
00216         RefCountedLinearizedTrackState linTrack = track->linearizedTrack();
00217         linTrack = linTrack->stateWithNewLinearizationPoint(state.position());
00218         return factory.vertexTrack(linTrack, state);
00219 }
00220 
00221 static std::vector<RefCountedVertexTrack> relinearizeTracks(
00222                         const std::vector<RefCountedVertexTrack> &tracks,
00223                         const VertexState &state)
00224 {
00225         VertexTrackFactory<5> vertexTrackFactory;
00226 
00227         std::vector<RefCountedVertexTrack> finalTracks;
00228         finalTracks.reserve(tracks.size());
00229 
00230         for(std::vector<RefCountedVertexTrack>::const_iterator iter =
00231                                 tracks.begin(); iter != tracks.end(); ++iter)
00232                 finalTracks.push_back(
00233                         relinearizeTrack(*iter, state, vertexTrackFactory));
00234 
00235         return finalTracks;
00236 }
00237 
00238 double GhostTrackVertexFinder::vertexCompat(const CachingVertex<5> &vtx1,
00239                            const CachingVertex<5> &vtx2,
00240                            const FinderInfo &info,
00241                            double scale1, double scale2)
00242 {
00243         Vector3 diff = conv(vtx2.position() - vtx1.position());
00244         Matrix3S cov = scale1 * vtx1.error().matrix_new() +
00245                        scale2 * vtx2.error().matrix_new();
00246 
00247         return sqr(ROOT::Math::Dot(diff, diff)) / ROOT::Math::Similarity(cov, diff);
00248 }
00249 
00250 static double trackVertexCompat(const CachingVertex<5> &vtx,
00251                                 const RefCountedVertexTrack &vertexTrack)
00252 {
00253         using namespace ROOT::Math;
00254 
00255         TransientTrack track = vertexTrack->linearizedTrack()->track();
00256         GlobalPoint point = vtx.position();
00257         AnalyticalImpactPointExtrapolator extrap(track.field());
00258         TrajectoryStateOnSurface tsos =
00259                         extrap.extrapolate(track.impactPointState(), point);
00260 
00261         if (!tsos.isValid())
00262                 return 1.0e6;
00263 
00264         GlobalPoint point1 = vtx.position();
00265         GlobalPoint point2 = tsos.globalPosition();
00266         Vector3 dir = conv(point2 - point1);
00267         Matrix3S error = vtx.error().matrix_new() +
00268                          tsos.cartesianError().matrix().Sub<Matrix3S>(0, 0);
00269         if (!error.Invert())
00270                 return 1.0e6;
00271 
00272         return ROOT::Math::Similarity(error, dir);
00273 }
00274 
00275 GhostTrackVertexFinder::GhostTrackVertexFinder() :
00276         maxFitChi2_(10.0),
00277         mergeThreshold_(3.0),
00278         primcut_(2.0),
00279         seccut_(5.0),
00280         fitType_(kAlwaysWithGhostTrack)
00281 {
00282 }
00283 
00284 GhostTrackVertexFinder::GhostTrackVertexFinder(
00285                 double maxFitChi2, double mergeThreshold, double primcut,
00286                 double seccut, FitType fitType) :
00287         maxFitChi2_(maxFitChi2),
00288         mergeThreshold_(mergeThreshold),
00289         primcut_(primcut),
00290         seccut_(seccut),
00291         fitType_(fitType)
00292 {
00293 }
00294 
00295 GhostTrackVertexFinder::~GhostTrackVertexFinder()
00296 {
00297 }
00298 
00299 GhostTrackFitter &GhostTrackVertexFinder::ghostTrackFitter() const
00300 {
00301         if (!ghostTrackFitter_.get())
00302                 ghostTrackFitter_.reset(new GhostTrackFitter);
00303 
00304         return *ghostTrackFitter_;
00305 }
00306 
00307 VertexFitter<5> &GhostTrackVertexFinder::vertexFitter(bool primary) const
00308 {
00309         std::auto_ptr<VertexFitter<5> > *ptr =
00310                         primary ? &primVertexFitter_ : &secVertexFitter_;
00311         if (!ptr->get())
00312                 ptr->reset(new AdaptiveVertexFitter(
00313                                         GeometricAnnealing(primary
00314                                                 ? primcut_ : seccut_)));
00315 
00316         return **ptr;
00317 }
00318 
00319 #ifdef DEBUG
00320 static void debugTrack(const TransientTrack &track, const TransientVertex *vtx = 0)
00321 {
00322         const FreeTrajectoryState &fts = track.initialFreeState();
00323         GlobalVector mom = fts.momentum();
00324         std::cout << "\tpt = " << mom.perp() << ", eta = " << mom.eta() << ", phi = " << mom.phi() << ", charge = " << fts.charge();
00325         if (vtx && vtx->hasTrackWeight())
00326                 std::cout << ", w = " << vtx->trackWeight(track) << std::endl;
00327         else
00328                 std::cout << std::endl;
00329 }
00330 
00331 static void debugTracks(const std::vector<TransientTrack> &tracks, const TransientVertex *vtx = 0)
00332 {
00333         TrackKinematics kin;
00334         for(std::vector<TransientTrack>::const_iterator iter = tracks.begin();
00335             iter != tracks.end(); ++iter) {
00336                 debugTrack(*iter, vtx);
00337                 try {
00338                         TrackBaseRef track = iter->trackBaseRef();
00339                         kin.add(*track);
00340                 } catch ( exception const & e ) {
00341                         // ignore, yeah this is debugging, so shut up ^^
00342                 }
00343         }
00344 
00345         std::cout << "mass = " << kin.vectorSum().M() << std::endl;
00346 }
00347 
00348 static void debugTracks(const std::vector<RefCountedVertexTrack> &tracks)
00349 {
00350         std::vector<TransientTrack> tracks_;
00351         for(std::vector<RefCountedVertexTrack>::const_iterator iter = tracks.begin();
00352             iter != tracks.end(); ++iter) {
00353                 std::cout << "+ " << (*iter)->linearizedTrack()->linearizationPoint() << std::endl;
00354                 tracks_.push_back((*iter)->linearizedTrack()->track());
00355         }
00356         debugTracks(tracks_);
00357 }
00358 
00359 static void debugVertex(const TransientVertex &vtx, const GhostTrackPrediction &pred)
00360 {
00361         double origin = 0.;
00362         std::cout << vtx.position() << "-> " << vtx.totalChiSquared() << " with " << vtx.degreesOfFreedom() << std::endl;
00363 
00364         double err = std::sqrt(vtxErrorLong(
00365                         vtx.positionError(),
00366                         pred.direction())
00367                                         / pred.rho2());
00368         std::cout << "* " << (pred.lambda(vtx.position()) * pred.rho() - origin)
00369                   << " +-" << err
00370                   << " => " << ((pred.lambda(vtx.position()) * pred.rho() - origin) / err)
00371                   << std::endl;
00372 
00373         std::vector<TransientTrack> tracks = vtx.originalTracks();
00374         debugTracks(tracks, &vtx);
00375 }
00376 #endif
00377 
00378 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(
00379                         const Vertex &primaryVertex,
00380                         const GlobalVector &direction,
00381                         double coneRadius,
00382                         const std::vector<TransientTrack> &tracks) const
00383 {
00384         return vertices(RecoVertex::convertPos(primaryVertex.position()),
00385                         RecoVertex::convertError(primaryVertex.error()),
00386                         direction, coneRadius, tracks);
00387 }
00388 
00389 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(
00390                         const Vertex &primaryVertex,
00391                         const GlobalVector &direction,
00392                         double coneRadius,
00393                         const reco::BeamSpot &beamSpot,
00394                         const std::vector<TransientTrack> &tracks) const
00395 {
00396         return vertices(RecoVertex::convertPos(primaryVertex.position()),
00397                         RecoVertex::convertError(primaryVertex.error()),
00398                         direction, coneRadius, beamSpot, tracks);
00399 }
00400 
00401 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(
00402                         const Vertex &primaryVertex,
00403                         const GlobalVector &direction,
00404                         double coneRadius,
00405                         const reco::BeamSpot &beamSpot,
00406                         const std::vector<TransientTrack> &primaries,
00407                         const std::vector<TransientTrack> &tracks) const
00408 {
00409         return vertices(RecoVertex::convertPos(primaryVertex.position()),
00410                         RecoVertex::convertError(primaryVertex.error()),
00411                         direction, coneRadius, beamSpot, primaries, tracks);
00412 }
00413 
00414 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(
00415                         const GlobalPoint &primaryPosition,
00416                         const GlobalError &primaryError,
00417                         const GlobalVector &direction,
00418                         double coneRadius,
00419                         const std::vector<TransientTrack> &tracks) const
00420 {
00421         GhostTrack ghostTrack =
00422                 ghostTrackFitter().fit(primaryPosition, primaryError,
00423                                        direction, coneRadius, tracks);
00424 
00425         CachingVertex<5> primary(primaryPosition, primaryError,
00426                                  std::vector<RefCountedVertexTrack>(), 0.);
00427 
00428         std::vector<TransientVertex> result = vertices(ghostTrack, primary);
00429 
00430 #ifdef DEBUG
00431         for(std::vector<TransientVertex>::const_iterator iter = result.begin();
00432             iter != result.end(); ++iter)
00433                 debugVertex(*iter, ghostTrack.prediction());
00434 #endif
00435 
00436         return result;
00437 }
00438 
00439 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(
00440                         const GlobalPoint &primaryPosition,
00441                         const GlobalError &primaryError,
00442                         const GlobalVector &direction,
00443                         double coneRadius,
00444                         const BeamSpot &beamSpot,
00445                         const std::vector<TransientTrack> &tracks) const
00446 {
00447         GhostTrack ghostTrack =
00448                 ghostTrackFitter().fit(primaryPosition, primaryError,
00449                                        direction, coneRadius, tracks);
00450 
00451         CachingVertex<5> primary(primaryPosition, primaryError,
00452                                  std::vector<RefCountedVertexTrack>(), 0.);
00453 
00454         std::vector<TransientVertex> result =
00455                                 vertices(ghostTrack, primary, beamSpot, true);
00456 
00457 #ifdef DEBUG
00458         for(std::vector<TransientVertex>::const_iterator iter = result.begin();
00459             iter != result.end(); ++iter)
00460                 debugVertex(*iter, ghostTrack.prediction());
00461 #endif
00462 
00463         return result;
00464 }
00465 
00466 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(
00467                         const GlobalPoint &primaryPosition,
00468                         const GlobalError &primaryError,
00469                         const GlobalVector &direction,
00470                         double coneRadius,
00471                         const BeamSpot &beamSpot,
00472                         const std::vector<TransientTrack> &primaries,
00473                         const std::vector<TransientTrack> &tracks) const
00474 {
00475         GhostTrack ghostTrack =
00476                 ghostTrackFitter().fit(primaryPosition, primaryError,
00477                                        direction, coneRadius, tracks);
00478 
00479         std::vector<RefCountedVertexTrack> primaryVertexTracks;
00480         if (!primaries.empty()) {
00481                 LinearizedTrackStateFactory linTrackFactory;
00482                 VertexTrackFactory<5> vertexTrackFactory;
00483 
00484                 VertexState state(primaryPosition, primaryError);
00485 
00486                 for(std::vector<TransientTrack>::const_iterator iter =
00487                         primaries.begin(); iter != primaries.end(); ++iter) {
00488 
00489                         RefCountedLinearizedTrackState linState =
00490                                 linTrackFactory.linearizedTrackState(
00491                                                 primaryPosition, *iter);
00492 
00493                         primaryVertexTracks.push_back(
00494                                         vertexTrackFactory.vertexTrack(
00495                                                         linState, state));
00496                 }
00497         }
00498 
00499         CachingVertex<5> primary(primaryPosition, primaryError,
00500                                  primaryVertexTracks, 0.);
00501 
00502         std::vector<TransientVertex> result =
00503                         vertices(ghostTrack, primary, beamSpot, true, true);
00504 
00505 #ifdef DEBUG
00506         for(std::vector<TransientVertex>::const_iterator iter = result.begin();
00507             iter != result.end(); ++iter)
00508                 debugVertex(*iter, ghostTrack.prediction());
00509 #endif
00510 
00511         return result;
00512 }
00513 
00514 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(
00515                         const GlobalPoint &primaryPosition,
00516                         const GlobalError &primaryError,
00517                         const GhostTrack &ghostTrack) const
00518 {
00519         CachingVertex<5> primary(primaryPosition, primaryError,
00520                                  std::vector<RefCountedVertexTrack>(), 0.);
00521 
00522         std::vector<TransientVertex> result = vertices(ghostTrack, primary);
00523 
00524 #ifdef DEBUG
00525         for(std::vector<TransientVertex>::const_iterator iter = result.begin();
00526             iter != result.end(); ++iter)
00527                 debugVertex(*iter, ghostTrack.prediction());
00528 #endif
00529 
00530         return result;
00531 }
00532 
00533 
00534 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(
00535                         const GlobalPoint &primaryPosition,
00536                         const GlobalError &primaryError,
00537                         const BeamSpot &beamSpot,
00538                         const GhostTrack &ghostTrack) const
00539 {
00540         CachingVertex<5> primary(primaryPosition, primaryError,
00541                                  std::vector<RefCountedVertexTrack>(), 0.);
00542 
00543         std::vector<TransientVertex> result =
00544                                 vertices(ghostTrack, primary, beamSpot, true);
00545 
00546 #ifdef DEBUG
00547         for(std::vector<TransientVertex>::const_iterator iter = result.begin();
00548             iter != result.end(); ++iter)
00549                 debugVertex(*iter, ghostTrack.prediction());
00550 #endif
00551 
00552         return result;
00553 }
00554 
00555 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(
00556                         const GlobalPoint &primaryPosition,
00557                         const GlobalError &primaryError,
00558                         const BeamSpot &beamSpot,
00559                         const std::vector<TransientTrack> &primaries,
00560                         const GhostTrack &ghostTrack) const
00561 {
00562         std::vector<RefCountedVertexTrack> primaryVertexTracks;
00563         if (!primaries.empty()) {
00564                 LinearizedTrackStateFactory linTrackFactory;
00565                 VertexTrackFactory<5> vertexTrackFactory;
00566 
00567                 VertexState state(primaryPosition, primaryError);
00568 
00569                 for(std::vector<TransientTrack>::const_iterator iter =
00570                         primaries.begin(); iter != primaries.end(); ++iter) {
00571 
00572                         RefCountedLinearizedTrackState linState =
00573                                 linTrackFactory.linearizedTrackState(
00574                                                 primaryPosition, *iter);
00575 
00576                         primaryVertexTracks.push_back(
00577                                         vertexTrackFactory.vertexTrack(
00578                                                         linState, state));
00579                 }
00580         }
00581 
00582         CachingVertex<5> primary(primaryPosition, primaryError,
00583                                  primaryVertexTracks, 0.);
00584 
00585         std::vector<TransientVertex> result =
00586                         vertices(ghostTrack, primary, beamSpot, true, true);
00587 
00588 #ifdef DEBUG
00589         for(std::vector<TransientVertex>::const_iterator iter = result.begin();
00590             iter != result.end(); ++iter)
00591                 debugVertex(*iter, ghostTrack.prediction());
00592 #endif
00593 
00594         return result;
00595 }
00596 
00597 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(
00598                         const Vertex &primaryVertex,
00599                         const GhostTrack &ghostTrack) const
00600 {
00601         return vertices(RecoVertex::convertPos(primaryVertex.position()),
00602                         RecoVertex::convertError(primaryVertex.error()),
00603                         ghostTrack);
00604 }
00605 
00606 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(
00607                         const Vertex &primaryVertex,
00608                         const BeamSpot &beamSpot,
00609                         const GhostTrack &ghostTrack) const
00610 {
00611         return vertices(RecoVertex::convertPos(primaryVertex.position()),
00612                         RecoVertex::convertError(primaryVertex.error()),
00613                         beamSpot, ghostTrack);
00614 }
00615 
00616 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(
00617                         const Vertex &primaryVertex,
00618                         const BeamSpot &beamSpot,
00619                         const std::vector<TransientTrack> &primaries,
00620                         const GhostTrack &ghostTrack) const
00621 {
00622         return vertices(RecoVertex::convertPos(primaryVertex.position()),
00623                         RecoVertex::convertError(primaryVertex.error()),
00624                         beamSpot, primaries, ghostTrack);
00625 }
00626 
00627 static GhostTrackPrediction dummyPrediction(
00628                 const Vertex &primaryVertex, const Track &ghostTrack)
00629 {
00630         return GhostTrackPrediction(
00631                         RecoVertex::convertPos(primaryVertex.position()),
00632                         RecoVertex::convertError(primaryVertex.error()),
00633                         GlobalVector(ghostTrack.px(), ghostTrack.py(),
00634                                      ghostTrack.pz()), 0.05);
00635 }
00636 
00637 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(
00638                         const Vertex &primaryVertex,
00639                         const Track &ghostTrack,
00640                         const std::vector<TransientTrack> &tracks,
00641                         const std::vector<float> &weights) const
00642 {
00643         GhostTrack ghostTrack_(ghostTrack, tracks, weights,
00644                                dummyPrediction(primaryVertex, ghostTrack),
00645                                RecoVertex::convertPos(
00646                                                 primaryVertex.position()),
00647                                true);
00648 
00649         CachingVertex<5> primary(
00650                         RecoVertex::convertPos(primaryVertex.position()),
00651                         RecoVertex::convertError(primaryVertex.error()),
00652                         std::vector<RefCountedVertexTrack>(), 0.);
00653 
00654         std::vector<TransientVertex> result = vertices(ghostTrack_, primary);
00655 
00656 #ifdef DEBUG
00657         for(std::vector<TransientVertex>::const_iterator iter = result.begin();
00658             iter != result.end(); ++iter)
00659                 debugVertex(*iter, ghostTrack_.prediction());
00660 #endif
00661 
00662         return result;
00663 }
00664 
00665 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(
00666                         const Vertex &primaryVertex,
00667                         const Track &ghostTrack,
00668                         const BeamSpot &beamSpot,
00669                         const std::vector<TransientTrack> &tracks,
00670                         const std::vector<float> &weights) const
00671 {
00672         GhostTrack ghostTrack_(ghostTrack, tracks, weights,
00673                                dummyPrediction(primaryVertex, ghostTrack),
00674                                RecoVertex::convertPos(primaryVertex.position()),
00675                                true);
00676 
00677         CachingVertex<5> primary(
00678                         RecoVertex::convertPos(primaryVertex.position()),
00679                         RecoVertex::convertError(primaryVertex.error()),
00680                         std::vector<RefCountedVertexTrack>(), 0.);
00681 
00682         std::vector<TransientVertex> result =
00683                         vertices(ghostTrack_, primary, beamSpot, true);
00684 
00685 #ifdef DEBUG
00686         for(std::vector<TransientVertex>::const_iterator iter = result.begin();
00687             iter != result.end(); ++iter)
00688                 debugVertex(*iter, ghostTrack_.prediction());
00689 #endif
00690 
00691         return result;
00692 }
00693 
00694 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(
00695                         const Vertex &primaryVertex,
00696                         const Track &ghostTrack,
00697                         const BeamSpot &beamSpot,
00698                         const std::vector<TransientTrack> &primaries,
00699                         const std::vector<TransientTrack> &tracks,
00700                         const std::vector<float> &weights) const
00701 {
00702         GhostTrack ghostTrack_(ghostTrack, tracks, weights,
00703                                dummyPrediction(primaryVertex, ghostTrack),
00704                                RecoVertex::convertPos(primaryVertex.position()),
00705                                true);
00706 
00707         std::vector<RefCountedVertexTrack> primaryVertexTracks;
00708         if (!primaries.empty()) {
00709                 LinearizedTrackStateFactory linTrackFactory;
00710                 VertexTrackFactory<5> vertexTrackFactory;
00711 
00712                 VertexState state(
00713                         RecoVertex::convertPos(primaryVertex.position()),
00714                         RecoVertex::convertError(primaryVertex.error()));
00715 
00716                 for(std::vector<TransientTrack>::const_iterator iter =
00717                         primaries.begin(); iter != primaries.end(); ++iter) {
00718 
00719                         RefCountedLinearizedTrackState linState =
00720                                 linTrackFactory.linearizedTrackState(
00721                                                 state.position(), *iter);
00722 
00723                         primaryVertexTracks.push_back(
00724                                         vertexTrackFactory.vertexTrack(
00725                                                         linState, state));
00726                 }
00727         }
00728 
00729         CachingVertex<5> primary(
00730                         RecoVertex::convertPos(primaryVertex.position()),
00731                         RecoVertex::convertError(primaryVertex.error()),
00732                         primaryVertexTracks, 0.);
00733 
00734         std::vector<TransientVertex> result =
00735                         vertices(ghostTrack_, primary, beamSpot, true, true);
00736 
00737 #ifdef DEBUG
00738         for(std::vector<TransientVertex>::const_iterator iter = result.begin();
00739             iter != result.end(); ++iter)
00740                 debugVertex(*iter, ghostTrack_.prediction());
00741 #endif
00742 
00743         return result;
00744 }
00745 
00746 static double fitChi2(const CachingVertex<5> &vtx)
00747 {
00748         return vtx.totalChiSquared() / vtx.degreesOfFreedom();
00749 }
00750 
00751 std::vector<CachingVertex<5> > GhostTrackVertexFinder::initialVertices(
00752                                                 const FinderInfo &info) const
00753 {
00754         std::vector<CachingVertex<5> > vertices;
00755         for(std::vector<GhostTrackState>::const_iterator iter =
00756                     info.states.begin(); iter != info. states.end(); ++iter) {
00757 
00758                 if (!iter->isValid())
00759                         continue;
00760 
00761                 GhostTrackState state(*iter);
00762                 state.linearize(info.pred);
00763 
00764                 if (!state.isValid())
00765                         continue;
00766 
00767                 CachingVertex<5> vtx = vertexAtState(info.ghostTrack,
00768                                                      info.pred, state);
00769 
00770                 if (vtx.isValid()) // && fitChi2(vtx) < maxFitChi2_)
00771                         vertices.push_back(vtx);
00772         }
00773 
00774         return vertices;
00775 }
00776 
00777 static void mergeTrackHelper(const std::vector<RefCountedVertexTrack> &tracks,
00778                              std::vector<RefCountedVertexTrack> &newTracks,
00779                              const VertexState &state,
00780                              const VtxTrackIs &ghostTrackFinder,
00781                              RefCountedVertexTrack &ghostTrack,
00782                              const VertexTrackFactory<5> &factory)
00783 {
00784         for(std::vector<RefCountedVertexTrack>::const_iterator iter =
00785                         tracks.begin(); iter != tracks.end(); ++iter) {
00786                 bool gt = ghostTrackFinder(*iter);
00787                 if (gt && ghostTrack)
00788                         continue;
00789 
00790                 RefCountedVertexTrack track =
00791                                 relinearizeTrack(*iter, state, factory);
00792 
00793                 if (gt)
00794                         ghostTrack = *iter;
00795                 else
00796                         newTracks.push_back(*iter);
00797         }
00798 }
00799 
00800 CachingVertex<5> GhostTrackVertexFinder::mergeVertices(
00801                                         const CachingVertex<5> &vertex1,
00802                                         const CachingVertex<5> &vertex2,
00803                                         const FinderInfo &info,
00804                                         bool isPrimary) const
00805 {
00806         VertexTrackFactory<5> vertexTrackFactory;
00807 
00808         VertexState state = stateMean(vertex1.vertexState(),
00809                                       vertex2.vertexState());
00810         std::vector<RefCountedVertexTrack> linTracks;
00811         VtxTrackIs isGhostTrack(info.ghostTrack);
00812         RefCountedVertexTrack vtxGhostTrack;
00813 
00814         mergeTrackHelper(vertex1.tracks(), linTracks, state,
00815                          isGhostTrack, vtxGhostTrack, vertexTrackFactory);
00816         mergeTrackHelper(vertex2.tracks(), linTracks, state,
00817                          isGhostTrack, vtxGhostTrack, vertexTrackFactory);
00818 
00819         if (vtxGhostTrack &&
00820             (fitType_ == kAlwaysWithGhostTrack || linTracks.size() < 2))
00821                 linTracks.push_back(
00822                         vertexTrackFactory.vertexTrack(
00823                                 vtxGhostTrack->linearizedTrack(),
00824                                 vtxGhostTrack->vertexState()));
00825 
00826         try {
00827                 CachingVertex<5> vtx;
00828                 if (info.hasBeamSpot && isPrimary)
00829                         vtx = vertexFitter(true).vertex(linTracks,
00830                                                 info.beamSpot.position(),
00831                                                 info.beamSpot.error());
00832                 else
00833                         vtx = vertexFitter(isPrimary).vertex(linTracks);
00834                 if (vtx.isValid())
00835                         return vtx;
00836         } catch(const VertexException &e) {
00837                 // fit failed
00838         }
00839 
00840         return CachingVertex<5>();
00841 }
00842 
00843 bool GhostTrackVertexFinder::recursiveMerge(
00844                                 std::vector<CachingVertex<5> > &vertices,
00845                                 const FinderInfo &info) const
00846 {
00847         typedef std::pair<unsigned int, unsigned int> Indices;
00848 
00849         std::multimap<float, Indices> compatMap;
00850         unsigned int n = vertices.size();
00851         for(unsigned int i = 0; i < n; i++) {
00852                 const CachingVertex<5> &v1 = vertices[i];
00853                 for(unsigned int j = i + 1; j < n; j++) {
00854                         const CachingVertex<5> &v2 = vertices[j];
00855 
00856                         float compat = vertexCompat(v1, v2, info,
00857                                         i == 0 ? primcut_ : seccut_, seccut_);
00858 
00859                         if (compat > mergeThreshold_)
00860                                 continue;
00861 
00862                         compatMap.insert(
00863                                 std::make_pair(compat, Indices(i, j)));
00864                 }
00865         }
00866 
00867         bool changesMade = false;
00868         bool repeat = true;
00869         while(repeat) {
00870                 repeat = false;
00871                 for(std::multimap<float, Indices>::const_iterator iter =
00872                         compatMap.begin(); iter != compatMap.end(); ++iter) {
00873                         unsigned int v1 = iter->second.first;
00874                         unsigned int v2 = iter->second.second;
00875 
00876                         CachingVertex<5> newVtx =
00877                                 mergeVertices(vertices[v1], vertices[v2],
00878                                               info, v1 == 0);
00879                         if (!newVtx.isValid() ||
00880                             (v1 != 0 && fitChi2(newVtx) > maxFitChi2_))
00881                                 continue;
00882 
00883                         std::swap(vertices[v1], newVtx);
00884                         vertices.erase(vertices.begin() + v2);
00885                         n--;
00886 
00887                         std::multimap<float, Indices> newCompatMap;
00888                         for(++iter; iter != compatMap.end(); ++iter) {
00889                                 if (iter->second.first == v1 ||
00890                                     iter->second.first == v2 ||
00891                                     iter->second.second == v1 ||
00892                                     iter->second.second == v2)
00893                                         continue;
00894 
00895                                 Indices indices = iter->second;
00896                                 indices.first -= indices.first > v2;
00897                                 indices.second -= indices.second > v2;
00898 
00899                                 newCompatMap.insert(std::make_pair(
00900                                                         iter->first, indices));
00901                         }
00902                         std::swap(compatMap, newCompatMap);
00903 
00904                         for(unsigned int i = 0; i < n; i++) {
00905                                 if (i == v1)
00906                                         continue;
00907 
00908                                 const CachingVertex<5> &other = vertices[i];
00909                                 float compat = vertexCompat(
00910                                                 vertices[v1], other, info,
00911                                                 v1 == 0 ? primcut_ : seccut_,
00912                                                 i == 0 ? primcut_ : seccut_);
00913 
00914                                 if (compat > mergeThreshold_)
00915                                         continue;
00916 
00917                                 compatMap.insert(
00918                                         std::make_pair(
00919                                                 compat,
00920                                                 Indices(std::min(i, v1),
00921                                                         std::max(i, v1))));
00922                         }
00923 
00924                         changesMade = true;
00925                         repeat = true;
00926                         break;
00927                 }
00928         }
00929 
00930         return changesMade;
00931 }
00932 
00933 bool GhostTrackVertexFinder::reassignTracks(
00934                                 std::vector<CachingVertex<5> > &vertices_,
00935                                 const FinderInfo &info) const
00936 {
00937         std::vector<std::pair<RefCountedVertexTrack,
00938                               std::vector<RefCountedVertexTrack> > >
00939                                         trackBundles(vertices_.size());
00940 
00941         VtxTrackIs isGhostTrack(info.ghostTrack);
00942 
00943         bool assignmentChanged = false;
00944         for(std::vector<CachingVertex<5> >::const_iterator iter =
00945                         vertices_.begin(); iter != vertices_.end(); ++iter) {
00946                 std::vector<RefCountedVertexTrack> vtxTracks = iter->tracks();
00947 
00948                 if (vtxTracks.empty()) {
00949                         LinearizedTrackStateFactory linTrackFactory;
00950                         VertexTrackFactory<5> vertexTrackFactory;
00951 
00952                         RefCountedLinearizedTrackState linState =
00953                                 linTrackFactory.linearizedTrackState(
00954                                         iter->position(), info.ghostTrack);
00955 
00956                         trackBundles[iter - vertices_.begin()].first =
00957                                         vertexTrackFactory.vertexTrack(
00958                                                 linState, iter->vertexState());
00959                 }
00960 
00961                 for(std::vector<RefCountedVertexTrack>::const_iterator track =
00962                         vtxTracks.begin(); track != vtxTracks.end(); ++track) {
00963 
00964                         if (isGhostTrack(*track)) {
00965                                 trackBundles[iter - vertices_.begin()]
00966                                                         .first = *track;
00967                                 continue;
00968                         }
00969 
00970                         if ((*track)->weight() < 1e-3) {
00971                                 trackBundles[iter - vertices_.begin()]
00972                                                 .second.push_back(*track);
00973                                 continue;
00974                         }
00975 
00976                         unsigned int idx = iter - vertices_.begin();
00977                         double best = 1.0e9;
00978                         for(std::vector<CachingVertex<5> >::const_iterator
00979                                                 vtx = vertices_.begin();
00980                             vtx != vertices_.end(); ++vtx) {
00981                                 if (info.pred.lambda(vtx->position()) <
00982                                     info.pred.lambda(vertices_[0].position()))
00983                                         continue;
00984 
00985                                 double compat =
00986                                         sqr(trackVertexCompat(*vtx, *track));
00987 
00988                                 compat /= (vtx == vertices_.begin()) ?
00989                                                         primcut_ : seccut_;
00990 
00991                                 if (compat < best) {
00992                                         best = compat;
00993                                         idx = vtx - vertices_.begin();
00994                                 }
00995                         }
00996 
00997                         if ((int)idx != iter - vertices_.begin())
00998                                 assignmentChanged = true;
00999 
01000                         trackBundles[idx].second.push_back(*track);
01001                 }
01002         }
01003 
01004         if (!assignmentChanged)
01005                 return false;
01006 
01007         VertexTrackFactory<5> vertexTrackFactory;
01008         std::vector<CachingVertex<5> > vertices;
01009         vertices.reserve(vertices_.size());
01010 
01011         for(std::vector<CachingVertex<5> >::const_iterator iter =
01012                         vertices_.begin(); iter != vertices_.end(); ++iter) {
01013                 const std::vector<RefCountedVertexTrack> &tracks =
01014                                 trackBundles[iter - vertices_.begin()].second;
01015                 if (tracks.empty())
01016                         continue;
01017 
01018                 CachingVertex<5> vtx;
01019 
01020                 if (tracks.size() == 1) {
01021                         const TransientTrack &track =
01022                                 tracks[0]->linearizedTrack()->track();
01023 
01024                         int idx = -1;
01025                         for(std::vector<GhostTrackState>::const_iterator iter =
01026                                                         info.states.begin();
01027                             iter != info.states.end(); ++iter) {
01028                                 if (iter->isTrack() && iter->track() == track) {
01029                                         idx = iter - info.states.begin();
01030                                         break;
01031                                 }
01032                         }
01033                         if (idx < 0)
01034                                 continue;
01035 
01036                         vtx = vertexAtState(info.ghostTrack, info.pred,
01037                                             info.states[idx]);
01038                         if (!vtx.isValid())
01039                                 continue;
01040                 } else {
01041                         std::vector<RefCountedVertexTrack> linTracks =
01042                                 relinearizeTracks(tracks, iter->vertexState());
01043 
01044                         if (fitType_ == kAlwaysWithGhostTrack)
01045                                 linTracks.push_back(relinearizeTrack(
01046                                         trackBundles[iter - vertices_.begin()].first,
01047                                         iter->vertexState(), vertexTrackFactory));
01048 
01049                         bool primary = iter == vertices_.begin();
01050                         try {
01051                                 if (primary && info.hasBeamSpot)
01052                                         vtx = vertexFitter(true).vertex(
01053                                                 linTracks,
01054                                                 info.beamSpot.position(),
01055                                                 info.beamSpot.error());
01056                                 else
01057                                         vtx = vertexFitter(primary).vertex(
01058                                                                 linTracks);
01059                         } catch(const VertexException &e) {
01060                                 // fit failed;
01061                         }
01062                         if (!vtx.isValid())
01063                                 return false;
01064                 }
01065 
01066                 vertices.push_back(vtx);
01067         };
01068 
01069         std::swap(vertices_, vertices);
01070         return true;
01071 }
01072 
01073 void GhostTrackVertexFinder::refitGhostTrack(
01074                                 std::vector<CachingVertex<5> > &vertices,
01075                                 FinderInfo &info) const
01076 {
01077         VtxTrackIs isGhostTrack(info.ghostTrack);
01078 
01079         std::vector<GhostTrackState> states;
01080         std::vector<unsigned int> oldStates;
01081         oldStates.reserve(info.states.size());
01082 
01083         for(std::vector<CachingVertex<5> >::const_iterator iter =
01084                         vertices.begin(); iter != vertices.end(); ++iter) {
01085                 std::vector<RefCountedVertexTrack> vtxTracks = iter->tracks();
01086 
01087                 oldStates.clear();
01088                 for(std::vector<RefCountedVertexTrack>::const_iterator track =
01089                         vtxTracks.begin(); track != vtxTracks.end(); ++track) {
01090 
01091                         if (isGhostTrack(*track) || (*track)->weight() < 1e-3)
01092                                 continue;
01093 
01094                         const TransientTrack &tt =
01095                                         (*track)->linearizedTrack()->track();
01096 
01097                         int idx = -1;
01098                         for(std::vector<GhostTrackState>::const_iterator iter =
01099                                                         info.states.begin();
01100                             iter != info.states.end(); ++iter) {
01101                                 if (iter->isTrack() && iter->track() == tt) {
01102                                         idx = iter - info.states.begin();
01103                                         break;
01104                                 }
01105                         }
01106 
01107                         if (idx >= 0)
01108                                 oldStates.push_back(idx);
01109                 }
01110 
01111                 if (oldStates.size() == 1)
01112                         states.push_back(info.states[oldStates[0]]);
01113                 else
01114                         states.push_back(GhostTrackState(iter->vertexState()));
01115         }
01116 
01117         KalmanGhostTrackUpdater updater;
01118         SequentialGhostTrackFitter fitter;
01119         double ndof, chi2;
01120         info.pred = fitter.fit(updater, info.prior, states, ndof, chi2);
01121         TransientTrack ghostTrack = transientGhostTrack(info.pred, info.field);
01122 
01123         std::swap(info.states, states);
01124         states.clear();
01125 
01126         std::vector<CachingVertex<5> > newVertices;
01127         for(std::vector<CachingVertex<5> >::const_iterator iter =
01128                         vertices.begin(); iter != vertices.end(); ++iter) {
01129                 std::vector<RefCountedVertexTrack> vtxTracks = iter->tracks();
01130 
01131                 int idx = -1;
01132                 bool redo = false;
01133                 for(std::vector<RefCountedVertexTrack>::iterator track =
01134                         vtxTracks.begin(); track != vtxTracks.end(); ++track) {
01135 
01136                         if (isGhostTrack(*track)) {
01137                                 LinearizedTrackStateFactory linTrackFactory;
01138                                 VertexTrackFactory<5> vertexTrackFactory;
01139 
01140                                 *track = vertexTrackFactory.vertexTrack(
01141                                         linTrackFactory.linearizedTrackState(
01142                                                 iter->position(), ghostTrack),
01143                                                 iter->vertexState());
01144                                 redo = true;
01145                                 continue;
01146                         }
01147 
01148                         const TransientTrack &tt =
01149                                         (*track)->linearizedTrack()->track();
01150 
01151                         if (idx >= 0) {
01152                                 idx = -1;
01153                                 break;
01154                         }
01155 
01156                         for(std::vector<GhostTrackState>::const_iterator it =
01157                                                         info.states.begin();
01158                             it != info.states.end(); ++it) {
01159                                 if (!it->isTrack())
01160                                         continue;
01161 
01162                                 if (it->track() == tt) {
01163                                         idx = it - info.states.begin();
01164                                         break;
01165                                 }
01166                         }
01167                 }
01168 
01169                 if (idx >= 0) {
01170                         CachingVertex<5> vtx =
01171                                 vertexAtState(ghostTrack, info.pred,
01172                                               info.states[idx]);
01173                         if (vtx.isValid())
01174                                 newVertices.push_back(vtx);
01175                 } else if (redo) {
01176                         bool primary = iter == vertices.begin();
01177                         CachingVertex<5> vtx;
01178                         if (primary && info.hasBeamSpot)
01179                                 vtx = vertexFitter(true).vertex(
01180                                                 vtxTracks,
01181                                                 info.beamSpot.position(),
01182                                                 info.beamSpot.error());
01183                         else
01184                                 vtx = vertexFitter(primary).vertex(vtxTracks);
01185                         if (vtx.isValid())
01186                                 newVertices.push_back(vtx);
01187                 } else
01188                         newVertices.push_back(*iter);
01189         }
01190 
01191         std::swap(newVertices, vertices);
01192         info.ghostTrack = ghostTrack;
01193 }
01194 
01195 // implementation
01196 
01197 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(
01198                                 const GhostTrack &ghostTrack,
01199                                 const CachingVertex<5> &primary,
01200                                 const BeamSpot &beamSpot,
01201                                 bool hasBeamSpot, bool hasPrimaries) const
01202 {
01203         FinderInfo info(primary, ghostTrack, beamSpot,
01204                         hasBeamSpot, hasPrimaries);
01205 
01206         std::vector<TransientVertex> result;
01207         if (info.states.empty())
01208                 return result;
01209 
01210         info.field = info.states[0].track().field();
01211         info.ghostTrack = transientGhostTrack(info.pred, info.field);
01212 
01213         std::vector<CachingVertex<5> > vertices = initialVertices(info);
01214         if (primary.isValid()) {
01215                 vertices.push_back(primary);
01216                 if (vertices.size() > 1)
01217                         std::swap(vertices.front(), vertices.back());
01218         }
01219 
01220         unsigned int reassigned = 0;
01221         while(reassigned < 3) {
01222                 if (vertices.size() < 2)
01223                         break;
01224 
01225 #ifdef DEBUG
01226                 for(std::vector<CachingVertex<5> >::const_iterator iter =
01227                         vertices.begin(); iter != vertices.end(); ++iter)
01228 
01229                         debugVertex(*iter, ghostTrack.prediction());
01230 
01231                 std::cout << "----- recursive merging: ---------" << std::endl;
01232 #endif
01233 
01234                 bool changed = recursiveMerge(vertices, info);
01235                 if ((!changed && !reassigned) || vertices.size() < 2)
01236                         break;
01237                 if (changed)
01238                         reassigned = 0;
01239 
01240                 if (fitType_ == kRefitGhostTrackWithVertices) {
01241                         refitGhostTrack(vertices, info);
01242                         changed = true;
01243                 }
01244 
01245                 try {
01246 #ifdef DEBUG
01247                         for(std::vector<CachingVertex<5> >::const_iterator
01248                                                 iter = vertices.begin();
01249                             iter != vertices.end(); ++iter)
01250                                 debugVertex(*iter, ghostTrack.prediction());
01251                         std::cout << "----- reassignment: ---------" << std::endl;
01252 #endif
01253                         if (reassignTracks(vertices, info)) {
01254                                 reassigned++;
01255                                 changed = true;
01256                         } else
01257                                 reassigned = 0;
01258                 } catch(const VertexException &e) {
01259                         // just keep vertices as they are
01260                 }
01261 
01262                 if (!changed)
01263                         break;
01264         }
01265 
01266         for(std::vector<CachingVertex<5> >::const_iterator iter =
01267                         vertices.begin(); iter != vertices.end(); ++iter) {
01268                 std::vector<RefCountedVertexTrack> tracks = iter->tracks();
01269                 std::vector<RefCountedVertexTrack> newTracks;
01270                 newTracks.reserve(tracks.size());
01271 
01272                 std::remove_copy_if(tracks.begin(), tracks.end(),
01273                                     std::back_inserter(newTracks),
01274                                     VtxTrackIs(info.ghostTrack));
01275 
01276                 if (newTracks.empty())
01277                         continue;
01278 
01279                 CachingVertex<5> vtx(iter->vertexState(), newTracks,
01280                                      iter->totalChiSquared());
01281                 result.push_back(vtx);
01282         }
01283 
01284         return result;
01285 }
01286