22 #include "CLHEP/Matrix/SymMatrix.h"
23 #include "CLHEP/Matrix/Vector.h"
57 const IInterface*
p) :
58 base_class(
t,
n,
p) {}
64 ATH_CHECK(m_ImpactPoint3dEstimator.retrieve());
68 return StatusCode::SUCCESS;
72 struct xAODVertex_pair {
80 std::pair<xAOD::VertexContainer*, xAOD::VertexAuxContainer*> InDetAdaptiveMultiSecVtxFinderTool::findVertex(
82 ATH_MSG_DEBUG(
" Number of input tracks before track selection: " << trackParticles->
size());
84 std::vector<Trk::ITrackLink*> selectedTracks;
93 vertexError.setZero();
94 null.setCovariancePosition(vertexError);
95 selectionPassed =
static_cast<bool>(m_trkFilter->accept(*itr, &
null));
96 if (selectionPassed) selectionPassed =
static_cast<bool>(m_SVtrkFilter->accept(*itr, &
null));
98 if (selectionPassed) {
99 trackPass(*itr) =
true;
102 par[0] = (itr)->hitPattern();
109 selectedTracks.push_back(linkTT);
113 ATH_MSG_DEBUG(
"Of " << trackParticles->size() <<
" tracks " << selectedTracks.size() <<
" survived the preselection.");
115 std::pair<xAOD::VertexContainer*, xAOD::VertexAuxContainer*> returnContainers = doVertexing(selectedTracks);
117 return returnContainers;
120 std::pair<xAOD::VertexContainer*, xAOD::VertexAuxContainer*> InDetAdaptiveMultiSecVtxFinderTool::doVertexing(
121 const std::vector<Trk::ITrackLink*>& trackVector) {
124 theVertexContainer->setStore(theVertexAuxContainer);
131 std::vector<xAODVertex_pair> myxAODVertices;
133 std::vector<Trk::TrackToVtxLink*> myTrackToVtxLinks;
136 std::vector<Trk::ITrackLink*> origTracks = trackVector;
137 std::vector<Trk::ITrackLink*> seedTracks = trackVector;
140 std::map<Trk::ITrackLink*, Trk::TrackToVtxLink*> TrackLinkOf;
148 TrackLinkOf[trkIter] = newTrkToVtxLink;
149 myTrackToVtxLinks.push_back(newTrkToVtxLink);
153 unsigned int seedtracknumber = seedTracks.size();
156 if (seedtracknumber == 0) {
ATH_MSG_DEBUG(
"New iteration. No tracks available after track selection for seeding."); }
159 ATH_MSG_DEBUG(
"Iteration number " << iteration <<
" and tracks left for seeding " << seedtracknumber);
161 std::vector<const Trk::TrackParameters*> perigeeList;
163 perigeeList.reserve(seedTracks.size());
168 std::unique_ptr<Trk::IMode3dInfo>
info;
170 seedVertex = m_SeedFinder->findSeed(m_privtx.x(), m_privtx.y(),
info, perigeeList);
172 ATH_MSG_DEBUG(
"Found seed at x: " << seedVertex.x() <<
" at y: " << seedVertex.y() <<
" at z: " << seedVertex.z());
175 theVertexContainer->
push_back(seededxAODVertex);
178 looseConstraintCovariance.setIdentity();
179 looseConstraintCovariance = looseConstraintCovariance * 1
e+8;
181 seededxAODVertex->
vxTrackAtVertex() = std::vector<Trk::VxTrackAtVertex>();
184 if (seedVertex.z() == 0) {
185 ATH_MSG_DEBUG(
"No good seed found: no further vertices in event");
186 ATH_MSG_DEBUG(
"Number of input tracks: " << perigeeList.size() <<
", but no good seed returned");
191 looseConstraintCovariance.setIdentity();
192 looseConstraintCovariance = looseConstraintCovariance * 1
e+8;
204 MvfFitInfo(*actualCandidate) =
206 isInitialized(*actualCandidate) =
false;
207 std::vector<Trk::VxTrackAtVertex*> vectorOfTracks(0);
208 VTAV(*actualCandidate) = vectorOfTracks;
212 float doe = findCompatibleTracks(seedVertex, trkIter);
213 if (doe < m_significanceCutSeeding) {
215 std::vector<xAOD::Vertex*>* actualvtxlink = actualLink->
vertices();
217 actualvtxlink->push_back(actualCandidate);
222 ATH_MSG_DEBUG(
" Considering n. " << VTAV(*actualCandidate).size() <<
" tracks for the fit. ");
223 if (actualCandidate !=
nullptr) {
224 if (VTAV(*actualCandidate).size() < 2) {
225 ATH_MSG_DEBUG(
"No tracks found near seed, while at least two tracks were expected.");
227 if (VTAV.isAvailable(*actualCandidate)) {
228 for (
auto *tav : VTAV(*actualCandidate)) {
229 if (tav ==
nullptr)
continue;
235 VTAV(*actualCandidate).clear();
237 if (MvfFitInfo.isAvailable(*actualCandidate) && MvfFitInfo(*actualCandidate) !=
nullptr) {
238 delete MvfFitInfo(*actualCandidate);
239 MvfFitInfo(*actualCandidate) =
nullptr;
241 delete actualCandidate;
242 actualCandidate =
nullptr;
249 m_VertexFitter->addVtxTofit(actualCandidate);
251 ATH_MSG_DEBUG(
"Deleting tracks with really good fit to vertex from seeding tracks.");
252 int nFound = removeTracksFromSeeds(actualCandidate, seedTracks);
254 ATH_MSG_DEBUG(
"Found and deleted " << nFound <<
" tracks from seeding tracks.");
256 ATH_MSG_DEBUG(
"All tracks used for fitting came from fiting tracks, removing closest from seeding.");
260 removeClosestTrack(seedVertex, seedTracks, nFound);
264 ATH_MSG_DEBUG(
"You still have not removed any tracks from seeds! Aborting.");
269 bool goodVertex = checkFit(actualCandidate);
272 ATH_MSG_DEBUG(
"Bad vertex, deleting the vertex and clearing all pointers");
276 if (actualCandidate) {
277 if (VTAV.isAvailable(*actualCandidate)) {
278 for (
auto *tav : VTAV(*actualCandidate)) {
279 if (tav ==
nullptr)
continue;
285 VTAV(*actualCandidate).clear();
287 if (MvfFitInfo.isAvailable(*actualCandidate) && MvfFitInfo(*actualCandidate) !=
nullptr) {
288 delete MvfFitInfo(*actualCandidate);
289 MvfFitInfo(*actualCandidate) =
nullptr;
292 delete actualCandidate;
293 actualCandidate =
nullptr;
301 myxAODVertices.emplace_back(0, actualCandidate);
303 seedtracknumber = seedTracks.size();
304 }
while (seedTracks.size() > 1 && iteration < m_maxIterations);
306 if (iteration >= m_maxIterations) {
308 << m_maxIterations <<
") reached; to reconstruct more vertices, set maxIterations to a higher value.");
311 ATH_MSG_DEBUG(
"Secondary vertex finding complete with " << iteration <<
" iterations and " << myxAODVertices.size()
312 <<
" vertices found.");
314 for (
const xAODVertex_pair& vtxIter : myxAODVertices) {
324 std::vector<Trk::VxTrackAtVertex>* tracksOfVertex = &(cand->
vxTrackAtVertex());
325 tracksOfVertex->clear();
328 if ((*MVFtrkIter).initialPerigee()) { (*MVFtrkIter).setPerigeeAtVertex(((*MVFtrkIter).initialPerigee())->clone()); }
329 tracksOfVertex->push_back(*MVFtrkIter);
333 for (
const xAODVertex_pair& vtxIter : myxAODVertices) {
339 MVFtrkIter =
nullptr;
342 delete MvfFitInfo(*cand);
348 std::vector<Trk::VxTrackAtVertex>* myVxTracksAtVtx = &((vxIter)->vxTrackAtVertex());
349 if (!myVxTracksAtVtx)
continue;
358 if (linkToXAODTP) { (vxIter)->addTrackAtVertex(*linkToXAODTP, (tracksIter).weight()); }
361 int ntrk = myVxTracksAtVtx->size();
365 bool isV0 = V0check(getVertexMomenta(vxIter), (&(*vxIter))->position());
376 if (!theVertexContainer->empty()) {
381 theVertexContainer->push_back(dummyxAODVertex);
384 dummyxAODVertex->
vxTrackAtVertex() = std::vector<Trk::VxTrackAtVertex>();
391 else if (theVertexContainer->empty()) {
393 theVertexContainer->push_back(dummyxAODVertex);
396 looseConstraintCovariance.setIdentity();
397 looseConstraintCovariance = looseConstraintCovariance * 1
e+8;
399 dummyxAODVertex->
vxTrackAtVertex() = std::vector<Trk::VxTrackAtVertex>();
408 for (
unsigned int i = 0;
i < theVertexContainer->size();
i++) {
420 ATH_MSG_DEBUG(
" Vtx: " <<
i <<
" x= " << (*theVertexContainer)[
i]->position().
x() <<
" y= "
421 << (*theVertexContainer)[
i]->position().
y() <<
" z= " << (*theVertexContainer)[
i]->position().
z()
422 <<
" ntracks= " << (*theVertexContainer)[
i]->vxTrackAtVertex().
size() <<
" chi2= "
423 << (*theVertexContainer)[
i]->
chiSquared() <<
" #dof = " << (*theVertexContainer)[
i]->numberDoF());
426 ATH_MSG_DEBUG(
"Done finding " << theVertexContainer->size() <<
" vertices and cleaning the container.");
427 ATH_MSG_DEBUG(
"Seeds good/bad/all : " << noVtx <<
"/" << kinkVtx <<
"/" << notSpec);
428 ATH_MSG_DEBUG(
"'Good' secondaries : " << secVtx <<
" and V0: " << V0vtx);
430 return std::make_pair(theVertexContainer, theVertexAuxContainer);
433 void InDetAdaptiveMultiSecVtxFinderTool::setPrimaryVertexPosition(
double vx,
double vy,
double vz) {
439 void InDetAdaptiveMultiSecVtxFinderTool::countTracksAndNdf(
xAOD::Vertex* myxAODVertex,
float&
ndf,
int& ntrk)
const {
449 if ((trkAtVtxIter)->
weight() > m_minWghtAtVtx) { ntrk += 1; }
454 const std::vector<Amg::Vector3D> InDetAdaptiveMultiSecVtxFinderTool::getVertexMomenta(
xAOD::Vertex* myxAODVertex)
const {
455 std::vector<Amg::Vector3D> TrkAtVtxMomenta;
457 std::vector<Trk::VxTrackAtVertex>* tracksAtVertex = &(myxAODVertex->
vxTrackAtVertex());
459 ATH_MSG_DEBUG(
" getVertexMomenta ... #Tracks associated at vertex : " << tracksAtVertex->size());
462 if ((tracksAtVertexIter).
weight() <= m_minWghtAtVtx)
continue;
470 double qp = 1. / (std::fabs(sv_perigee->parameters()[
Trk::qOverP]));
478 return TrkAtVtxMomenta;
481 bool InDetAdaptiveMultiSecVtxFinderTool::checkFit(
xAOD::Vertex* actualCandidate)
const {
488 if ((trkAtVtxIter)->
weight() > m_minWghtAtVtx) { ntracks += 1; }
491 ATH_MSG_DEBUG(
" xAOD::Vertex : " << (actualCandidate !=
nullptr ? 1 : 0) <<
", #dof = " <<
ndf
492 <<
", #tracks (weight>0.01) = " << ntracks);
494 return (actualCandidate !=
nullptr &&
ndf > 0 && ntracks >= 2);
497 int InDetAdaptiveMultiSecVtxFinderTool::removeTracksFromSeeds(
xAOD::Vertex* actualCandidate,
498 std::vector<Trk::ITrackLink*>& seedTracks)
const {
504 bool goodVertex = checkFit(actualCandidate);
516 if ((*seedtrkiter)->parameters() == (trkAtVtxIter)->trackOrParticleLink()->parameters() &&
517 (trkAtVtxIter)->
weight() > m_minWghtAtVtx) {
518 foundTrack = seedtrkiter;
522 if (foundTrack != seedEnd) {
523 seedTracks.erase(foundTrack);
527 seedBegin = seedTracks.begin();
528 seedEnd = seedTracks.end();
539 std::unique_ptr<Trk::PlaneSurface> mySurface =
540 m_ImpactPoint3dEstimator->Estimate3dIP((trkIter)->
parameters(), &seedVertex,
distance);
543 ATH_MSG_DEBUG(
" ImpactPoint3dEstimator failed to find minimum distance between track and vertex seed: " <<
err.p);
549 double doe = 99999999.9;
552 if (myPerigee && myPerigee->covariance()) {
568 void InDetAdaptiveMultiSecVtxFinderTool::removeClosestTrack(
Amg::Vector3D& seedVertex, std::vector<Trk::ITrackLink*>& seedTracks,
576 std::unique_ptr<Trk::PlaneSurface> mySurface =
577 m_ImpactPoint3dEstimator->Estimate3dIP((trkIter)->
parameters(), &seedVertex,
distance);
579 ATH_MSG_DEBUG(
" ImpactPoint3dEstimator failed to find minimum distance between this track and vertex seed: " <<
err.p);
584 if (
distance > 0 && !nearestTrack) {
586 nearestTrack = trkIter;
590 nearestTrack = trkIter;
599 if (foundTrack != seedEnd) {
600 seedTracks.erase(foundTrack);
601 seedBegin = seedTracks.begin();
602 seedEnd = seedTracks.end();
611 bool InDetAdaptiveMultiSecVtxFinderTool::V0check(
const std::vector<Amg::Vector3D>& momenta,
const Amg::Vector3D& posi)
const {
612 int ntrk = momenta.size();
619 std::vector<double> Pv(ntrk);
620 double vx = 0., vy = 0., vz = 0., eK0 = 0.;
621 double pi2 = 139.57018 * 139.57018;
623 for (
int t = 0;
t < ntrk;
t++) {
629 Pv[
t] = trk.x() * trk.x() + trk.y() * trk.y() + trk.z() * trk.z();
630 eK0 += std::sqrt(Pv[
t] +
pi2);
633 double mnt2 = vx * vx + vy * vy + vz * vz;
634 double mass = eK0 * eK0 - mnt2;
635 mass = 0.001 * (std::sqrt(std::abs(
mass)));
640 double modir = vmoment.dot(vdif) / std::sqrt(mnt2);
643 double a0z = (vdif + vmoment * vmoment.dot(vdif) / (mnt2 + 0.00001)).
z();
644 double Rxy = vdif.perp();
646 ATH_MSG_DEBUG(
" V0kine : a0z = " << a0z <<
" Rxy = " << Rxy <<
" direction " << modir);
653 if (a0z > 15. || Rxy > 500.) {
return false; }
658 double eGam = std::sqrt(Pv[0] + 0.511 * 0.511) + std::sqrt(Pv[1] + 0.511 * 0.511);
659 double mGam = eGam * eGam - mnt2;
661 double prtn2 = 938.27205 * 938.27205;
662 double eLam = Pv[0] > Pv[1] ? std::sqrt(Pv[0] + prtn2) + std::sqrt(Pv[1] +
pi2) : std::sqrt(Pv[0] +
pi2) + std::sqrt(Pv[1] + prtn2);
663 double mLam = eLam * eLam - mnt2;
665 ATH_MSG_DEBUG(
" V0 masses : " <<
mass <<
" " << std::sqrt(std::abs(mGam)) <<
" " << std::sqrt(std::abs(mLam)));
667 return ((fabs(
mass - 497.614) < 100.)
668 || (mGam > 0 && std::sqrt(mGam) < 40.)
669 || (mLam > 0 && std::abs(std::sqrt(mLam) - 1115.683) < 200.)