ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHActsVertexPropagator.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHActsVertexPropagator.cc
3 
6 #include <phool/getClass.h>
7 #include <phool/phool.h>
8 #include <phool/PHDataNode.h>
9 #include <phool/PHNode.h>
10 #include <phool/PHNodeIterator.h>
11 #include <phool/PHObject.h>
12 #include <phool/PHTimer.h>
13 
18 
19 #include <Acts/Geometry/GeometryIdentifier.hpp>
25 
26 #include <ActsExamples/Plugins/BField/ScalableBField.hpp>
27 
29  : SubsysReco(name)
30  , m_trajectories(nullptr)
31 {}
32 
34 {
36 }
37 
39 {
40  int returnval = getNodes(topNode);
41  return returnval;
42 }
44 {
45 
46  if(m_vertexMap->size() == 0)
47  { setTrackVertexTo0(); }
48 
49  std::vector<unsigned int> deletedKeys;
50  for(const auto& [trackKey, trajectory] : *m_trajectories)
51  {
52  auto svtxTrack = m_trackMap->get(trackKey);
53  if(!svtxTrack)
54  {
57  deletedKeys.push_back(trackKey);
58  continue;
59  }
60 
61  if(Verbosity() > 2)
62  { svtxTrack->identify(); }
63 
64  const auto &[trackTips, mj] = trajectory.trajectory();
65 
66  if(trackTips.size() > 1 and Verbosity() > 0)
67  {
68  std::cout << PHWHERE
69  << "More than 1 track tip per track. Should never happen..."
70  << std::endl;
71  }
72 
73  for(const auto& trackTip : trackTips)
74  {
75  const auto& boundParams = trajectory.trackParameters(trackTip);
76 
77  auto propresult = propagateTrack(boundParams, svtxTrack->get_vertex_id());
78  if(propresult.ok())
79  {
80 
81  auto paramsAtVertex = std::move(**propresult);
82  updateSvtxTrack(svtxTrack,paramsAtVertex);
83  }
84  }
85  }
86 
88  for(auto& key : deletedKeys)
89  {
90  m_trajectories->erase(key);
91  }
92 
94 }
95 
97  const Acts::BoundTrackParameters& params)
98 {
99  auto position = params.position(m_tGeometry->geoContext);
100 
101  if(Verbosity() > 2)
102  {
103  std::cout << "Updating position track parameters from " << track->get_x()
104  << ", " << track->get_y() << ", " << track->get_z() << " to "
105  << position.transpose() / 10.
106  << std::endl;
107  }
108 
112 
113  ActsTransformations rotater;
114  rotater.setVerbosity(Verbosity());
115  if(params.covariance())
116  {
117  auto rotatedCov = rotater.rotateActsCovToSvtxTrack(params, m_tGeometry->geoContext);
118 
120  for(int i = 0; i < 3; i++) {
121  for(int j = 0; j < 3; j++) {
122  track->set_error(i,j, rotatedCov(i,j));
123  if(i < 2 and j < 2)
124  { track->set_acts_covariance(i,j, params.covariance().value()(i,j)); }
125  }
126  }
127  }
128 
129  updateTrackDCA(track);
130 
131 }
132 
134 {
135  Acts::Vector3D pos(track->get_x(),
136  track->get_y(),
137  track->get_z());
138  Acts::Vector3D mom(track->get_px(),
139  track->get_py(),
140  track->get_pz());
141 
142  auto vtxid = track->get_vertex_id();
143  auto svtxVertex = m_vertexMap->get(vtxid);
144  Acts::Vector3D vertex(svtxVertex->get_x(),
145  svtxVertex->get_y(),
146  svtxVertex->get_z());
147 
148  pos -= vertex;
149 
151  for(int i = 0; i < 3; ++i)
152  {
153  for(int j = 0; j < 3; ++j)
154  {
155  posCov(i, j) = track->get_error(i, j);
156  }
157  }
158 
159  Acts::Vector3D r = mom.cross(Acts::Vector3D(0.,0.,1.));
160  float phi = atan2(r(1), r(0));
161 
164  rot(0,0) = cos(phi);
165  rot(0,1) = -sin(phi);
166  rot(0,2) = 0;
167  rot(1,0) = sin(phi);
168  rot(1,1) = cos(phi);
169  rot(1,2) = 0;
170  rot(2,0) = 0;
171  rot(2,1) = 0;
172  rot(2,2) = 1;
173 
174  rot_T = rot.transpose();
175 
176  Acts::Vector3D pos_R = rot * pos;
177  Acts::ActsSymMatrixD<3> rotCov = rot * posCov * rot_T;
178 
179  const auto dca3Dxy = pos_R(0);
180  const auto dca3Dz = pos_R(2);
181  const auto dca3DxyCov = rotCov(0,0);
182  const auto dca3DzCov = rotCov(2,2);
183 
184  track->set_dca3d_xy(dca3Dxy);
185  track->set_dca3d_z(dca3Dz);
186  track->set_dca3d_xy_error(sqrt(dca3DxyCov));
187  track->set_dca3d_z_error(sqrt(dca3DzCov));
188 
189 }
191  const Acts::BoundTrackParameters& params,
192  const unsigned int vtxid)
193 {
194 
196  auto actsVertex = getVertex(vtxid);
197  auto perigee = Acts::Surface::makeShared<Acts::PerigeeSurface>(actsVertex);
198 
199  return std::visit(
200  [params, perigee, this]
201  (auto && inputField) ->BoundTrackParamPtrResult {
202  using InputMagneticField =
203  typename std::decay_t<decltype(inputField)>::element_type;
204  using MagneticField = Acts::SharedBField<InputMagneticField>;
207 
208  MagneticField field(inputField);
209  Stepper stepper(field);
210  Propagator propagator(stepper);
211 
213  if(Verbosity() > 3)
214  { logLevel = Acts::Logging::VERBOSE; }
215 
216  auto logger = Acts::getDefaultLogger("PHActsVertexPropagator",
217  logLevel);
218 
221  Acts::LoggerWrapper{*logger});
222 
223  auto result = propagator.propagate(params, *perigee,
224  options);
225  if(result.ok())
226  { return std::move((*result).endParameters); }
227 
228  return result.error();
229  },
230  m_tGeometry->magField);
231 }
232 
234 {
235  auto svtxVertex = m_vertexMap->get(vtxid);
236  return Acts::Vector3D(svtxVertex->get_x() * Acts::UnitConstants::cm,
237  svtxVertex->get_y() * Acts::UnitConstants::cm,
238  svtxVertex->get_z() * Acts::UnitConstants::cm);
239 }
240 
242 {
244  auto vertex = std::make_unique<SvtxVertex_v1>();
245  vertex->set_chisq(0.);
246  vertex->set_ndof(0);
247  vertex->set_t0(0);
248  vertex->set_id(0);
249  vertex->set_x(0);
250  vertex->set_y(0);
251  vertex->set_z(0);
252  for(int i=0; i<3; i++) {
253  for(int j=0; j<3; j++) {
254  vertex->set_error(i,j, 20.);
255  }
256  }
257 
258  m_vertexMap->insert(vertex.release());
259 
260  for(auto& [key, track] : *m_trackMap)
261  {
262  track->set_vertex_id(0);
263  }
264 
265 }
266 
268 {
270 }
271 
272 
274 {
275  m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode, "ActsTrackingGeometry");
276  if(!m_tGeometry)
277  {
278  std::cout << PHWHERE << "Acts tracking geometry not on node tree, exiting."
279  << std::endl;
281  }
282 
283  m_trajectories = findNode::getClass<std::map<const unsigned int, Trajectory>>(topNode, "ActsTrajectories");
284 
285  if(!m_trajectories)
286  {
287  std::cout << PHWHERE << "No acts trajectories on node tree, exiting. "
288  << std::endl;
290  }
291 
292  m_vertexMap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
293  if(!m_vertexMap)
294  {
295  std::cout << PHWHERE << "No svtx vertex map, exiting." << std::endl;
297  }
298 
299  m_trackMap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
300  if(!m_trackMap)
301  {
302  std::cout << PHWHERE << "No svtx track map, exiting. " << std::endl;
304  }
305 
307 }