ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHActsToSvtxTracks.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHActsToSvtxTracks.cc
1 #include "PHActsToSvtxTracks.h"
3 
7 #include <phool/PHDataNode.h>
8 #include <phool/PHNode.h>
9 #include <phool/PHNodeIterator.h>
10 #include <phool/PHObject.h>
11 #include <phool/getClass.h>
12 #include <phool/phool.h>
13 #include <phool/PHTimer.h>
14 
16 #include <Acts/Utilities/Units.hpp>
17 
26 
28 #include <cmath>
29 #include <iostream>
30 #include <memory>
31 #include <utility>
32 
33 #include <TMatrixDSym.h>
34 
35 
37  : SubsysReco(name)
38  , m_actsFitResults(nullptr)
39 {
40  Verbosity(0);
41 }
42 
44 {
45 
46  if (Verbosity() > 10)
47  {
48  std::cout << "Finished PHActsToSvtxTracks" << std::endl;
49  }
50 
52 }
53 
55 {
57 }
58 
60 {
63 
66 
68 }
69 
71 {
72  if (Verbosity() > 1)
73  {
74  std::cout << "Start process_event in PHActsToSvtxTracks" << std::endl;
75  }
76 
77  for(auto& [trackKey, trajectory] : *m_actsFitResults)
78  {
79  createSvtxTrack(trackKey, trajectory);
80  }
81 
82  if(Verbosity() > 0)
83  std::cout << "Finished PHActsToSvtxTracks process_event"
84  << std::endl;
85 
87 }
88 
89 
91 {
93 }
94 
95 void PHActsToSvtxTracks::createSvtxTrack(const unsigned int trackKey,
96  Trajectory traj)
97 {
100  SvtxTrack *svtxTrack = m_svtxTrackMap->find(trackKey)->second;
101 
102  const auto vertexId = svtxTrack->get_vertex_id();
103 
104  auto svtxVertex = m_svtxVertexMap->find(vertexId)->second;
105 
106  Acts::Vector3D vertex(svtxVertex->get_x() * Acts::UnitConstants::cm,
107  svtxVertex->get_y() * Acts::UnitConstants::cm,
108  svtxVertex->get_z() * Acts::UnitConstants::cm);
109 
110  const auto &[trackTips, mj] = traj.trajectory();
111 
113  const auto& trackTip = trackTips.front();
114  const auto& params = traj.trackParameters(trackTip);
115 
116  svtxTrack->clear_states();
117  svtxTrack->clear_cluster_keys();
118 
119  float pathlength = 0.0;
120  SvtxTrackState_v1 out( pathlength);
121  out.set_x(0.0);
122  out.set_y(0.0);
123  out.set_z(0.0);
124  svtxTrack->insert_state(&out);
125 
126  auto trajState =
127  Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip);
128 
129  svtxTrack->set_x(params.position(m_tGeometry->geoContext)(0)
131  svtxTrack->set_y(params.position(m_tGeometry->geoContext)(1)
133  svtxTrack->set_z(params.position(m_tGeometry->geoContext)(2)
135 
136  svtxTrack->set_px(params.momentum()(0));
137  svtxTrack->set_py(params.momentum()(1));
138  svtxTrack->set_pz(params.momentum()(2));
139 
140  svtxTrack->set_charge(params.charge());
141  svtxTrack->set_chisq(trajState.chi2Sum);
142  svtxTrack->set_ndf(trajState.NDF);
143 
144  float dca3Dxy = NAN;
145  float dca3Dz = NAN;
146  float dca3DxyCov = NAN;
147  float dca3DzCov = NAN;
148 
149  auto rotater = std::make_unique<ActsTransformations>();
150  rotater->setVerbosity(Verbosity());
151 
152  if(params.covariance())
153  {
154  auto rotatedCov =
155  rotater->rotateActsCovToSvtxTrack(params,
157 
158  for(int i = 0; i < 6; i++)
159  for(int j = 0; j < 6; j++)
160  svtxTrack->set_error(i, j, rotatedCov(i,j));
161 
162  rotater->calculateDCA(params, vertex, rotatedCov,
164  dca3Dxy, dca3Dz,
165  dca3DxyCov, dca3DzCov);
166 
167  }
168 
169  svtxTrack->set_dca3d_xy(dca3Dxy / Acts::UnitConstants::cm);
170  svtxTrack->set_dca3d_z(dca3Dz / Acts::UnitConstants::cm);
171 
173  svtxTrack->set_dca3d_xy_error(sqrt(dca3DxyCov));
174  svtxTrack->set_dca3d_z_error(sqrt(dca3DzCov));
175 
176  rotater->fillSvtxTrackStates(traj, trackTip, svtxTrack,
178  return;
179 }
180 
182 {
183 
184  PHNodeIterator iter(topNode);
185 
186  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
187 
188  if (!dstNode)
189  {
190  std::cerr << "DST node is missing, quitting" << std::endl;
191  throw std::runtime_error("Failed to find DST node in PHActsTracks::createNodes");
192  }
193 
194  PHCompositeNode *svtxNode =
195  dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "SVTX"));
196 
197  if (!svtxNode)
198  {
199  svtxNode = new PHCompositeNode("SVTX");
200  dstNode->addNode(svtxNode);
201  }
202 
203  m_svtxTrackMap = findNode::getClass<SvtxTrackMap>(topNode,
204  m_svtxMapName.c_str());
205  if(!m_svtxTrackMap)
206  {
208  PHIODataNode<PHObject> *trackNode =
210  "PHObject");
211  dstNode->addNode(trackNode);
212  }
213 
215 }
216 
218 {
219 
220 
221  m_svtxTrackMap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
222 
223  if (!m_svtxTrackMap)
224  {
225  std::cout << PHWHERE << "SvtxTrackMap not found on node tree. Exiting."
226  << std::endl;
228  }
229 
230  m_svtxVertexMap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
231  if(!m_svtxVertexMap)
232  {
233  std::cout << PHWHERE << "SvtxVertexMap not found on node tree. Exiting."
234  << std::endl;
236  }
237 
238  m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode, "ActsTrackingGeometry");
239 
240  if(!m_tGeometry)
241  {
242  std::cout << PHWHERE << "ActsTrackingGeometry not on node tree, exiting."
243  << std::endl;
245  }
246 
247  m_actsFitResults = findNode::getClass<std::map<const unsigned int, Trajectory>>(topNode, "ActsFitResults");
248 
249  if(!m_actsFitResults)
250  {
251  std::cout << PHWHERE << "ActsFitResults not on node tree, exiting."
252  << std::endl;
254  }
255 
257 }
258