ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHActsVertexFitter.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHActsVertexFitter.cc
1 #include "PHActsVertexFitter.h"
2 
5 #include <phool/PHIODataNode.h>
6 #include <phool/PHObject.h>
7 #include <phool/getClass.h>
8 #include <phool/phool.h>
9 
16 
22 
23 #include <ActsExamples/Plugins/BField/BFieldOptions.hpp>
24 #include <ActsExamples/Plugins/BField/ScalableBField.hpp>
25 
39 
40 #include <iostream>
41 
43  : SubsysReco(name)
44  , m_actsFitResults(nullptr)
45  , m_tGeometry(nullptr)
46 {
47 }
48 
50 {
51  if (Verbosity() > 1)
52  std::cout << "PHActsVertexFitter::Init" << std::endl;
53 
55 }
56 
58 {
59  if (Verbosity() > 1)
60  std::cout << "PHActsVertexFitter::End " << std::endl;
62 }
63 
65 {
67 
69 }
70 
72 {
73  if (getNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
75 
78 
80 }
81 
83 {
84  auto logLevel = Acts::Logging::FATAL;
85  if (Verbosity() > 0)
86  {
87  std::cout << "Beginning PHActsVertexFitter::process_event number "
88  << m_event << std::endl;
89  if (Verbosity() > 5)
90  logLevel = Acts::Logging::VERBOSE;
91  }
92 
93  const auto vertexTrackMap = getTracks();
94 
95  for (const auto &[vertexId, trackVec] : vertexTrackMap)
96  {
97  const auto vertex = fitVertex(trackVec, logLevel);
98 
99  createActsSvtxVertex(vertexId, vertex);
100 
102  updateSvtxVertex(vertexId, vertex);
103  }
104 
105  if (Verbosity() > 1)
106  std::cout << "Finished PHActsVertexFitter::process_event"
107  << std::endl;
108 
110 }
111 
112 void PHActsVertexFitter::updateSvtxVertex(const unsigned int vertexId,
113  ActsVertex vertex)
114 {
115  auto svtxVertex = m_vertexMap->get(vertexId);
116 
117  if (Verbosity() > 1)
118  std::cout << "Updating SvtxVertex id " << vertexId << std::endl;
119 
120  svtxVertex->set_x(vertex.position().x() / Acts::UnitConstants::cm);
121  svtxVertex->set_y(vertex.position().y() / Acts::UnitConstants::cm);
122  svtxVertex->set_z(vertex.position().z() / Acts::UnitConstants::cm);
123 
124  for (int i = 0; i < 3; ++i)
125  {
126  for (int j = 0; j < 3; ++j)
127  {
128  svtxVertex->set_error(i, j,
129  vertex.covariance()(i, j) / Acts::UnitConstants::cm2);
130  }
131  }
132 
133  const auto &[chi2, ndf] = vertex.fitQuality();
134  svtxVertex->set_ndof(ndf);
135  svtxVertex->set_chisq(chi2);
136  svtxVertex->set_t0(vertex.time());
137 }
138 
139 void PHActsVertexFitter::createActsSvtxVertex(const unsigned int vertexId,
140  ActsVertex vertex)
141 {
142 #if __cplusplus < 201402L
143  auto svtxVertex = boost::make_unique<SvtxVertex_v1>();
144 #else
145  auto svtxVertex = std::make_unique<SvtxVertex_v1>();
146 #endif
147 
148  if (Verbosity() > 2)
149  {
150  std::cout << "Creating vertex for id " << vertexId
151  << std::endl;
152  }
153 
154  svtxVertex->set_x(vertex.position().x() / Acts::UnitConstants::cm);
155  svtxVertex->set_y(vertex.position().y() / Acts::UnitConstants::cm);
156  svtxVertex->set_z(vertex.position().z() / Acts::UnitConstants::cm);
157 
158  for (int i = 0; i < 3; ++i)
159  {
160  for (int j = 0; j < 3; ++j)
161  {
162  svtxVertex->set_error(i, j,
163  vertex.covariance()(i, j) / Acts::UnitConstants::cm2);
164  }
165  }
166 
167  const auto &[chi2, ndf] = vertex.fitQuality();
168  svtxVertex->set_ndof(ndf);
169  svtxVertex->set_chisq(chi2);
170  svtxVertex->set_t0(vertex.time());
171  svtxVertex->set_id(vertexId);
172 
173  m_actsVertexMap->insert(svtxVertex.release());
174 }
175 
177 {
180  return std::visit([tracks, logLevel, this](auto &inputField) {
182  using InputMagneticField =
183  typename std::decay_t<decltype(inputField)>::element_type;
184  using MagneticField = Acts::SharedBField<InputMagneticField>;
187  using PropagatorOptions = Acts::PropagatorOptions<>;
190  using VertexFitter =
192  using VertexFitterOptions = Acts::VertexingOptions<TrackParameters>;
193 
194  auto logger = Acts::getDefaultLogger("PHActsVertexFitter",
195  logLevel);
196 
198  MagneticField bField(inputField);
199  auto propagator = std::make_shared<Propagator>(Stepper(bField));
200  PropagatorOptions propagatorOpts(m_tGeometry->geoContext,
202  Acts::LoggerWrapper(*logger));
203 
204  typename VertexFitter::Config vertexFitterCfg;
205  VertexFitter fitter(vertexFitterCfg);
207 
208  typename Linearizer::Config linConfig(bField, propagator);
209  Linearizer linearizer(linConfig);
210 
212  VertexFitterOptions vfOptions(m_tGeometry->geoContext,
214 
216  auto fitRes = fitter.fit(tracks, linearizer,
217  vfOptions, state);
218 
219  Acts::Vertex<TrackParameters> fittedVertex;
220 
221  if (fitRes.ok())
222  {
223  fittedVertex = *fitRes;
224  if (Verbosity() > 3)
225  {
226  std::cout << "Fitted vertex position "
227  << fittedVertex.position().x()
228  << ", "
229  << fittedVertex.position().y()
230  << ", "
231  << fittedVertex.position().z()
232  << std::endl;
233  }
234  }
235  else
236  {
237  if (Verbosity() > 3)
238  {
239  std::cout << "Acts vertex fit error: "
240  << fitRes.error().message()
241  << std::endl;
242  }
243  }
244 
245  return fittedVertex;
246  },
248 }
249 
251 {
252  VertexTrackMap trackPtrs;
253 
254  for (const auto &[key, track] : *m_trackMap)
255  {
256  const unsigned int vertexId = track->get_vertex_id();
257  const auto trackParam = makeTrackParam(track);
258  auto trackVecPos = trackPtrs.find(vertexId);
259 
260  if (trackVecPos == trackPtrs.end())
261  {
262  BoundTrackParamVec trackVec;
263 
264  trackVec.push_back(trackParam);
265  auto pair = std::make_pair(vertexId, trackVec);
266  trackPtrs.insert(pair);
267  }
268  else
269  {
270  trackVecPos->second.push_back(trackParam);
271  }
272  }
273 
274  if (Verbosity() > 3)
275  {
276  for (const auto &[vertexId, trackVec] : trackPtrs)
277  {
278  std::cout << "Fitting vertexId : " << vertexId
279  << " with the following number of tracks "
280  << trackVec.size()
281  << std::endl;
282 
283  for (const auto param : trackVec)
284  {
285  std::cout << "Track position: ("
286  << param->position(m_tGeometry->geoContext)(0)
287  << ", " << param->position(m_tGeometry->geoContext)(1) << ", "
288  << param->position(m_tGeometry->geoContext)(2) << ")"
289  << std::endl;
290  }
291  }
292  }
293 
294  return trackPtrs;
295 }
296 
298 {
299  const Acts::Vector4D trackPos(
300  track->get_x() * Acts::UnitConstants::cm,
301  track->get_y() * Acts::UnitConstants::cm,
302  track->get_z() * Acts::UnitConstants::cm,
304 
305  const Acts::Vector3D trackMom(track->get_px(),
306  track->get_py(),
307  track->get_pz());
308 
309  const int trackQ = track->get_charge() * Acts::UnitConstants::e;
310  const double p = track->get_p();
312 
313  for (int i = 0; i < 6; ++i)
314  {
315  for (int j = 0; j < 6; ++j)
316  {
317  cov(i, j) = track->get_error(i, j);
318  }
319  }
320 
321  auto perigee = Acts::Surface::makeShared<Acts::PerigeeSurface>(
323  track->get_y() * Acts::UnitConstants::cm,
324  track->get_z() * Acts::UnitConstants::cm));
325 
326  const auto param = new Acts::BoundTrackParameters(
327  perigee, m_tGeometry->geoContext,
328  trackPos, trackMom, p, trackQ, cov);
329 
330  return param;
331 }
332 
334 {
335  m_actsFitResults = findNode::getClass<std::map<const unsigned int, Trajectory>>(topNode, "ActsFitResults");
336  if (!m_actsFitResults)
337  {
338  std::cout << PHWHERE << "Acts Trajectories not found on node tree, exiting."
339  << std::endl;
341  }
342 
343  m_trackMap = findNode::getClass<SvtxTrackMap>(topNode,
344  "SvtxTrackMap");
345  if (!m_trackMap)
346  {
347  std::cout << PHWHERE << "No SvtxTrackMap on node tree. Bailing."
348  << std::endl;
350  }
351 
352  m_vertexMap = findNode::getClass<SvtxVertexMap>(topNode,
353  "SvtxVertexMap");
354  if (!m_vertexMap)
355  {
356  std::cout << PHWHERE << "No SvtxVertexMap on node tree, bailing"
357  << std::endl;
359  }
360 
361  m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode, "ActsTrackingGeometry");
362  if (!m_tGeometry)
363  {
364  std::cout << PHWHERE << "ActsTrackingGeometry not on node tree. Exiting"
365  << std::endl;
367  }
368 
370 }
371 
373 {
374  PHNodeIterator iter(topNode);
375 
376  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
377 
378  if (!dstNode)
379  {
380  std::cerr << "DST node is missing, quitting" << std::endl;
381  throw std::runtime_error("Failed to find DST node in PHActsTracks::createNodes");
382  }
383 
384  PHCompositeNode *svtxNode =
385  dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "SVTX"));
386 
387  if (!svtxNode)
388  {
389  svtxNode = new PHCompositeNode("SVTX");
390  dstNode->addNode(svtxNode);
391  }
392 
394  findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMapActs");
395  if (!m_actsVertexMap)
396  {
398  PHIODataNode<PHObject> *node =
400  "SvtxVertexMapActs", "PHObject");
401  svtxNode->addNode(node);
402  }
403 
405 }