File indexing completed on 2024-04-06 12:29:08
0001 #include <algorithm>
0002 #include <iterator>
0003 #include <memory>
0004
0005 #include <map>
0006 #include <set>
0007 #include <vector>
0008
0009 #include <Math/SVector.h>
0010 #include <Math/SMatrix.h>
0011 #include <Math/MatrixFunctions.h>
0012
0013 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0014 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0015 #include "DataFormats/GeometryCommonDetAlgo/interface/GlobalError.h"
0016
0017 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0018
0019 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0020 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0021 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0022 #include "TrackingTools/GeomPropagators/interface/AnalyticalImpactPointExtrapolator.h"
0023 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0024 #include "TrackingTools/TransientTrack/interface/TransientTrackFromFTSFactory.h"
0025
0026 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0027 #include "RecoVertex/VertexPrimitives/interface/VertexFitter.h"
0028 #include "RecoVertex/VertexPrimitives/interface/VertexState.h"
0029 #include "RecoVertex/VertexPrimitives/interface/CachingVertex.h"
0030 #include "RecoVertex/VertexPrimitives/interface/LinearizedTrackState.h"
0031 #include "RecoVertex/VertexPrimitives/interface/ConvertError.h"
0032 #include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
0033 #include "RecoVertex/VertexTools/interface/LinearizedTrackStateFactory.h"
0034 #include "RecoVertex/VertexTools/interface/VertexTrackFactory.h"
0035 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
0036 #include "RecoVertex/VertexTools/interface/GeometricAnnealing.h"
0037 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
0038 #include "RecoVertex/AdaptiveVertexFit/interface/AdaptiveVertexFitter.h"
0039
0040 #include "RecoVertex/GhostTrackFitter/interface/GhostTrack.h"
0041 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackState.h"
0042 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackFitter.h"
0043 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackPrediction.h"
0044 #include "RecoVertex/GhostTrackFitter/interface/SequentialGhostTrackFitter.h"
0045 #include "RecoVertex/GhostTrackFitter/interface/KalmanGhostTrackUpdater.h"
0046
0047 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackVertexFinder.h"
0048
0049 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0050
0051
0052
0053 #ifdef DEBUG
0054 #include "RecoBTag/SecondaryVertex/interface/TrackKinematics.h"
0055 #endif
0056
0057 using namespace reco;
0058
0059 namespace {
0060 using namespace ROOT::Math;
0061
0062 typedef SVector<double, 3> Vector3;
0063
0064 typedef SMatrix<double, 3, 3, MatRepSym<double, 3> > Matrix3S;
0065 typedef SMatrix<double, 2, 2, MatRepSym<double, 2> > Matrix2S;
0066 typedef SMatrix<double, 2, 3> Matrix23;
0067
0068 typedef ReferenceCountingPointer<VertexTrack<5> > RefCountedVertexTrack;
0069 typedef ReferenceCountingPointer<LinearizedTrackState<5> > RefCountedLinearizedTrackState;
0070
0071 struct VtxTrackIs {
0072 VtxTrackIs(const TransientTrack &track) : track(track) {}
0073 bool operator()(const RefCountedVertexTrack &vtxTrack) const {
0074 return vtxTrack->linearizedTrack()->track() == track;
0075 }
0076
0077 const TransientTrack &track;
0078 };
0079
0080 inline Vector3 conv(const GlobalVector &vec) {
0081 Vector3 result;
0082 result[0] = vec.x();
0083 result[1] = vec.y();
0084 result[2] = vec.z();
0085 return result;
0086 }
0087
0088 inline double sqr(double arg) { return arg * arg; }
0089 }
0090
0091 struct GhostTrackVertexFinder::FinderInfo {
0092 FinderInfo(const CachingVertex<5> &primaryVertex,
0093 const GhostTrack &ghostTrack,
0094 const BeamSpot &beamSpot,
0095 bool hasBeamSpot,
0096 bool hasPrimaries)
0097 : primaryVertex(primaryVertex),
0098 pred(ghostTrack.prediction()),
0099 prior(ghostTrack.prior()),
0100 states(ghostTrack.states()),
0101 beamSpot(beamSpot),
0102 hasBeamSpot(hasBeamSpot),
0103 hasPrimaries(hasPrimaries),
0104 field(nullptr) {}
0105
0106 const CachingVertex<5> &primaryVertex;
0107 GhostTrackPrediction pred;
0108 GhostTrackPrediction prior;
0109 std::vector<GhostTrackState> states;
0110 VertexState beamSpot;
0111 bool hasBeamSpot;
0112 bool hasPrimaries;
0113 const MagneticField *field;
0114 TransientTrack ghostTrack;
0115 };
0116
0117 static TransientTrack transientGhostTrack(const GhostTrackPrediction &pred, const MagneticField *field) {
0118 return TransientTrackFromFTSFactory().build(pred.fts(field));
0119 }
0120
0121 static double vtxErrorLong(const GlobalError &error, const GlobalVector &dir) {
0122 return ROOT::Math::Similarity(conv(dir), error.matrix());
0123 }
0124
0125 static GlobalPoint vtxMean(const GlobalPoint &p1, const GlobalError &e1, const GlobalPoint &p2, const GlobalError &e2) {
0126 GlobalVector diff = p2 - p1;
0127
0128 if UNLIKELY (p1 == p2) {
0129 edm::LogWarning("GhostTrackVertexFinder") << "Averaging with self at " << p1;
0130 return p1;
0131 }
0132
0133 double err1 = vtxErrorLong(e1, diff);
0134 double err2 = vtxErrorLong(e2, diff);
0135
0136 double weight = err1 / (err1 + err2);
0137
0138 return p1 + weight * diff;
0139 }
0140
0141 static VertexState stateMean(const VertexState &v1, const VertexState &v2) {
0142 return VertexState(vtxMean(v1.position(), v1.error(), v2.position(), v2.error()), v1.error() + v2.error());
0143 }
0144
0145 static bool covarianceUpdate(
0146 Matrix3S &cov, const Vector3 &residual, const Matrix3S &error, double &chi2, double theta, double phi) {
0147 using namespace ROOT::Math;
0148
0149 Matrix23 jacobian;
0150 jacobian(0, 0) = std::cos(phi) * std::cos(theta);
0151 jacobian(0, 1) = std::sin(phi) * std::cos(theta);
0152 jacobian(0, 2) = -std::sin(theta);
0153 jacobian(1, 0) = -std::sin(phi);
0154 jacobian(1, 1) = std::cos(phi);
0155
0156 Matrix2S measErr = Similarity(jacobian, error);
0157 Matrix2S combErr = Similarity(jacobian, cov) + measErr;
0158 if (!measErr.Invert() || !combErr.Invert())
0159 return false;
0160
0161 cov -= SimilarityT(jacobian * cov, combErr);
0162 chi2 += Similarity(jacobian * residual, measErr);
0163
0164 return true;
0165 }
0166
0167 static CachingVertex<5> vertexAtState(const TransientTrack &ghostTrack,
0168 const GhostTrackPrediction &pred,
0169 const GhostTrackState &state) {
0170 LinearizedTrackStateFactory linTrackFactory;
0171 VertexTrackFactory<5> vertexTrackFactory;
0172
0173 if (!state.isValid())
0174 return CachingVertex<5>();
0175
0176 GlobalPoint pca1 = pred.position(state.lambda());
0177 GlobalError err1 = pred.positionError(state.lambda());
0178
0179 GlobalPoint pca2 = state.globalPosition();
0180 GlobalError err2 = state.cartesianError();
0181
0182 GlobalPoint point = vtxMean(pca1, err1, pca2, err2);
0183
0184 const TransientTrack &recTrack = state.track();
0185
0186 RefCountedLinearizedTrackState linState[2] = {linTrackFactory.linearizedTrackState(point, ghostTrack),
0187 linTrackFactory.linearizedTrackState(point, recTrack)};
0188 if (!linState[0]->isValid() || !linState[1]->isValid())
0189 return CachingVertex<5>();
0190
0191 Matrix3S cov = SMatrixIdentity();
0192 cov *= 10000;
0193 double chi2 = 0.;
0194 if (!covarianceUpdate(cov,
0195 conv(pca1 - point),
0196 err1.matrix(),
0197 chi2,
0198 linState[0]->predictedStateParameters()[1],
0199 linState[0]->predictedStateParameters()[2]) ||
0200 !covarianceUpdate(cov,
0201 conv(pca2 - point),
0202 err2.matrix(),
0203 chi2,
0204 linState[1]->predictedStateParameters()[1],
0205 linState[1]->predictedStateParameters()[2]))
0206 return CachingVertex<5>();
0207
0208 GlobalError error(cov);
0209 VertexState vtxState(point, error);
0210
0211 std::vector<RefCountedVertexTrack> linTracks(2);
0212 linTracks[0] = vertexTrackFactory.vertexTrack(linState[0], vtxState);
0213 linTracks[1] = vertexTrackFactory.vertexTrack(linState[1], vtxState);
0214
0215 return CachingVertex<5>(point, error, linTracks, chi2);
0216 }
0217
0218 static RefCountedVertexTrack relinearizeTrack(const RefCountedVertexTrack &track,
0219 const VertexState &state,
0220 const VertexTrackFactory<5> factory) {
0221 RefCountedLinearizedTrackState linTrack = track->linearizedTrack();
0222 linTrack = linTrack->stateWithNewLinearizationPoint(state.position());
0223 return factory.vertexTrack(linTrack, state);
0224 }
0225
0226 static std::vector<RefCountedVertexTrack> relinearizeTracks(const std::vector<RefCountedVertexTrack> &tracks,
0227 const VertexState &state) {
0228 VertexTrackFactory<5> vertexTrackFactory;
0229
0230 std::vector<RefCountedVertexTrack> finalTracks;
0231 finalTracks.reserve(tracks.size());
0232
0233 for (std::vector<RefCountedVertexTrack>::const_iterator iter = tracks.begin(); iter != tracks.end(); ++iter)
0234 finalTracks.push_back(relinearizeTrack(*iter, state, vertexTrackFactory));
0235
0236 return finalTracks;
0237 }
0238
0239 double GhostTrackVertexFinder::vertexCompat(
0240 const CachingVertex<5> &vtx1, const CachingVertex<5> &vtx2, const FinderInfo &info, double scale1, double scale2) {
0241 Vector3 diff = conv(vtx2.position() - vtx1.position());
0242 Matrix3S cov = scale1 * vtx1.error().matrix() + scale2 * vtx2.error().matrix();
0243
0244 return sqr(ROOT::Math::Dot(diff, diff)) / ROOT::Math::Similarity(cov, diff);
0245 }
0246
0247 static double trackVertexCompat(const CachingVertex<5> &vtx, const RefCountedVertexTrack &vertexTrack) {
0248 using namespace ROOT::Math;
0249
0250 TransientTrack track = vertexTrack->linearizedTrack()->track();
0251 GlobalPoint point = vtx.position();
0252 AnalyticalImpactPointExtrapolator extrap(track.field());
0253 TrajectoryStateOnSurface tsos = extrap.extrapolate(track.impactPointState(), point);
0254
0255 if (!tsos.isValid())
0256 return 1.0e6;
0257
0258 GlobalPoint point1 = vtx.position();
0259 GlobalPoint point2 = tsos.globalPosition();
0260 Vector3 dir = conv(point2 - point1);
0261 Matrix3S error = vtx.error().matrix() + tsos.cartesianError().matrix().Sub<Matrix3S>(0, 0);
0262 if (!error.Invert())
0263 return 1.0e6;
0264
0265 return ROOT::Math::Similarity(error, dir);
0266 }
0267
0268 GhostTrackVertexFinder::GhostTrackVertexFinder()
0269 : maxFitChi2_(10.0), mergeThreshold_(3.0), primcut_(2.0), seccut_(5.0), fitType_(kAlwaysWithGhostTrack) {}
0270
0271 GhostTrackVertexFinder::GhostTrackVertexFinder(
0272 double maxFitChi2, double mergeThreshold, double primcut, double seccut, FitType fitType)
0273 : maxFitChi2_(maxFitChi2), mergeThreshold_(mergeThreshold), primcut_(primcut), seccut_(seccut), fitType_(fitType) {}
0274
0275 GhostTrackVertexFinder::~GhostTrackVertexFinder() {}
0276
0277 GhostTrackFitter &GhostTrackVertexFinder::ghostTrackFitter() const {
0278 if (!ghostTrackFitter_.get())
0279 ghostTrackFitter_ = std::make_unique<GhostTrackFitter>();
0280
0281 return *ghostTrackFitter_;
0282 }
0283
0284 VertexFitter<5> &GhostTrackVertexFinder::vertexFitter(bool primary) const {
0285 std::unique_ptr<VertexFitter<5> > *ptr = primary ? &primVertexFitter_ : &secVertexFitter_;
0286 if (!ptr->get())
0287 *ptr = std::make_unique<AdaptiveVertexFitter>(GeometricAnnealing(primary ? primcut_ : seccut_));
0288
0289 return **ptr;
0290 }
0291
0292 #ifdef DEBUG
0293 static void debugTrack(const TransientTrack &track, const TransientVertex *vtx = 0) {
0294 const FreeTrajectoryState &fts = track.initialFreeState();
0295 GlobalVector mom = fts.momentum();
0296 std::cout << "\tpt = " << mom.perp() << ", eta = " << mom.eta() << ", phi = " << mom.phi()
0297 << ", charge = " << fts.charge();
0298 if (vtx && vtx->hasTrackWeight())
0299 std::cout << ", w = " << vtx->trackWeight(track) << std::endl;
0300 else
0301 std::cout << std::endl;
0302 }
0303
0304 static void debugTracks(const std::vector<TransientTrack> &tracks, const TransientVertex *vtx = 0) {
0305 TrackKinematics kin;
0306 for (std::vector<TransientTrack>::const_iterator iter = tracks.begin(); iter != tracks.end(); ++iter) {
0307 debugTrack(*iter, vtx);
0308 try {
0309 TrackBaseRef track = iter->trackBaseRef();
0310 kin.add(*track);
0311 } catch (exception const &e) {
0312
0313 }
0314 }
0315
0316 std::cout << "mass = " << kin.vectorSum().M() << std::endl;
0317 }
0318
0319 static void debugTracks(const std::vector<RefCountedVertexTrack> &tracks) {
0320 std::vector<TransientTrack> tracks_;
0321 for (std::vector<RefCountedVertexTrack>::const_iterator iter = tracks.begin(); iter != tracks.end(); ++iter) {
0322 std::cout << "+ " << (*iter)->linearizedTrack()->linearizationPoint() << std::endl;
0323 tracks_.push_back((*iter)->linearizedTrack()->track());
0324 }
0325 debugTracks(tracks_);
0326 }
0327
0328 static void debugVertex(const TransientVertex &vtx, const GhostTrackPrediction &pred) {
0329 double origin = 0.;
0330 std::cout << vtx.position() << "-> " << vtx.totalChiSquared() << " with " << vtx.degreesOfFreedom() << std::endl;
0331
0332 double err = std::sqrt(vtxErrorLong(vtx.positionError(), pred.direction()) / pred.rho2());
0333 std::cout << "* " << (pred.lambda(vtx.position()) * pred.rho() - origin) << " +-" << err << " => "
0334 << ((pred.lambda(vtx.position()) * pred.rho() - origin) / err) << std::endl;
0335
0336 std::vector<TransientTrack> tracks = vtx.originalTracks();
0337 debugTracks(tracks, &vtx);
0338 }
0339 #endif
0340
0341 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const Vertex &primaryVertex,
0342 const GlobalVector &direction,
0343 double coneRadius,
0344 const std::vector<TransientTrack> &tracks) const {
0345 return vertices(RecoVertex::convertPos(primaryVertex.position()),
0346 RecoVertex::convertError(primaryVertex.error()),
0347 direction,
0348 coneRadius,
0349 tracks);
0350 }
0351
0352 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const Vertex &primaryVertex,
0353 const GlobalVector &direction,
0354 double coneRadius,
0355 const reco::BeamSpot &beamSpot,
0356 const std::vector<TransientTrack> &tracks) const {
0357 return vertices(RecoVertex::convertPos(primaryVertex.position()),
0358 RecoVertex::convertError(primaryVertex.error()),
0359 direction,
0360 coneRadius,
0361 beamSpot,
0362 tracks);
0363 }
0364
0365 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const Vertex &primaryVertex,
0366 const GlobalVector &direction,
0367 double coneRadius,
0368 const reco::BeamSpot &beamSpot,
0369 const std::vector<TransientTrack> &primaries,
0370 const std::vector<TransientTrack> &tracks) const {
0371 return vertices(RecoVertex::convertPos(primaryVertex.position()),
0372 RecoVertex::convertError(primaryVertex.error()),
0373 direction,
0374 coneRadius,
0375 beamSpot,
0376 primaries,
0377 tracks);
0378 }
0379
0380 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const GlobalPoint &primaryPosition,
0381 const GlobalError &primaryError,
0382 const GlobalVector &direction,
0383 double coneRadius,
0384 const std::vector<TransientTrack> &tracks) const {
0385 GhostTrack ghostTrack = ghostTrackFitter().fit(primaryPosition, primaryError, direction, coneRadius, tracks);
0386
0387 CachingVertex<5> primary(primaryPosition, primaryError, std::vector<RefCountedVertexTrack>(), 0.);
0388
0389 std::vector<TransientVertex> result = vertices(ghostTrack, primary);
0390
0391 #ifdef DEBUG
0392 for (std::vector<TransientVertex>::const_iterator iter = result.begin(); iter != result.end(); ++iter)
0393 debugVertex(*iter, ghostTrack.prediction());
0394 #endif
0395
0396 return result;
0397 }
0398
0399 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const GlobalPoint &primaryPosition,
0400 const GlobalError &primaryError,
0401 const GlobalVector &direction,
0402 double coneRadius,
0403 const BeamSpot &beamSpot,
0404 const std::vector<TransientTrack> &tracks) const {
0405 GhostTrack ghostTrack = ghostTrackFitter().fit(primaryPosition, primaryError, direction, coneRadius, tracks);
0406
0407 CachingVertex<5> primary(primaryPosition, primaryError, std::vector<RefCountedVertexTrack>(), 0.);
0408
0409 std::vector<TransientVertex> result = vertices(ghostTrack, primary, beamSpot, true);
0410
0411 #ifdef DEBUG
0412 for (std::vector<TransientVertex>::const_iterator iter = result.begin(); iter != result.end(); ++iter)
0413 debugVertex(*iter, ghostTrack.prediction());
0414 #endif
0415
0416 return result;
0417 }
0418
0419 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const GlobalPoint &primaryPosition,
0420 const GlobalError &primaryError,
0421 const GlobalVector &direction,
0422 double coneRadius,
0423 const BeamSpot &beamSpot,
0424 const std::vector<TransientTrack> &primaries,
0425 const std::vector<TransientTrack> &tracks) const {
0426 GhostTrack ghostTrack = ghostTrackFitter().fit(primaryPosition, primaryError, direction, coneRadius, tracks);
0427
0428 std::vector<RefCountedVertexTrack> primaryVertexTracks;
0429 if (!primaries.empty()) {
0430 LinearizedTrackStateFactory linTrackFactory;
0431 VertexTrackFactory<5> vertexTrackFactory;
0432
0433 VertexState state(primaryPosition, primaryError);
0434
0435 for (std::vector<TransientTrack>::const_iterator iter = primaries.begin(); iter != primaries.end(); ++iter) {
0436 RefCountedLinearizedTrackState linState = linTrackFactory.linearizedTrackState(primaryPosition, *iter);
0437
0438 primaryVertexTracks.push_back(vertexTrackFactory.vertexTrack(linState, state));
0439 }
0440 }
0441
0442 CachingVertex<5> primary(primaryPosition, primaryError, primaryVertexTracks, 0.);
0443
0444 std::vector<TransientVertex> result = vertices(ghostTrack, primary, beamSpot, true, true);
0445
0446 #ifdef DEBUG
0447 for (std::vector<TransientVertex>::const_iterator iter = result.begin(); iter != result.end(); ++iter)
0448 debugVertex(*iter, ghostTrack.prediction());
0449 #endif
0450
0451 return result;
0452 }
0453
0454 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const GlobalPoint &primaryPosition,
0455 const GlobalError &primaryError,
0456 const GhostTrack &ghostTrack) const {
0457 CachingVertex<5> primary(primaryPosition, primaryError, std::vector<RefCountedVertexTrack>(), 0.);
0458
0459 std::vector<TransientVertex> result = vertices(ghostTrack, primary);
0460
0461 #ifdef DEBUG
0462 for (std::vector<TransientVertex>::const_iterator iter = result.begin(); iter != result.end(); ++iter)
0463 debugVertex(*iter, ghostTrack.prediction());
0464 #endif
0465
0466 return result;
0467 }
0468
0469 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const GlobalPoint &primaryPosition,
0470 const GlobalError &primaryError,
0471 const BeamSpot &beamSpot,
0472 const GhostTrack &ghostTrack) const {
0473 CachingVertex<5> primary(primaryPosition, primaryError, std::vector<RefCountedVertexTrack>(), 0.);
0474
0475 std::vector<TransientVertex> result = vertices(ghostTrack, primary, beamSpot, true);
0476
0477 #ifdef DEBUG
0478 for (std::vector<TransientVertex>::const_iterator iter = result.begin(); iter != result.end(); ++iter)
0479 debugVertex(*iter, ghostTrack.prediction());
0480 #endif
0481
0482 return result;
0483 }
0484
0485 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const GlobalPoint &primaryPosition,
0486 const GlobalError &primaryError,
0487 const BeamSpot &beamSpot,
0488 const std::vector<TransientTrack> &primaries,
0489 const GhostTrack &ghostTrack) const {
0490 std::vector<RefCountedVertexTrack> primaryVertexTracks;
0491 if (!primaries.empty()) {
0492 LinearizedTrackStateFactory linTrackFactory;
0493 VertexTrackFactory<5> vertexTrackFactory;
0494
0495 VertexState state(primaryPosition, primaryError);
0496
0497 for (std::vector<TransientTrack>::const_iterator iter = primaries.begin(); iter != primaries.end(); ++iter) {
0498 RefCountedLinearizedTrackState linState = linTrackFactory.linearizedTrackState(primaryPosition, *iter);
0499
0500 primaryVertexTracks.push_back(vertexTrackFactory.vertexTrack(linState, state));
0501 }
0502 }
0503
0504 CachingVertex<5> primary(primaryPosition, primaryError, primaryVertexTracks, 0.);
0505
0506 std::vector<TransientVertex> result = vertices(ghostTrack, primary, beamSpot, true, true);
0507
0508 #ifdef DEBUG
0509 for (std::vector<TransientVertex>::const_iterator iter = result.begin(); iter != result.end(); ++iter)
0510 debugVertex(*iter, ghostTrack.prediction());
0511 #endif
0512
0513 return result;
0514 }
0515
0516 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const Vertex &primaryVertex,
0517 const GhostTrack &ghostTrack) const {
0518 return vertices(
0519 RecoVertex::convertPos(primaryVertex.position()), RecoVertex::convertError(primaryVertex.error()), ghostTrack);
0520 }
0521
0522 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const Vertex &primaryVertex,
0523 const BeamSpot &beamSpot,
0524 const GhostTrack &ghostTrack) const {
0525 return vertices(RecoVertex::convertPos(primaryVertex.position()),
0526 RecoVertex::convertError(primaryVertex.error()),
0527 beamSpot,
0528 ghostTrack);
0529 }
0530
0531 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const Vertex &primaryVertex,
0532 const BeamSpot &beamSpot,
0533 const std::vector<TransientTrack> &primaries,
0534 const GhostTrack &ghostTrack) const {
0535 return vertices(RecoVertex::convertPos(primaryVertex.position()),
0536 RecoVertex::convertError(primaryVertex.error()),
0537 beamSpot,
0538 primaries,
0539 ghostTrack);
0540 }
0541
0542 static GhostTrackPrediction dummyPrediction(const Vertex &primaryVertex, const Track &ghostTrack) {
0543 return GhostTrackPrediction(RecoVertex::convertPos(primaryVertex.position()),
0544 RecoVertex::convertError(primaryVertex.error()),
0545 GlobalVector(ghostTrack.px(), ghostTrack.py(), ghostTrack.pz()),
0546 0.05);
0547 }
0548
0549 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const Vertex &primaryVertex,
0550 const Track &ghostTrack,
0551 const std::vector<TransientTrack> &tracks,
0552 const std::vector<float> &weights) const {
0553 GhostTrack ghostTrack_(ghostTrack,
0554 tracks,
0555 weights,
0556 dummyPrediction(primaryVertex, ghostTrack),
0557 RecoVertex::convertPos(primaryVertex.position()),
0558 true);
0559
0560 CachingVertex<5> primary(RecoVertex::convertPos(primaryVertex.position()),
0561 RecoVertex::convertError(primaryVertex.error()),
0562 std::vector<RefCountedVertexTrack>(),
0563 0.);
0564
0565 std::vector<TransientVertex> result = vertices(ghostTrack_, primary);
0566
0567 #ifdef DEBUG
0568 for (std::vector<TransientVertex>::const_iterator iter = result.begin(); iter != result.end(); ++iter)
0569 debugVertex(*iter, ghostTrack_.prediction());
0570 #endif
0571
0572 return result;
0573 }
0574
0575 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const Vertex &primaryVertex,
0576 const Track &ghostTrack,
0577 const BeamSpot &beamSpot,
0578 const std::vector<TransientTrack> &tracks,
0579 const std::vector<float> &weights) const {
0580 GhostTrack ghostTrack_(ghostTrack,
0581 tracks,
0582 weights,
0583 dummyPrediction(primaryVertex, ghostTrack),
0584 RecoVertex::convertPos(primaryVertex.position()),
0585 true);
0586
0587 CachingVertex<5> primary(RecoVertex::convertPos(primaryVertex.position()),
0588 RecoVertex::convertError(primaryVertex.error()),
0589 std::vector<RefCountedVertexTrack>(),
0590 0.);
0591
0592 std::vector<TransientVertex> result = vertices(ghostTrack_, primary, beamSpot, true);
0593
0594 #ifdef DEBUG
0595 for (std::vector<TransientVertex>::const_iterator iter = result.begin(); iter != result.end(); ++iter)
0596 debugVertex(*iter, ghostTrack_.prediction());
0597 #endif
0598
0599 return result;
0600 }
0601
0602 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const Vertex &primaryVertex,
0603 const Track &ghostTrack,
0604 const BeamSpot &beamSpot,
0605 const std::vector<TransientTrack> &primaries,
0606 const std::vector<TransientTrack> &tracks,
0607 const std::vector<float> &weights) const {
0608 GhostTrack ghostTrack_(ghostTrack,
0609 tracks,
0610 weights,
0611 dummyPrediction(primaryVertex, ghostTrack),
0612 RecoVertex::convertPos(primaryVertex.position()),
0613 true);
0614
0615 std::vector<RefCountedVertexTrack> primaryVertexTracks;
0616 if (!primaries.empty()) {
0617 LinearizedTrackStateFactory linTrackFactory;
0618 VertexTrackFactory<5> vertexTrackFactory;
0619
0620 VertexState state(RecoVertex::convertPos(primaryVertex.position()),
0621 RecoVertex::convertError(primaryVertex.error()));
0622
0623 for (std::vector<TransientTrack>::const_iterator iter = primaries.begin(); iter != primaries.end(); ++iter) {
0624 RefCountedLinearizedTrackState linState = linTrackFactory.linearizedTrackState(state.position(), *iter);
0625
0626 primaryVertexTracks.push_back(vertexTrackFactory.vertexTrack(linState, state));
0627 }
0628 }
0629
0630 CachingVertex<5> primary(RecoVertex::convertPos(primaryVertex.position()),
0631 RecoVertex::convertError(primaryVertex.error()),
0632 primaryVertexTracks,
0633 0.);
0634
0635 std::vector<TransientVertex> result = vertices(ghostTrack_, primary, beamSpot, true, true);
0636
0637 #ifdef DEBUG
0638 for (std::vector<TransientVertex>::const_iterator iter = result.begin(); iter != result.end(); ++iter)
0639 debugVertex(*iter, ghostTrack_.prediction());
0640 #endif
0641
0642 return result;
0643 }
0644
0645 static double fitChi2(const CachingVertex<5> &vtx) { return vtx.totalChiSquared() / vtx.degreesOfFreedom(); }
0646
0647 std::vector<CachingVertex<5> > GhostTrackVertexFinder::initialVertices(const FinderInfo &info) const {
0648 std::vector<CachingVertex<5> > vertices;
0649 for (std::vector<GhostTrackState>::const_iterator iter = info.states.begin(); iter != info.states.end(); ++iter) {
0650 if (!iter->isValid())
0651 continue;
0652
0653 GhostTrackState state(*iter);
0654 state.linearize(info.pred);
0655
0656 if (!state.isValid())
0657 continue;
0658
0659 CachingVertex<5> vtx = vertexAtState(info.ghostTrack, info.pred, state);
0660
0661 if (vtx.isValid())
0662 vertices.push_back(vtx);
0663 }
0664
0665 return vertices;
0666 }
0667
0668 static void mergeTrackHelper(const std::vector<RefCountedVertexTrack> &tracks,
0669 std::vector<RefCountedVertexTrack> &newTracks,
0670 const VertexState &state,
0671 const VtxTrackIs &ghostTrackFinder,
0672 RefCountedVertexTrack &ghostTrack,
0673 const VertexTrackFactory<5> &factory) {
0674 for (std::vector<RefCountedVertexTrack>::const_iterator iter = tracks.begin(); iter != tracks.end(); ++iter) {
0675 bool gt = ghostTrackFinder(*iter);
0676 if (gt && ghostTrack)
0677 continue;
0678
0679 RefCountedVertexTrack track = relinearizeTrack(*iter, state, factory);
0680
0681 if (gt)
0682 ghostTrack = *iter;
0683 else
0684 newTracks.push_back(*iter);
0685 }
0686 }
0687
0688 CachingVertex<5> GhostTrackVertexFinder::mergeVertices(const CachingVertex<5> &vertex1,
0689 const CachingVertex<5> &vertex2,
0690 const FinderInfo &info,
0691 bool isPrimary) const {
0692 VertexTrackFactory<5> vertexTrackFactory;
0693
0694 VertexState state = stateMean(vertex1.vertexState(), vertex2.vertexState());
0695 std::vector<RefCountedVertexTrack> linTracks;
0696 VtxTrackIs isGhostTrack(info.ghostTrack);
0697 RefCountedVertexTrack vtxGhostTrack;
0698
0699 mergeTrackHelper(vertex1.tracks(), linTracks, state, isGhostTrack, vtxGhostTrack, vertexTrackFactory);
0700 mergeTrackHelper(vertex2.tracks(), linTracks, state, isGhostTrack, vtxGhostTrack, vertexTrackFactory);
0701
0702 if (vtxGhostTrack && (fitType_ == kAlwaysWithGhostTrack || linTracks.size() < 2))
0703 linTracks.push_back(vertexTrackFactory.vertexTrack(vtxGhostTrack->linearizedTrack(), vtxGhostTrack->vertexState()));
0704
0705 try {
0706 CachingVertex<5> vtx;
0707 if (info.hasBeamSpot && isPrimary)
0708 vtx = vertexFitter(true).vertex(linTracks, info.beamSpot.position(), info.beamSpot.error());
0709 else
0710 vtx = vertexFitter(isPrimary).vertex(linTracks);
0711 if (vtx.isValid())
0712 return vtx;
0713 } catch (const VertexException &e) {
0714
0715 }
0716
0717 return CachingVertex<5>();
0718 }
0719
0720 bool GhostTrackVertexFinder::recursiveMerge(std::vector<CachingVertex<5> > &vertices, const FinderInfo &info) const {
0721 typedef std::pair<unsigned int, unsigned int> Indices;
0722
0723 std::multimap<float, Indices> compatMap;
0724 unsigned int n = vertices.size();
0725 for (unsigned int i = 0; i < n; i++) {
0726 const CachingVertex<5> &v1 = vertices[i];
0727 for (unsigned int j = i + 1; j < n; j++) {
0728 const CachingVertex<5> &v2 = vertices[j];
0729
0730 float compat = vertexCompat(v1, v2, info, i == 0 ? primcut_ : seccut_, seccut_);
0731
0732 if (compat > mergeThreshold_)
0733 continue;
0734
0735 compatMap.insert(std::make_pair(compat, Indices(i, j)));
0736 }
0737 }
0738
0739 bool changesMade = false;
0740 bool repeat = true;
0741 while (repeat) {
0742 repeat = false;
0743 for (std::multimap<float, Indices>::const_iterator iter = compatMap.begin(); iter != compatMap.end(); ++iter) {
0744 unsigned int v1 = iter->second.first;
0745 unsigned int v2 = iter->second.second;
0746
0747 CachingVertex<5> newVtx = mergeVertices(vertices[v1], vertices[v2], info, v1 == 0);
0748 if (!newVtx.isValid() || (v1 != 0 && fitChi2(newVtx) > maxFitChi2_))
0749 continue;
0750
0751 std::swap(vertices[v1], newVtx);
0752 vertices.erase(vertices.begin() + v2);
0753 n--;
0754
0755 std::multimap<float, Indices> newCompatMap;
0756 for (++iter; iter != compatMap.end(); ++iter) {
0757 if (iter->second.first == v1 || iter->second.first == v2 || iter->second.second == v1 ||
0758 iter->second.second == v2)
0759 continue;
0760
0761 Indices indices = iter->second;
0762 indices.first -= indices.first > v2;
0763 indices.second -= indices.second > v2;
0764
0765 newCompatMap.insert(std::make_pair(iter->first, indices));
0766 }
0767 std::swap(compatMap, newCompatMap);
0768
0769 for (unsigned int i = 0; i < n; i++) {
0770 if (i == v1)
0771 continue;
0772
0773 const CachingVertex<5> &other = vertices[i];
0774 float compat =
0775 vertexCompat(vertices[v1], other, info, v1 == 0 ? primcut_ : seccut_, i == 0 ? primcut_ : seccut_);
0776
0777 if (compat > mergeThreshold_)
0778 continue;
0779
0780 compatMap.insert(std::make_pair(compat, Indices(std::min(i, v1), std::max(i, v1))));
0781 }
0782
0783 changesMade = true;
0784 repeat = true;
0785 break;
0786 }
0787 }
0788
0789 return changesMade;
0790 }
0791
0792 bool GhostTrackVertexFinder::reassignTracks(std::vector<CachingVertex<5> > &vertices_, const FinderInfo &info) const {
0793 std::vector<std::pair<RefCountedVertexTrack, std::vector<RefCountedVertexTrack> > > trackBundles(vertices_.size());
0794
0795 VtxTrackIs isGhostTrack(info.ghostTrack);
0796
0797 bool assignmentChanged = false;
0798 for (std::vector<CachingVertex<5> >::const_iterator iter = vertices_.begin(); iter != vertices_.end(); ++iter) {
0799 std::vector<RefCountedVertexTrack> vtxTracks = iter->tracks();
0800
0801 if (vtxTracks.empty()) {
0802 LinearizedTrackStateFactory linTrackFactory;
0803 VertexTrackFactory<5> vertexTrackFactory;
0804
0805 RefCountedLinearizedTrackState linState = linTrackFactory.linearizedTrackState(iter->position(), info.ghostTrack);
0806
0807 trackBundles[iter - vertices_.begin()].first = vertexTrackFactory.vertexTrack(linState, iter->vertexState());
0808 }
0809
0810 for (std::vector<RefCountedVertexTrack>::const_iterator track = vtxTracks.begin(); track != vtxTracks.end();
0811 ++track) {
0812 if (isGhostTrack(*track)) {
0813 trackBundles[iter - vertices_.begin()].first = *track;
0814 continue;
0815 }
0816
0817 if ((*track)->weight() < 1e-3) {
0818 trackBundles[iter - vertices_.begin()].second.push_back(*track);
0819 continue;
0820 }
0821
0822 unsigned int idx = iter - vertices_.begin();
0823 double best = 1.0e9;
0824 for (std::vector<CachingVertex<5> >::const_iterator vtx = vertices_.begin(); vtx != vertices_.end(); ++vtx) {
0825 if (info.pred.lambda(vtx->position()) < info.pred.lambda(vertices_[0].position()))
0826 continue;
0827
0828 double compat = sqr(trackVertexCompat(*vtx, *track));
0829
0830 compat /= (vtx == vertices_.begin()) ? primcut_ : seccut_;
0831
0832 if (compat < best) {
0833 best = compat;
0834 idx = vtx - vertices_.begin();
0835 }
0836 }
0837
0838 if ((int)idx != iter - vertices_.begin())
0839 assignmentChanged = true;
0840
0841 trackBundles[idx].second.push_back(*track);
0842 }
0843 }
0844
0845 if (!assignmentChanged)
0846 return false;
0847
0848 VertexTrackFactory<5> vertexTrackFactory;
0849 std::vector<CachingVertex<5> > vertices;
0850 vertices.reserve(vertices_.size());
0851
0852 for (std::vector<CachingVertex<5> >::const_iterator iter = vertices_.begin(); iter != vertices_.end(); ++iter) {
0853 const std::vector<RefCountedVertexTrack> &tracks = trackBundles[iter - vertices_.begin()].second;
0854 if (tracks.empty())
0855 continue;
0856
0857 CachingVertex<5> vtx;
0858
0859 if (tracks.size() == 1) {
0860 const TransientTrack &track = tracks[0]->linearizedTrack()->track();
0861
0862 int idx = -1;
0863 for (std::vector<GhostTrackState>::const_iterator iter = info.states.begin(); iter != info.states.end(); ++iter) {
0864 if (iter->isTrack() && iter->track() == track) {
0865 idx = iter - info.states.begin();
0866 break;
0867 }
0868 }
0869 if (idx < 0)
0870 continue;
0871
0872 vtx = vertexAtState(info.ghostTrack, info.pred, info.states[idx]);
0873 if (!vtx.isValid())
0874 continue;
0875 } else {
0876 std::vector<RefCountedVertexTrack> linTracks = relinearizeTracks(tracks, iter->vertexState());
0877
0878 if (fitType_ == kAlwaysWithGhostTrack)
0879 linTracks.push_back(
0880 relinearizeTrack(trackBundles[iter - vertices_.begin()].first, iter->vertexState(), vertexTrackFactory));
0881
0882 bool primary = iter == vertices_.begin();
0883 try {
0884 if (primary && info.hasBeamSpot)
0885 vtx = vertexFitter(true).vertex(linTracks, info.beamSpot.position(), info.beamSpot.error());
0886 else
0887 vtx = vertexFitter(primary).vertex(linTracks);
0888 } catch (const VertexException &e) {
0889
0890 }
0891 if (!vtx.isValid())
0892 return false;
0893 }
0894
0895 vertices.push_back(vtx);
0896 };
0897
0898 std::swap(vertices_, vertices);
0899 return true;
0900 }
0901
0902 void GhostTrackVertexFinder::refitGhostTrack(std::vector<CachingVertex<5> > &vertices, FinderInfo &info) const {
0903 VtxTrackIs isGhostTrack(info.ghostTrack);
0904
0905 std::vector<GhostTrackState> states;
0906 std::vector<unsigned int> oldStates;
0907 oldStates.reserve(info.states.size());
0908
0909 for (std::vector<CachingVertex<5> >::const_iterator iter = vertices.begin(); iter != vertices.end(); ++iter) {
0910 std::vector<RefCountedVertexTrack> vtxTracks = iter->tracks();
0911
0912 oldStates.clear();
0913 for (std::vector<RefCountedVertexTrack>::const_iterator track = vtxTracks.begin(); track != vtxTracks.end();
0914 ++track) {
0915 if (isGhostTrack(*track) || (*track)->weight() < 1e-3)
0916 continue;
0917
0918 const TransientTrack &tt = (*track)->linearizedTrack()->track();
0919
0920 int idx = -1;
0921 for (std::vector<GhostTrackState>::const_iterator iter = info.states.begin(); iter != info.states.end(); ++iter) {
0922 if (iter->isTrack() && iter->track() == tt) {
0923 idx = iter - info.states.begin();
0924 break;
0925 }
0926 }
0927
0928 if (idx >= 0)
0929 oldStates.push_back(idx);
0930 }
0931
0932 if (oldStates.size() == 1)
0933 states.push_back(info.states[oldStates[0]]);
0934 else
0935 states.push_back(GhostTrackState(iter->vertexState()));
0936 }
0937
0938 KalmanGhostTrackUpdater updater;
0939 SequentialGhostTrackFitter fitter;
0940 double ndof, chi2;
0941 info.pred = fitter.fit(updater, info.prior, states, ndof, chi2);
0942 TransientTrack ghostTrack = transientGhostTrack(info.pred, info.field);
0943
0944 std::swap(info.states, states);
0945 states.clear();
0946
0947 std::vector<CachingVertex<5> > newVertices;
0948 for (std::vector<CachingVertex<5> >::const_iterator iter = vertices.begin(); iter != vertices.end(); ++iter) {
0949 std::vector<RefCountedVertexTrack> vtxTracks = iter->tracks();
0950
0951 int idx = -1;
0952 bool redo = false;
0953 for (std::vector<RefCountedVertexTrack>::iterator track = vtxTracks.begin(); track != vtxTracks.end(); ++track) {
0954 if (isGhostTrack(*track)) {
0955 LinearizedTrackStateFactory linTrackFactory;
0956 VertexTrackFactory<5> vertexTrackFactory;
0957
0958 *track = vertexTrackFactory.vertexTrack(linTrackFactory.linearizedTrackState(iter->position(), ghostTrack),
0959 iter->vertexState());
0960 redo = true;
0961 continue;
0962 }
0963
0964 const TransientTrack &tt = (*track)->linearizedTrack()->track();
0965
0966 if (idx >= 0) {
0967 idx = -1;
0968 break;
0969 }
0970
0971 for (std::vector<GhostTrackState>::const_iterator it = info.states.begin(); it != info.states.end(); ++it) {
0972 if (!it->isTrack())
0973 continue;
0974
0975 if (it->track() == tt) {
0976 idx = it - info.states.begin();
0977 break;
0978 }
0979 }
0980 }
0981
0982 if (idx >= 0) {
0983 CachingVertex<5> vtx = vertexAtState(ghostTrack, info.pred, info.states[idx]);
0984 if (vtx.isValid())
0985 newVertices.push_back(vtx);
0986 } else if (redo) {
0987 bool primary = iter == vertices.begin();
0988 CachingVertex<5> vtx;
0989 if (primary && info.hasBeamSpot)
0990 vtx = vertexFitter(true).vertex(vtxTracks, info.beamSpot.position(), info.beamSpot.error());
0991 else
0992 vtx = vertexFitter(primary).vertex(vtxTracks);
0993 if (vtx.isValid())
0994 newVertices.push_back(vtx);
0995 } else
0996 newVertices.push_back(*iter);
0997 }
0998
0999 std::swap(newVertices, vertices);
1000 info.ghostTrack = ghostTrack;
1001 }
1002
1003
1004
1005 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const GhostTrack &ghostTrack,
1006 const CachingVertex<5> &primary,
1007 const BeamSpot &beamSpot,
1008 bool hasBeamSpot,
1009 bool hasPrimaries) const {
1010 FinderInfo info(primary, ghostTrack, beamSpot, hasBeamSpot, hasPrimaries);
1011
1012 std::vector<TransientVertex> result;
1013 if (info.states.empty())
1014 return result;
1015
1016 info.field = info.states[0].track().field();
1017 info.ghostTrack = transientGhostTrack(info.pred, info.field);
1018
1019 std::vector<CachingVertex<5> > vertices = initialVertices(info);
1020 if (primary.isValid()) {
1021 vertices.push_back(primary);
1022 if (vertices.size() > 1)
1023 std::swap(vertices.front(), vertices.back());
1024 }
1025
1026 unsigned int reassigned = 0;
1027 while (reassigned < 3) {
1028 if (vertices.size() < 2)
1029 break;
1030
1031 #ifdef DEBUG
1032 for (std::vector<CachingVertex<5> >::const_iterator iter = vertices.begin(); iter != vertices.end(); ++iter)
1033
1034 debugVertex(*iter, ghostTrack.prediction());
1035
1036 std::cout << "----- recursive merging: ---------" << std::endl;
1037 #endif
1038
1039 bool changed = recursiveMerge(vertices, info);
1040 if ((!changed && !reassigned) || vertices.size() < 2)
1041 break;
1042 if (changed)
1043 reassigned = 0;
1044
1045 if (fitType_ == kRefitGhostTrackWithVertices) {
1046 refitGhostTrack(vertices, info);
1047 changed = true;
1048 }
1049
1050 try {
1051 #ifdef DEBUG
1052 for (std::vector<CachingVertex<5> >::const_iterator iter = vertices.begin(); iter != vertices.end(); ++iter)
1053 debugVertex(*iter, ghostTrack.prediction());
1054 std::cout << "----- reassignment: ---------" << std::endl;
1055 #endif
1056 if (reassignTracks(vertices, info)) {
1057 reassigned++;
1058 changed = true;
1059 } else
1060 reassigned = 0;
1061 } catch (const VertexException &e) {
1062
1063 }
1064
1065 if (!changed)
1066 break;
1067 }
1068
1069 for (std::vector<CachingVertex<5> >::const_iterator iter = vertices.begin(); iter != vertices.end(); ++iter) {
1070 std::vector<RefCountedVertexTrack> tracks = iter->tracks();
1071 std::vector<RefCountedVertexTrack> newTracks;
1072 newTracks.reserve(tracks.size());
1073
1074 std::remove_copy_if(tracks.begin(), tracks.end(), std::back_inserter(newTracks), VtxTrackIs(info.ghostTrack));
1075
1076 if (newTracks.empty())
1077 continue;
1078
1079 CachingVertex<5> vtx(iter->vertexState(), newTracks, iter->totalChiSquared());
1080 result.push_back(vtx);
1081 }
1082
1083 return result;
1084 }