11 #if __cplusplus < 201402L
12 #include <boost/make_unique.hpp>
22 #include <ActsExamples/Plugins/BField/BFieldOptions.hpp>
23 #include <ActsExamples/Plugins/BField/ScalableBField.hpp>
54 , m_actsVertexMap(nullptr)
55 , m_trajectories(nullptr)
76 std::cout <<
"Starting event " <<
m_event <<
" in PHActsVertexFinder"
93 for(
auto track : trackPointers)
104 std::cout <<
"Finished PHActsVertexFinder::process_event" << std::endl;
122 std::cout <<
"Acts Final vertex finder succeeeded " <<
m_goodFits
123 <<
" out of " <<
m_totalFits <<
" events processed"
133 auto vertId =
track->get_vertex_id();
139 const auto trackZ =
track->get_z();
140 double closestVertZ = 9999;
145 double dz = fabs(trackZ - vertex->get_z());
146 if(dz < closestVertZ)
153 auto vertex = m_svtxVertexMap->get(vertId);
154 vertex->insert_track(key);
155 track->set_vertex_id(vertId);
163 std::vector<const Acts::BoundTrackParameters*> trackPtrs;
178 clusIter != svtxTrack->end_cluster_keys();
181 auto key = *clusIter;
191 if((nmaps + nintt + ntpc) < 35)
193 if(svtxTrack->get_pt() < 0.5)
198 const auto& [trackTips, mj] = traj.trajectory();
199 for(
const auto& trackTip : trackTips)
203 trackPtrs.push_back(params);
205 keyMap.insert(std::make_pair(params, key));
212 std::cout <<
"Finding vertices for the following number of tracks "
213 << trackPtrs.size() << std::endl;
215 for(
const auto& [param, key] : keyMap)
217 std::cout <<
"Track ID : " << key <<
" Track position: ("
235 return std::visit([tracks,
this](
auto &inputField) {
237 using InputMagneticField =
238 typename std::decay_t<decltype(inputField)>::element_type;
248 using ImpactPointEstimator =
255 static_assert(Acts::VertexFinderConcept<VertexSeeder>,
256 "VertexSeeder does not fulfill vertex finder concept.");
257 static_assert(Acts::VertexFinderConcept<VertexFinder>,
258 "VertexFinder does not fulfill vertex finder concept.");
265 MagneticField
bField(inputField);
266 auto propagator = std::make_shared<Propagator>(
Stepper(bField));
270 VertexFitter vertexFitter(std::move(vertexFitterConfig));
273 Linearizer linearizer(std::move(linearizerConfig));
276 ImpactPointEstimator ipEst(std::move(ipEstConfig));
279 VertexSeeder seeder(std::move(seederConfig));
282 std::move(linearizer),
283 std::move(seeder), ipEst);
285 finderConfig.reassignTracksAfterFirstFit =
true;
292 auto result = finder.find(tracks, finderOptions, state);
298 auto vertexCollection = *result;
303 std::cout <<
"Acts IVF found " << vertexCollection.size()
304 <<
" vertices in event" << std::endl;
307 for(
const auto& vertex : vertexCollection)
309 vertexVector.push_back(vertex);
316 std::cout <<
"Acts vertex finder returned error: "
317 << result.error().message() << std::endl;
334 unsigned int key = 0;
336 if(vertices.size() > 0)
339 else if (vertices.size() == 0)
341 auto svtxVertex = std::make_unique<SvtxVertex_v1>();
342 svtxVertex->set_x(0);
343 svtxVertex->set_y(0);
344 svtxVertex->set_z(0);
345 for(
int i = 0; i < 3; ++i)
347 for(
int j = 0; j < 3; ++j)
349 svtxVertex->set_error(i, j, 100.);
355 for(
auto& vertex : vertices)
357 const auto &[chi2, ndf] = vertex.fitQuality();
358 const auto numTracks = vertex.tracks().size();
362 std::cout <<
"Found vertex at (" << vertex.position().x()
363 <<
", " << vertex.position().y() <<
", "
364 << vertex.position().z() <<
")" << std::endl;
365 std::cout <<
"Vertex has ntracks = " << numTracks
366 <<
" with chi2/ndf " << chi2 / ndf << std::endl;
371 auto pair = std::make_pair(key, vertex);
375 #if __cplusplus < 201402L
376 auto svtxVertex = boost::make_unique<SvtxVertex_v1>();
378 auto svtxVertex = std::make_unique<SvtxVertex_v1>();
385 svtxVertex->set_x(vertexX);
386 svtxVertex->set_y(vertexY);
387 svtxVertex->set_z(vertexZ);
388 for(
int i = 0; i < 3; ++i)
390 for(
int j = 0; j < 3; ++j)
392 svtxVertex->set_error(i, j,
397 for(
const auto&
track : vertex.tracks())
399 const auto originalParams =
track.originalParams;
400 const auto trackKey = keyMap.find(originalParams)->second;
401 svtxVertex->insert_track(trackKey);
404 svtxTrack->set_vertex_id(key);
410 svtxVertex->set_chisq(chi2);
411 svtxVertex->set_ndof(ndf);
412 svtxVertex->set_t0(vertex.time());
413 svtxVertex->set_id(key);
422 std::cout <<
"Identify vertices in vertex map" << std::endl;
443 svtxTrack->get_pz());
448 for(
int i = 0; i < 3; ++i)
450 for(
int j = 0; j < 3; ++j)
452 posCov(i, j) = svtxTrack->get_error(i, j);
457 float phi = atan2(
r(1),
r(0));
462 rot(0,1) = -sin(phi);
471 rot_T = rot.transpose();
476 const auto dca3Dxy = pos_R(0);
477 const auto dca3Dz = pos_R(2);
478 const auto dca3DxyCov = rotCov(0,0);
479 const auto dca3DzCov = rotCov(2,2);
481 svtxTrack->set_dca3d_xy(dca3Dxy);
482 svtxTrack->set_dca3d_z(dca3Dz);
483 svtxTrack->set_dca3d_xy_error(sqrt(dca3DxyCov));
484 svtxTrack->set_dca3d_z_error(sqrt(dca3DzCov));
496 std::cerr <<
"DST node is missing, quitting" << std::endl;
497 throw std::runtime_error(
"Failed to find DST node in PHActsTracks::createNodes");
510 findNode::getClass<VertexMap>(topNode,
"ActsVertexMap");
523 findNode::getClass<SvtxVertexMap>(topNode,
"SvtxVertexMapActs");
529 "SvtxVertexMapActs",
"PHObject");
533 m_svtxVertexMap = findNode::getClass<SvtxVertexMap>(topNode,
"SvtxVertexMap");
548 m_trajectories = findNode::getClass<std::map<const unsigned int, Trajectory>>(topNode,
"ActsTrajectories");
551 std::cout <<
PHWHERE <<
"No trajectories on node tree, bailing."
556 m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,
557 "ActsTrackingGeometry");
560 std::cout <<
PHWHERE <<
"ActsTrackingGeometry not on node tree. Exiting"
569 std::cout <<
PHWHERE <<
"No SvtxTrackMap on node tree, exiting."