7 #include <Math/SVector.h>
8 #include <Math/SMatrix.h>
9 #include <Math/MatrixFunctions.h>
56 using namespace ROOT::Math;
58 typedef SVector<double, 3> Vector3;
60 typedef SMatrix<double, 3, 3, MatRepSym<double, 3> > Matrix3S;
61 typedef SMatrix<double, 2, 2, MatRepSym<double, 2> > Matrix2S;
62 typedef SMatrix<double, 2, 3> Matrix23;
65 RefCountedVertexTrack;
67 RefCountedLinearizedTrackState;
69 struct VtxTrackIs :
public std::unary_function<
70 RefCountedVertexTrack, bool> {
72 bool operator()(
const RefCountedVertexTrack &vtxTrack)
const
73 {
return vtxTrack->linearizedTrack()->
track() == track; }
87 static inline double sqr(
double arg) {
return arg *
arg; }
94 primaryVertex(primaryVertex),
pred(ghostTrack.prediction()),
96 beamSpot(beamSpot), hasBeamSpot(hasBeamSpot),
97 hasPrimaries(hasPrimaries),
field(0) {}
127 double weight = err1 / (err1 + err2);
129 return p1 + weight *
diff;
132 static VertexState
stateMean(
const VertexState &v1,
const VertexState &v2)
134 return VertexState(
vtxMean(v1.position(), v1.error(),
135 v2.position(), v2.error()),
136 v1.error() + v2.error());
140 const Matrix3S &
error,
double &chi2,
143 using namespace ROOT::Math;
152 Matrix2S measErr = Similarity(jacobian, error);
153 Matrix2S combErr = Similarity(jacobian, cov) + measErr;
154 if (!measErr.Invert() || !combErr.Invert())
157 cov -= SimilarityT(jacobian * cov, combErr);
158 chi2 += Similarity(jacobian * residual, measErr);
183 RefCountedLinearizedTrackState linState[2] = {
187 if ( !linState[0]->isValid() || !linState[1]->isValid() )
190 Matrix3S cov = SMatrixIdentity();
194 linState[0]->predictedStateParameters()[1],
195 linState[0]->predictedStateParameters()[2]) ||
197 linState[1]->predictedStateParameters()[1],
198 linState[1]->predictedStateParameters()[2]))
202 VertexState vtxState(point, error);
204 std::vector<RefCountedVertexTrack> linTracks(2);
205 linTracks[0] = vertexTrackFactory.
vertexTrack(linState[0], vtxState);
206 linTracks[1] = vertexTrackFactory.
vertexTrack(linState[1], vtxState);
212 const RefCountedVertexTrack &track,
213 const VertexState &state,
216 RefCountedLinearizedTrackState linTrack = track->linearizedTrack();
217 linTrack = linTrack->stateWithNewLinearizationPoint(state.position());
222 const std::vector<RefCountedVertexTrack> &
tracks,
223 const VertexState &state)
227 std::vector<RefCountedVertexTrack> finalTracks;
228 finalTracks.reserve(tracks.size());
230 for(std::vector<RefCountedVertexTrack>::const_iterator
iter =
231 tracks.begin();
iter != tracks.end(); ++
iter)
232 finalTracks.push_back(
241 double scale1,
double scale2)
247 return sqr(ROOT::Math::Dot(diff, diff)) / ROOT::Math::Similarity(cov, diff);
251 const RefCountedVertexTrack &vertexTrack)
253 using namespace ROOT::Math;
259 extrap.extrapolate(track.impactPointState(),
point);
266 Vector3
dir =
conv(point2 - point1);
268 tsos.cartesianError().matrix().Sub<Matrix3S>(0, 0);
272 return ROOT::Math::Similarity(
error,
dir);
277 mergeThreshold_(3.0),
280 fitType_(kAlwaysWithGhostTrack)
285 double maxFitChi2,
double mergeThreshold,
double primcut,
286 double seccut,
FitType fitType) :
287 maxFitChi2_(maxFitChi2),
288 mergeThreshold_(mergeThreshold),
309 std::auto_ptr<VertexFitter<5> > *ptr =
324 std::cout <<
"\tpt = " << mom.perp() <<
", eta = " << mom.eta() <<
", phi = " << mom.phi() <<
", charge = " << fts.
charge();
325 if (vtx && vtx->hasTrackWeight())
326 std::cout <<
", w = " << vtx->trackWeight(track) << std::endl;
334 for(std::vector<TransientTrack>::const_iterator
iter = tracks.begin();
336 debugTrack(*
iter, vtx);
348 static void debugTracks(
const std::vector<RefCountedVertexTrack> &tracks)
350 std::vector<TransientTrack>
tracks_;
351 for(std::vector<RefCountedVertexTrack>::const_iterator
iter = tracks.begin();
353 std::cout <<
"+ " << (*iter)->linearizedTrack()->linearizationPoint() << std::endl;
354 tracks_.push_back((*iter)->linearizedTrack()->track());
356 debugTracks(tracks_);
374 debugTracks(tracks, &vtx);
382 const std::vector<TransientTrack> &tracks)
const
386 direction, coneRadius, tracks);
394 const std::vector<TransientTrack> &tracks)
const
398 direction, coneRadius, beamSpot, tracks);
406 const std::vector<TransientTrack> &primaries,
407 const std::vector<TransientTrack> &tracks)
const
411 direction, coneRadius, beamSpot, primaries, tracks);
419 const std::vector<TransientTrack> &tracks)
const
423 direction, coneRadius, tracks);
426 std::vector<RefCountedVertexTrack>(), 0.);
428 std::vector<TransientVertex> result =
vertices(ghostTrack, primary);
431 for(std::vector<TransientVertex>::const_iterator
iter = result.begin();
445 const std::vector<TransientTrack> &tracks)
const
449 direction, coneRadius, tracks);
452 std::vector<RefCountedVertexTrack>(), 0.);
454 std::vector<TransientVertex> result =
455 vertices(ghostTrack, primary, beamSpot,
true);
458 for(std::vector<TransientVertex>::const_iterator
iter = result.begin();
472 const std::vector<TransientTrack> &primaries,
473 const std::vector<TransientTrack> &tracks)
const
477 direction, coneRadius, tracks);
479 std::vector<RefCountedVertexTrack> primaryVertexTracks;
480 if (!primaries.empty()) {
484 VertexState state(primaryPosition, primaryError);
486 for(std::vector<TransientTrack>::const_iterator
iter =
487 primaries.begin();
iter != primaries.end(); ++
iter) {
489 RefCountedLinearizedTrackState linState =
491 primaryPosition, *
iter);
493 primaryVertexTracks.push_back(
500 primaryVertexTracks, 0.);
502 std::vector<TransientVertex> result =
503 vertices(ghostTrack, primary, beamSpot,
true,
true);
506 for(std::vector<TransientVertex>::const_iterator
iter = result.begin();
520 std::vector<RefCountedVertexTrack>(), 0.);
522 std::vector<TransientVertex> result =
vertices(ghostTrack, primary);
525 for(std::vector<TransientVertex>::const_iterator
iter = result.begin();
541 std::vector<RefCountedVertexTrack>(), 0.);
543 std::vector<TransientVertex> result =
544 vertices(ghostTrack, primary, beamSpot,
true);
547 for(std::vector<TransientVertex>::const_iterator
iter = result.begin();
559 const std::vector<TransientTrack> &primaries,
562 std::vector<RefCountedVertexTrack> primaryVertexTracks;
563 if (!primaries.empty()) {
567 VertexState state(primaryPosition, primaryError);
569 for(std::vector<TransientTrack>::const_iterator
iter =
570 primaries.begin();
iter != primaries.end(); ++
iter) {
572 RefCountedLinearizedTrackState linState =
574 primaryPosition, *
iter);
576 primaryVertexTracks.push_back(
583 primaryVertexTracks, 0.);
585 std::vector<TransientVertex> result =
586 vertices(ghostTrack, primary, beamSpot,
true,
true);
589 for(std::vector<TransientVertex>::const_iterator
iter = result.begin();
613 beamSpot, ghostTrack);
619 const std::vector<TransientTrack> &primaries,
624 beamSpot, primaries, ghostTrack);
634 ghostTrack.
pz()), 0.05);
640 const std::vector<TransientTrack> &tracks,
641 const std::vector<float> &
weights)
const
643 GhostTrack ghostTrack_(ghostTrack, tracks, weights,
652 std::vector<RefCountedVertexTrack>(), 0.);
654 std::vector<TransientVertex> result =
vertices(ghostTrack_, primary);
657 for(std::vector<TransientVertex>::const_iterator
iter = result.begin();
669 const std::vector<TransientTrack> &tracks,
670 const std::vector<float> &
weights)
const
672 GhostTrack ghostTrack_(ghostTrack, tracks, weights,
680 std::vector<RefCountedVertexTrack>(), 0.);
682 std::vector<TransientVertex> result =
683 vertices(ghostTrack_, primary, beamSpot,
true);
686 for(std::vector<TransientVertex>::const_iterator
iter = result.begin();
698 const std::vector<TransientTrack> &primaries,
699 const std::vector<TransientTrack> &tracks,
700 const std::vector<float> &
weights)
const
702 GhostTrack ghostTrack_(ghostTrack, tracks, weights,
707 std::vector<RefCountedVertexTrack> primaryVertexTracks;
708 if (!primaries.empty()) {
716 for(std::vector<TransientTrack>::const_iterator
iter =
717 primaries.begin();
iter != primaries.end(); ++
iter) {
719 RefCountedLinearizedTrackState linState =
721 state.position(), *
iter);
723 primaryVertexTracks.push_back(
732 primaryVertexTracks, 0.);
734 std::vector<TransientVertex> result =
735 vertices(ghostTrack_, primary, beamSpot,
true,
true);
738 for(std::vector<TransientVertex>::const_iterator
iter = result.begin();
754 std::vector<CachingVertex<5> >
vertices;
755 for(std::vector<GhostTrackState>::const_iterator
iter =
758 if (!
iter->isValid())
771 vertices.push_back(vtx);
778 std::vector<RefCountedVertexTrack> &newTracks,
779 const VertexState &state,
780 const VtxTrackIs &ghostTrackFinder,
784 for(std::vector<RefCountedVertexTrack>::const_iterator
iter =
785 tracks.begin();
iter != tracks.end(); ++
iter) {
786 bool gt = ghostTrackFinder(*
iter);
787 if (gt && ghostTrack)
790 RefCountedVertexTrack track =
796 newTracks.push_back(*
iter);
804 bool isPrimary)
const
810 std::vector<RefCountedVertexTrack> linTracks;
812 RefCountedVertexTrack vtxGhostTrack;
815 isGhostTrack, vtxGhostTrack, vertexTrackFactory);
817 isGhostTrack, vtxGhostTrack, vertexTrackFactory);
823 vtxGhostTrack->linearizedTrack(),
824 vtxGhostTrack->vertexState()));
847 typedef std::pair<unsigned int, unsigned int>
Indices;
849 std::multimap<float, Indices> compatMap;
851 for(
unsigned int i = 0;
i <
n;
i++) {
853 for(
unsigned int j =
i + 1;
j <
n;
j++) {
867 bool changesMade =
false;
871 for(std::multimap<float, Indices>::const_iterator
iter =
872 compatMap.begin();
iter != compatMap.end(); ++
iter) {
873 unsigned int v1 =
iter->second.first;
874 unsigned int v2 =
iter->second.second;
887 std::multimap<float, Indices> newCompatMap;
889 if (
iter->second.first == v1 ||
890 iter->second.first == v2 ||
891 iter->second.second == v1 ||
892 iter->second.second == v2)
895 Indices indices =
iter->second;
896 indices.first -= indices.first > v2;
897 indices.second -= indices.second > v2;
899 newCompatMap.insert(std::make_pair(
900 iter->first, indices));
904 for(
unsigned int i = 0;
i <
n;
i++) {
937 std::vector<std::pair<RefCountedVertexTrack,
938 std::vector<RefCountedVertexTrack> > >
939 trackBundles(vertices_.size());
943 bool assignmentChanged =
false;
945 vertices_.begin();
iter != vertices_.end(); ++
iter) {
946 std::vector<RefCountedVertexTrack> vtxTracks =
iter->tracks();
948 if (vtxTracks.empty()) {
952 RefCountedLinearizedTrackState linState =
956 trackBundles[
iter - vertices_.begin()].first =
958 linState,
iter->vertexState());
961 for(std::vector<RefCountedVertexTrack>::const_iterator track =
962 vtxTracks.begin(); track != vtxTracks.end(); ++track) {
964 if (isGhostTrack(*track)) {
965 trackBundles[
iter - vertices_.begin()]
970 if ((*track)->weight() < 1
e-3) {
971 trackBundles[
iter - vertices_.begin()]
972 .second.push_back(*track);
976 unsigned int idx =
iter - vertices_.begin();
979 vtx = vertices_.begin();
980 vtx != vertices_.end(); ++vtx) {
988 compat /= (vtx == vertices_.begin()) ?
993 idx = vtx - vertices_.begin();
997 if ((
int)idx !=
iter - vertices_.begin())
998 assignmentChanged =
true;
1000 trackBundles[
idx].second.push_back(*track);
1004 if (!assignmentChanged)
1008 std::vector<CachingVertex<5> >
vertices;
1009 vertices.reserve(vertices_.size());
1012 vertices_.begin();
iter != vertices_.end(); ++
iter) {
1013 const std::vector<RefCountedVertexTrack> &tracks =
1014 trackBundles[
iter - vertices_.begin()].second;
1020 if (tracks.size() == 1) {
1022 tracks[0]->linearizedTrack()->track();
1025 for(std::vector<GhostTrackState>::const_iterator
iter =
1028 if (
iter->isTrack() &&
iter->track() == track) {
1041 std::vector<RefCountedVertexTrack> linTracks =
1046 trackBundles[
iter - vertices_.begin()].first,
1047 iter->vertexState(), vertexTrackFactory));
1049 bool primary =
iter == vertices_.begin();
1066 vertices.push_back(vtx);
1079 std::vector<GhostTrackState>
states;
1080 std::vector<unsigned int> oldStates;
1081 oldStates.reserve(info.
states.size());
1085 std::vector<RefCountedVertexTrack> vtxTracks =
iter->tracks();
1088 for(std::vector<RefCountedVertexTrack>::const_iterator track =
1089 vtxTracks.begin(); track != vtxTracks.end(); ++track) {
1091 if (isGhostTrack(*track) || (*track)->weight() < 1
e-3)
1095 (*track)->linearizedTrack()->track();
1098 for(std::vector<GhostTrackState>::const_iterator
iter =
1101 if (
iter->isTrack() &&
iter->track() ==
tt) {
1108 oldStates.push_back(idx);
1111 if (oldStates.size() == 1)
1112 states.push_back(info.
states[oldStates[0]]);
1120 info.
pred = fitter.
fit(updater, info.
prior, states, ndof, chi2);
1126 std::vector<CachingVertex<5> > newVertices;
1129 std::vector<RefCountedVertexTrack> vtxTracks =
iter->tracks();
1133 for(std::vector<RefCountedVertexTrack>::iterator track =
1134 vtxTracks.begin(); track != vtxTracks.end(); ++track) {
1136 if (isGhostTrack(*track)) {
1143 iter->vertexState());
1149 (*track)->linearizedTrack()->track();
1156 for(std::vector<GhostTrackState>::const_iterator it =
1158 it != info.
states.end(); ++it) {
1162 if (it->track() ==
tt) {
1163 idx = it - info.
states.begin();
1174 newVertices.push_back(vtx);
1186 newVertices.push_back(vtx);
1188 newVertices.push_back(*
iter);
1201 bool hasBeamSpot,
bool hasPrimaries)
const
1204 hasBeamSpot, hasPrimaries);
1206 std::vector<TransientVertex>
result;
1215 vertices.push_back(primary);
1216 if (vertices.size() > 1)
1217 std::swap(vertices.front(), vertices.back());
1220 unsigned int reassigned = 0;
1221 while(reassigned < 3) {
1222 if (vertices.size() < 2)
1227 vertices.begin();
iter != vertices.end(); ++
iter)
1231 std::cout <<
"----- recursive merging: ---------" << std::endl;
1235 if ((!changed && !reassigned) || vertices.size() < 2)
1248 iter = vertices.begin();
1251 std::cout <<
"----- reassignment: ---------" << std::endl;
1267 vertices.begin();
iter != vertices.end(); ++
iter) {
1268 std::vector<RefCountedVertexTrack> tracks =
iter->tracks();
1269 std::vector<RefCountedVertexTrack> newTracks;
1270 newTracks.reserve(tracks.size());
1272 std::remove_copy_if(tracks.begin(), tracks.end(),
1273 std::back_inserter(newTracks),
1276 if (newTracks.empty())
1280 iter->totalChiSquared());
1281 result.push_back(vtx);
GlobalError positionError() const
std::auto_ptr< VertexFitter< 5 > > primVertexFitter_
reco::Vertex::Point convertPos(const GlobalPoint &p)
~GhostTrackVertexFinder()
GhostTrackPrediction pred
const std::vector< reco::PFCandidatePtr > & tracks_
static bool covarianceUpdate(Matrix3S &cov, const Vector3 &residual, const Matrix3S &error, double &chi2, double theta, double phi)
double lambda(const GlobalPoint &point) const
const TransientTrack & track() const
static HepMC::IO_HEPEVT conv
void add(const reco::Track &track, double weight=1.0)
VertexState const & vertexState() const
Sin< T >::type sin(const T &t)
float totalChiSquared() const
static RefCountedVertexTrack relinearizeTrack(const RefCountedVertexTrack &track, const VertexState &state, const VertexTrackFactory< 5 > factory)
reco::Vertex::Error convertError(const GlobalError &ge)
Geom::Theta< T > theta() const
FreeTrajectoryState fts(const MagneticField *fieldProvider) const
static double vertexCompat(const CachingVertex< 5 > &vtx1, const CachingVertex< 5 > &vtx2, const FinderInfo &info, double scale1=1.0, double scale2=1.0)
bool reassignTracks(std::vector< CachingVertex< 5 > > &vertices, const FinderInfo &info) const
bool recursiveMerge(std::vector< CachingVertex< 5 > > &vertices, const FinderInfo &info) const
GlobalError cartesianError() const
TrackCharge charge() const
double px() const
x coordinate of momentum vector
const MagneticField * field
static TransientTrack transientGhostTrack(const GhostTrackPrediction &pred, const MagneticField *field)
static double fitChi2(const CachingVertex< 5 > &vtx)
std::auto_ptr< VertexFitter< 5 > > secVertexFitter_
const Point & position() const
position
static VertexState stateMean(const VertexState &v1, const VertexState &v2)
const AlgebraicSymMatrix33 & matrix_new() const
GhostTrack fit(const GlobalPoint &priorPosition, const GlobalError &priorError, const GlobalVector &direction, double coneRadius, const std::vector< TransientTrack > &tracks) const
const CachingVertex< 5 > & primaryVertex
std::vector< RefCountedVertexTrack > const & tracks() const
std::vector< TransientVertex > vertices(const reco::Vertex &primaryVertex, const GlobalVector &direction, double coneRadius, const std::vector< TransientTrack > &tracks) const
static std::vector< RefCountedVertexTrack > relinearizeTracks(const std::vector< RefCountedVertexTrack > &tracks, const VertexState &state)
const T & max(const T &a, const T &b)
virtual CachingVertex< N > vertex(const std::vector< reco::TransientTrack > &tracks) const =0
std::vector< reco::TransientTrack > const & originalTracks() const
float degreesOfFreedom() const
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
std::auto_ptr< GhostTrackFitter > ghostTrackFitter_
GlobalPoint position() const
Cos< T >::type cos(const T &t)
RefCountedLinearizedTrackState linearizedTrackState(const GlobalPoint &linP, const reco::TransientTrack &track) const
float totalChiSquared() const
static GhostTrackPrediction dummyPrediction(const Vertex &primaryVertex, const Track &ghostTrack)
static GlobalPoint vtxMean(const GlobalPoint &p1, const GlobalError &e1, const GlobalPoint &p2, const GlobalError &e2)
float degreesOfFreedom() const
GlobalVector momentum() const
double pz() const
z coordinate of momentum vector
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< float >, ROOT::Math::GlobalCoordinateSystemTag > GlobalVector
vector in glovbal coordinate system
FreeTrajectoryState initialFreeState() const
const Track & track() const
CachingVertex< 5 > mergeVertices(const CachingVertex< 5 > &vertex1, const CachingVertex< 5 > &vertex2, const FinderInfo &info, bool isPrimary) const
bool linearize(const GhostTrackPrediction &pred, bool initial=false, double lambda=0.)
tuple idx
DEBUGGING if hasattr(process,"trackMonIterativeTracking2012"): print "trackMonIterativeTracking2012 D...
reco::TransientTrack build(const FreeTrajectoryState &fts) const
std::vector< CachingVertex< 5 > > initialVertices(const FinderInfo &info) const
Error error() const
return SMatrix
FinderInfo(const CachingVertex< 5 > &primaryVertex, const GhostTrack &ghostTrack, const BeamSpot &beamSpot, bool hasBeamSpot, bool hasPrimaries)
VertexFitter< 5 > & vertexFitter(bool primary) const
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< float >, ROOT::Math::GlobalCoordinateSystemTag > GlobalPoint
point in global coordinate system
GlobalPoint position() const
GhostTrackPrediction fit(const GhostTrackFitter::PredictionUpdater &updater, const GhostTrackPrediction &prior, std::vector< GhostTrackState > &states, double &ndof, double &chi2)
static double trackVertexCompat(const CachingVertex< 5 > &vtx, const RefCountedVertexTrack &vertexTrack)
GhostTrackFitter & ghostTrackFitter() const
const GhostTrackPrediction & prediction() const
const math::XYZTLorentzVector & vectorSum() const
Square< F >::type sqr(const F &f)
GlobalPoint position(double lambda=0.) const
void refitGhostTrack(std::vector< CachingVertex< 5 > > &vertices, FinderInfo &info) const
std::vector< GhostTrackState > states
static double vtxErrorLong(const GlobalError &error, const GlobalVector &dir)
static void mergeTrackHelper(const std::vector< RefCountedVertexTrack > &tracks, std::vector< RefCountedVertexTrack > &newTracks, const VertexState &state, const VtxTrackIs &ghostTrackFinder, RefCountedVertexTrack &ghostTrack, const VertexTrackFactory< 5 > &factory)
TransientTrack ghostTrack
static CachingVertex< 5 > vertexAtState(const TransientTrack &ghostTrack, const GhostTrackPrediction &pred, const GhostTrackState &state)
GlobalPoint globalPosition() const
GlobalError error() const
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
GhostTrackPrediction prior
double py() const
y coordinate of momentum vector
Global3DVector GlobalVector
GlobalError positionError(double lambda=0.) const
const GlobalVector direction() const