ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHActsVertexFinder.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHActsVertexFinder.cc
1 #include "PHActsVertexFinder.h"
3 
6 #include <phool/PHIODataNode.h>
7 #include <phool/PHObject.h>
8 #include <phool/getClass.h>
9 #include <phool/phool.h>
10 
11 #if __cplusplus < 201402L
12 #include <boost/make_unique.hpp>
13 #endif
14 
21 
22 #include <ActsExamples/Plugins/BField/BFieldOptions.hpp>
23 #include <ActsExamples/Plugins/BField/ScalableBField.hpp>
24 
35 #include <Acts/Utilities/Units.hpp>
45 
48 
49 #include <memory>
50 #include <iostream>
51 
53  : PHInitVertexing(name)
54  , m_actsVertexMap(nullptr)
55  , m_trajectories(nullptr)
56 {
57 }
58 
60 {
61  int ret = getNodes(topNode);
63  return ret;
64 
65  ret = createNodes(topNode);
67  return ret;
68 
70 }
71 
73 {
74  if(Verbosity() > 0)
75  {
76  std::cout << "Starting event " << m_event << " in PHActsVertexFinder"
77  << std::endl;
78  }
79 
81  KeyMap keyMap;
82 
84  auto trackPointers = getTracks(keyMap);
85 
86  auto vertices = findVertices(trackPointers);
87 
88  fillVertexMap(vertices, keyMap);
89 
91 
93  for(auto track : trackPointers)
94  {
95  delete track;
96  }
97 
98  for(auto [key, svtxVertex] : *m_svtxVertexMap)
99  {
100  m_svtxVertexMapActs->insert(dynamic_cast<SvtxVertex*>(svtxVertex->CloneMe()));
101  }
102 
103  if(Verbosity() > 0)
104  std::cout << "Finished PHActsVertexFinder::process_event" << std::endl;
105 
106  m_event++;
107 
109 }
110 
112 {
113  m_actsVertexMap->clear();
115 
117 
118 }
119 
121 {
122  std::cout << "Acts Final vertex finder succeeeded " << m_goodFits
123  << " out of " << m_totalFits << " events processed"
124  << std::endl;
125 
127 }
128 
130 {
131  for(auto& [key, track] : *m_svtxTrackMap)
132  {
133  auto vertId = track->get_vertex_id();
134  if(!m_svtxVertexMap->get(vertId))
135  {
138 
139  const auto trackZ = track->get_z();
140  double closestVertZ = 9999;
141  vertId = UINT_MAX;
142 
143  for(auto& [vertexKey, vertex] : *m_svtxVertexMap)
144  {
145  double dz = fabs(trackZ - vertex->get_z());
146  if(dz < closestVertZ)
147  {
148  vertId = vertexKey;
149  closestVertZ = dz;
150  }
151  }
152 
153  auto vertex = m_svtxVertexMap->get(vertId);
154  vertex->insert_track(key);
155  track->set_vertex_id(vertId);
156  }
157  }
158 
159 }
160 
162 {
163  std::vector<const Acts::BoundTrackParameters*> trackPtrs;
164 
165  for(const auto& [key, traj] : *m_trajectories)
166  {
167  const auto svtxTrack = m_svtxTrackMap->get(key);
168 
169  // Track was removed by the track cleaner, skip it
170  if(!svtxTrack)
171  { continue; }
172 
173  int nmaps = 0;
174  int nintt = 0;
175  int ntpc = 0;
176 
177  for(SvtxTrack::ConstClusterKeyIter clusIter = svtxTrack->begin_cluster_keys();
178  clusIter != svtxTrack->end_cluster_keys();
179  ++clusIter)
180  {
181  auto key = *clusIter;
182  auto trkrId = TrkrDefs::getTrkrId(key);
183  if(trkrId == TrkrDefs::TrkrId::mvtxId)
184  { nmaps++; }
185  if(trkrId == TrkrDefs::TrkrId::inttId)
186  { nintt++; }
187  if(trkrId == TrkrDefs::TrkrId::tpcId)
188  { ntpc++; }
189  }
190 
191  if((nmaps + nintt + ntpc) < 35)
192  { continue; }
193  if(svtxTrack->get_pt() < 0.5)
194  { continue; }
195  if(nmaps < 1)
196  { continue; }
197 
198  const auto& [trackTips, mj] = traj.trajectory();
199  for(const auto& trackTip : trackTips)
200  {
201  const Acts::BoundTrackParameters* params =
202  new Acts::BoundTrackParameters(traj.trackParameters(trackTip));
203  trackPtrs.push_back(params);
204 
205  keyMap.insert(std::make_pair(params, key));
206 
207  }
208  }
209 
210  if(Verbosity() > 3)
211  {
212  std::cout << "Finding vertices for the following number of tracks "
213  << trackPtrs.size() << std::endl;
214 
215  for(const auto& [param, key] : keyMap)
216  {
217  std::cout << "Track ID : " << key << " Track position: ("
218  << param->position(m_tGeometry->geoContext).transpose()
219  << std::endl;
220  }
221  }
222 
223  return trackPtrs;
224 
225 }
226 
228 {
229  m_totalFits++;
230 
231  auto field = m_tGeometry->magField;
232 
235  return std::visit([tracks, this](auto &inputField) {
237  using InputMagneticField =
238  typename std::decay_t<decltype(inputField)>::element_type;
239  using MagneticField = Acts::SharedBField<InputMagneticField>;
240 
243  using PropagatorOptions = Acts::PropagatorOptions<>;
246  using VertexFitter =
248  using ImpactPointEstimator =
250  using VertexSeeder = Acts::ZScanVertexFinder<VertexFitter>;
251  using VertexFinder =
253  using VertexFinderOptions = Acts::VertexingOptions<TrackParameters>;
254 
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.");
259 
260  auto logLevel = Acts::Logging::FATAL;
261  if(Verbosity() > 4)
262  logLevel = Acts::Logging::VERBOSE;
263  auto logger = Acts::getDefaultLogger("PHActsVertexFinder", logLevel);
264 
265  MagneticField bField(inputField);
266  auto propagator = std::make_shared<Propagator>(Stepper(bField));
267 
269  typename VertexFitter::Config vertexFitterConfig;
270  VertexFitter vertexFitter(std::move(vertexFitterConfig));
271 
272  typename Linearizer::Config linearizerConfig(bField, propagator);
273  Linearizer linearizer(std::move(linearizerConfig));
274 
275  typename ImpactPointEstimator::Config ipEstConfig(bField, propagator);
276  ImpactPointEstimator ipEst(std::move(ipEstConfig));
277 
278  typename VertexSeeder::Config seederConfig(ipEst);
279  VertexSeeder seeder(std::move(seederConfig));
280 
281  typename VertexFinder::Config finderConfig(std::move(vertexFitter),
282  std::move(linearizer),
283  std::move(seeder), ipEst);
284  finderConfig.maxVertices = m_maxVertices;
285  finderConfig.reassignTracksAfterFirstFit = true;
286  VertexFinder finder(finderConfig, std::move(logger));
287 
289  VertexFinderOptions finderOptions(m_tGeometry->geoContext,
291 
292  auto result = finder.find(tracks, finderOptions, state);
293 
294  VertexVector vertexVector;
295 
296  if(result.ok())
297  {
298  auto vertexCollection = *result;
299  m_goodFits++;
300 
301  if(Verbosity() > 1)
302  {
303  std::cout << "Acts IVF found " << vertexCollection.size()
304  << " vertices in event" << std::endl;
305  }
306 
307  for(const auto& vertex : vertexCollection)
308  {
309  vertexVector.push_back(vertex);
310  }
311  }
312  else
313  {
314  if(Verbosity() > 1)
315  {
316  std::cout << "Acts vertex finder returned error: "
317  << result.error().message() << std::endl;
318  }
319  }
320 
321  return vertexVector;
322 
323  }
324  , field
325  );
326 }
327 
328 
329 
330 
332  KeyMap& keyMap)
333 {
334  unsigned int key = 0;
335 
336  if(vertices.size() > 0)
338 
339  else if (vertices.size() == 0)
340  {
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)
346  {
347  for(int j = 0; j < 3; ++j)
348  {
349  svtxVertex->set_error(i, j, 100.);
350  }
351  }
352  m_svtxVertexMap->insert(svtxVertex.release());
353  }
354 
355  for(auto& vertex : vertices)
356  {
357  const auto &[chi2, ndf] = vertex.fitQuality();
358  const auto numTracks = vertex.tracks().size();
359 
360  if(Verbosity() > 1)
361  {
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;
367  }
368 
369 
371  auto pair = std::make_pair(key, vertex);
372  m_actsVertexMap->insert(pair);
373 
375  #if __cplusplus < 201402L
376  auto svtxVertex = boost::make_unique<SvtxVertex_v1>();
377  #else
378  auto svtxVertex = std::make_unique<SvtxVertex_v1>();
379  #endif
380 
381  const auto vertexX = vertex.position().x() / Acts::UnitConstants::cm;
382  const auto vertexY = vertex.position().y() / Acts::UnitConstants::cm;
383  const auto vertexZ = vertex.position().z() / Acts::UnitConstants::cm;
384 
385  svtxVertex->set_x(vertexX);
386  svtxVertex->set_y(vertexY);
387  svtxVertex->set_z(vertexZ);
388  for(int i = 0; i < 3; ++i)
389  {
390  for(int j = 0; j < 3; ++j)
391  {
392  svtxVertex->set_error(i, j,
393  vertex.covariance()(i,j) / Acts::UnitConstants::cm2);
394  }
395  }
396 
397  for(const auto& track : vertex.tracks())
398  {
399  const auto originalParams = track.originalParams;
400  const auto trackKey = keyMap.find(originalParams)->second;
401  svtxVertex->insert_track(trackKey);
402 
403  const auto svtxTrack = m_svtxTrackMap->find(trackKey)->second;
404  svtxTrack->set_vertex_id(key);
405  updateTrackDCA(trackKey, Acts::Vector3D(vertexX,
406  vertexY,
407  vertexZ));
408  }
409 
410  svtxVertex->set_chisq(chi2);
411  svtxVertex->set_ndof(ndf);
412  svtxVertex->set_t0(vertex.time());
413  svtxVertex->set_id(key);
414 
415  m_svtxVertexMap->insert(svtxVertex.release());
416 
417  ++key;
418  }
419 
420  if(Verbosity() > 2)
421  {
422  std::cout << "Identify vertices in vertex map" << std::endl;
423  for(auto& [key, vert] : *m_svtxVertexMap)
424  {
425  vert->identify();
426  }
427  }
428 
429  return;
430 }
431 
432 void PHActsVertexFinder::updateTrackDCA(const unsigned int trackKey,
433  const Acts::Vector3D vertex)
434 {
435 
436  auto svtxTrack = m_svtxTrackMap->find(trackKey)->second;
437 
438  Acts::Vector3D pos(svtxTrack->get_x(),
439  svtxTrack->get_y(),
440  svtxTrack->get_z());
441  Acts::Vector3D mom(svtxTrack->get_px(),
442  svtxTrack->get_py(),
443  svtxTrack->get_pz());
444 
445  pos -= vertex;
446 
448  for(int i = 0; i < 3; ++i)
449  {
450  for(int j = 0; j < 3; ++j)
451  {
452  posCov(i, j) = svtxTrack->get_error(i, j);
453  }
454  }
455 
456  Acts::Vector3D r = mom.cross(Acts::Vector3D(0.,0.,1.));
457  float phi = atan2(r(1), r(0));
458 
461  rot(0,0) = cos(phi);
462  rot(0,1) = -sin(phi);
463  rot(0,2) = 0;
464  rot(1,0) = sin(phi);
465  rot(1,1) = cos(phi);
466  rot(1,2) = 0;
467  rot(2,0) = 0;
468  rot(2,1) = 0;
469  rot(2,2) = 1;
470 
471  rot_T = rot.transpose();
472 
473  Acts::Vector3D pos_R = rot * pos;
474  Acts::ActsSymMatrixD<3> rotCov = rot * posCov * rot_T;
475 
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);
480 
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));
485 
486 }
487 
489 {
490  PHNodeIterator iter(topNode);
491 
492  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
493 
494  if (!dstNode)
495  {
496  std::cerr << "DST node is missing, quitting" << std::endl;
497  throw std::runtime_error("Failed to find DST node in PHActsTracks::createNodes");
498  }
499 
500  PHCompositeNode *svtxNode =
501  dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "SVTX"));
502 
503  if (!svtxNode)
504  {
505  svtxNode = new PHCompositeNode("SVTX");
506  dstNode->addNode(svtxNode);
507  }
508 
509  m_actsVertexMap =
510  findNode::getClass<VertexMap>(topNode, "ActsVertexMap");
511  if(!m_actsVertexMap)
512  {
514 
515  PHDataNode<VertexMap> *node =
517  "ActsVertexMap");
518 
519  svtxNode->addNode(node);
520  }
521 
523  findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMapActs");
525  {
527  PHIODataNode<PHObject> *node =
529  "SvtxVertexMapActs", "PHObject");
530  svtxNode->addNode(node);
531  }
532 
533  m_svtxVertexMap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
534  if(!m_svtxVertexMap)
535  {
538  "SvtxVertexMap",
539  "PHObject");
540  svtxNode->addNode(node);
541  }
542 
544 }
545 
547 {
548  m_trajectories = findNode::getClass<std::map<const unsigned int, Trajectory>>(topNode, "ActsTrajectories");
549  if(!m_trajectories)
550  {
551  std::cout << PHWHERE << "No trajectories on node tree, bailing."
552  << std::endl;
554  }
555 
556  m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,
557  "ActsTrackingGeometry");
558  if(!m_tGeometry)
559  {
560  std::cout << PHWHERE << "ActsTrackingGeometry not on node tree. Exiting"
561  << std::endl;
563  }
564 
565  m_svtxTrackMap = findNode::getClass<SvtxTrackMap>(topNode,
566  "SvtxTrackMap");
567  if(!m_svtxTrackMap)
568  {
569  std::cout << PHWHERE << "No SvtxTrackMap on node tree, exiting."
570  << std::endl;
572  }
573 
575 }