ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ActsTransformations.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ActsTransformations.cc
1 #include "ActsTransformations.h"
2 #include "SvtxTrackState_v1.h"
4 
5 #include <chrono>
6 using namespace std::chrono;
7 
8 
9 namespace
10 {
11  template<class T> inline constexpr T square(const T& x) {return x*x;}
12  template<class T> T radius(const T& x, const T& y)
13  { return std::sqrt(square(x) + square(y));}
14 }
15 
17  const SvtxTrack *track,
18  Acts::GeometryContext /*geoCtxt*/) const
19 {
20  Acts::BoundSymMatrix svtxCovariance = Acts::BoundSymMatrix::Zero();
21 
22  for(int i = 0; i < 6; ++i)
23  {
24  for(int j = 0; j < 6; ++j)
25  {
26  svtxCovariance(i,j) = track->get_error(i,j);
28  if(i < 3 && j < 3)
29  svtxCovariance(i,j) *= Acts::UnitConstants::cm2;
30  else if (i < 3)
31  svtxCovariance(i,j) *= Acts::UnitConstants::cm;
32  else if (j < 3)
33  svtxCovariance(i,j) *= Acts::UnitConstants::cm;
34 
35  }
36  }
37 
38  printMatrix("svtx covariance, acts units: ", svtxCovariance);
39 
40  const double px = track->get_px();
41  const double py = track->get_py();
42  const double pz = track->get_pz();
43  const double p = sqrt(px*px + py*py + pz*pz);
44 
45  const double uPx = px / p;
46  const double uPy = py / p;
47  const double uPz = pz / p;
48 
49  //Acts version
50  const double cosTheta = uPz;
51  const double sinTheta = sqrt(uPx * uPx + uPy * uPy);
52  const double invSinTheta = 1. / sinTheta;
53  const double cosPhi = uPx * invSinTheta; // equivalent to x/r
54  const double sinPhi = uPy * invSinTheta; // equivalent to y/r
55 
60 
62  Acts::ActsMatrixD<6,8> sphenixRot;
63  sphenixRot.setZero();
65  sphenixRot(0,0) = 1;
66  sphenixRot(1,1) = 1;
67  sphenixRot(2,2) = 1;
68  sphenixRot(3,4) = 1./p;
69  sphenixRot(4,5) = 1./p;
70  sphenixRot(5,6) = 1./p;
71  sphenixRot(3,7) = uPx * p * p;
72  sphenixRot(4,7) = uPy * p * p;
73  sphenixRot(5,7) = uPz * p * p;
74 
75  auto rotatedMatrix
76  = sphenixRot.transpose() * svtxCovariance * sphenixRot;
77 
79  Acts::BoundToFreeMatrix jacobianLocalToGlobal = Acts::BoundToFreeMatrix::Zero();
80  jacobianLocalToGlobal(0, Acts::eBoundLoc0) = -sinPhi;
81  jacobianLocalToGlobal(0, Acts::eBoundLoc1) = -cosPhi * cosTheta;
82  jacobianLocalToGlobal(1, Acts::eBoundLoc0) = cosPhi;
83  jacobianLocalToGlobal(1, Acts::eBoundLoc1) = -sinPhi * cosTheta;
84  jacobianLocalToGlobal(2, Acts::eBoundLoc1) = sinTheta;
85  jacobianLocalToGlobal(3, Acts::eBoundTime) = 1;
86  jacobianLocalToGlobal(4, Acts::eBoundPhi) = -sinTheta * sinPhi;
87  jacobianLocalToGlobal(4, Acts::eBoundTheta) = cosTheta * cosPhi;
88  jacobianLocalToGlobal(5, Acts::eBoundPhi) = sinTheta * cosPhi;
89  jacobianLocalToGlobal(5, Acts::eBoundTheta) = cosTheta * sinPhi;
90  jacobianLocalToGlobal(6, Acts::eBoundTheta) = -sinTheta;
91  jacobianLocalToGlobal(7, Acts::eBoundQOverP) = 1;
92 
95  Acts::BoundSymMatrix actsLocalCov =
96  jacobianLocalToGlobal.transpose() * rotatedMatrix * jacobianLocalToGlobal;
97 
98  printMatrix("Rotated to Acts local cov : ",actsLocalCov);
99  return actsLocalCov;
100 
101 }
102 
103 
105  const Acts::BoundTrackParameters params,
106  Acts::GeometryContext /*geoCtxt*/) const
107 {
108 
109  auto covarianceMatrix = *params.covariance();
110 
111  printMatrix("Initial Acts covariance: ", covarianceMatrix);
112 
113  const double px = params.momentum()(0);
114  const double py = params.momentum()(1);
115  const double pz = params.momentum()(2);
116  const double p = params.momentum().norm();
117 
118  const double uPx = px / p;
119  const double uPy = py / p;
120  const double uPz = pz / p;
121 
122  const double cosTheta = uPz;
123  const double sinTheta = sqrt(uPx * uPx + uPy * uPy);
124  const double invSinTheta = 1. / sinTheta;
125  const double cosPhi = uPx * invSinTheta;
126  const double sinPhi = uPy * invSinTheta;
127 
128  Acts::BoundToFreeMatrix jacobianLocalToGlobal = Acts::BoundToFreeMatrix::Zero();
129  jacobianLocalToGlobal(0, Acts::eBoundLoc0) = -sinPhi;
130  jacobianLocalToGlobal(0, Acts::eBoundLoc1) = -cosPhi * cosTheta;
131  jacobianLocalToGlobal(1, Acts::eBoundLoc0) = cosPhi;
132  jacobianLocalToGlobal(1, Acts::eBoundLoc1) = -sinPhi * cosTheta;
133  jacobianLocalToGlobal(2, Acts::eBoundLoc1) = sinTheta;
134  jacobianLocalToGlobal(3, Acts::eBoundTime) = 1;
135  jacobianLocalToGlobal(4, Acts::eBoundPhi) = -sinTheta * sinPhi;
136  jacobianLocalToGlobal(4, Acts::eBoundTheta) = cosTheta * cosPhi;
137  jacobianLocalToGlobal(5, Acts::eBoundPhi) = sinTheta * cosPhi;
138  jacobianLocalToGlobal(5, Acts::eBoundTheta) = cosTheta * sinPhi;
139  jacobianLocalToGlobal(6, Acts::eBoundTheta) = -sinTheta;
140  jacobianLocalToGlobal(7, Acts::eBoundQOverP) = 1;
141 
143  auto rotatedMatrix
144  = jacobianLocalToGlobal * covarianceMatrix * jacobianLocalToGlobal.transpose();
145 
148  Acts::ActsMatrixD<6,8> sphenixRot;
149  sphenixRot.setZero();
150 
152  sphenixRot(0,0) = 1;
153  sphenixRot(1,1) = 1;
154  sphenixRot(2,2) = 1;
155  sphenixRot(3,4) = p;
156  sphenixRot(4,5) = p;
157  sphenixRot(5,6) = p;
158  sphenixRot(3,7) = uPx * p * p;
159  sphenixRot(4,7) = uPy * p * p;
160  sphenixRot(5,7) = uPz * p * p;
161 
162  Acts::BoundSymMatrix globalCov = Acts::BoundSymMatrix::Zero();
163  globalCov = sphenixRot * rotatedMatrix * sphenixRot.transpose();
164 
165  printMatrix("Global sPHENIX cov : ", globalCov);
166 
168  for(int i = 0; i < 6; ++i)
169  {
170  for(int j = 0; j < 6; ++j)
171  {
172  if(i < 3 && j < 3)
173  globalCov(i,j) /= Acts::UnitConstants::cm2;
174  else if (i < 3)
175  globalCov(i,j) /= Acts::UnitConstants::cm;
176  else if (j < 3)
177  globalCov(i,j) /= Acts::UnitConstants::cm;
178 
179  }
180  }
181 
182  printMatrix("Global sphenix cov after unit conv: " , globalCov);
183 
184  return globalCov;
185 }
186 
187 
189  Acts::BoundSymMatrix matrix) const
190 {
191 
192  if(m_verbosity > 10)
193  {
194  std::cout << std::endl << message.c_str() << std::endl;
195  for(int i = 0 ; i < matrix.rows(); ++i)
196  {
197  std::cout << std::endl;
198  for(int j = 0; j < matrix.cols(); ++j)
199  {
200  std::cout << matrix(i,j) << ", ";
201  }
202  }
203 
204  std::cout << std::endl;
205  }
206 }
207 
208 
209 
211  Acts::Vector3D vertex,
213  Acts::GeometryContext geoCtxt,
214  float &dca3Dxy,
215  float &dca3Dz,
216  float &dca3DxyCov,
217  float &dca3DzCov) const
218 {
219  Acts::Vector3D pos = param.position(geoCtxt);
220  Acts::Vector3D mom = param.momentum();
221 
223  pos -= vertex;
224 
226  for(int i = 0; i < 3; ++i)
227  {
228  for(int j = 0; j < 3; ++j)
229  {
230  posCov(i,j) = cov(i,j);
231  }
232  }
233 
234  Acts::Vector3D r = mom.cross(Acts::Vector3D(0.,0.,1.));
235  float phi = atan2(r(1), r(0));
236 
239  rot(0,0) = cos(phi);
240  rot(0,1) = -sin(phi);
241  rot(0,2) = 0;
242  rot(1,0) = sin(phi);
243  rot(1,1) = cos(phi);
244  rot(1,2) = 0;
245  rot(2,0) = 0;
246  rot(2,1) = 0;
247  rot(2,2) = 1;
248 
249  rot_T = rot.transpose();
250 
251  Acts::Vector3D pos_R = rot * pos;
252  Acts::ActsSymMatrixD<3> rotCov = rot * posCov * rot_T;
253 
254  dca3Dxy = pos_R(0);
255  dca3Dz = pos_R(2);
256  dca3DxyCov = rotCov(0,0);
257  dca3DzCov = rotCov(2,2);
258 
259 }
260 
262  ActsSurfaceMaps* surfMaps,
264 {
265  const Acts::Vector3D doublePos = getGlobalPosition(cluster,surfMaps,tGeometry);
266  return Acts::Vector3F(doublePos(0), doublePos(1), doublePos(2));
267 }
269  ActsSurfaceMaps* surfMaps,
271 {
272  Acts::Vector3D glob;
273 
274  const auto trkrid = TrkrDefs::getTrkrId(cluster->getClusKey());
275 
276  auto surface = getSurface(cluster, surfMaps);
277  if(!surface)
278  {
279  std::cerr << "Couldn't identify cluster surface. Returning NAN"
280  << std::endl;
281  glob(0) = NAN;
282  glob(1) = NAN;
283  glob(2) = NAN;
284  return glob;
285  }
286 
287  Acts::Vector2D local(cluster->getLocalX(), cluster->getLocalY());
288  Acts::Vector3D global;
290  if(trkrid != TrkrDefs::tpcId)
291  {
292  global = surface->localToGlobal(tGeometry->geoContext,
294  Acts::Vector3D(1,1,1));
295  global /= Acts::UnitConstants::cm;
296  return global;
297  }
298 
301  auto surfCenter = surface->center(tGeometry->geoContext);
302  surfCenter /= Acts::UnitConstants::cm;
303  double surfPhiCenter = atan2(surfCenter(1), surfCenter(0));
304  double surfRadius = radius(surfCenter(0), surfCenter(1));
305  double surfRPhiCenter = surfPhiCenter * surfRadius;
306 
307  double clusRPhi = local(0) + surfRPhiCenter;
308  double gclusz = local(1) + surfCenter(2);
309 
310  double clusphi = clusRPhi / surfRadius;
311  double gclusx = surfRadius * cos(clusphi);
312  double gclusy = surfRadius * sin(clusphi);
313  global(0) = gclusx;
314  global(1) = gclusy;
315  global(2) = gclusz;
316 
317  return global;
318 }
319 
321  ActsSurfaceMaps *surfMaps) const
322 {
323  const auto cluskey = cluster->getClusKey();
324  const auto surfkey = cluster->getSubSurfKey();
325  const auto trkrid = TrkrDefs::getTrkrId(cluskey);
327 
328  switch( trkrid )
329  {
330  case TrkrDefs::TrkrId::micromegasId: return getMMSurface( hitsetkey, surfMaps );
331  case TrkrDefs::TrkrId::tpcId: return getTpcSurface(hitsetkey, surfkey, surfMaps);
334  {
335  return getSiliconSurface(hitsetkey, surfMaps);
336  }
337  }
338 
339  // unreachable
340  return nullptr;
341 }
342 
344  ActsSurfaceMaps *maps) const
345 {
346  auto surfMap = maps->siliconSurfaceMap;
347  auto iter = surfMap.find(hitsetkey);
348  if(iter != surfMap.end())
349  {
350  return iter->second;
351  }
352 
354  return nullptr;
355 
356 }
357 
359  TrkrDefs::subsurfkey surfkey,
360  ActsSurfaceMaps* maps) const
361 {
362  unsigned int layer = TrkrDefs::getLayer(hitsetkey);
363  const auto iter = maps->tpcSurfaceMap.find(layer);
364  if(iter != maps->tpcSurfaceMap.end())
365  {
366  auto surfvec = iter->second;
367  return surfvec.at(surfkey);
368  }
369 
371  return nullptr;
372 }
373 
374 
376  ActsSurfaceMaps *maps) const
377 {
378  const auto iter = maps->mmSurfaceMap.find( hitsetkey );
379  return (iter == maps->mmSurfaceMap.end()) ? nullptr:iter->second;
380 }
381 
383  const size_t &trackTip,
384  SvtxTrack *svtxTrack,
386 {
387 
388  const auto &[trackTips, mj] = traj.trajectory();
389 
390  mj.visitBackwards(trackTip, [&](const auto &state)
391  {
392 
394  const auto typeFlags = state.typeFlags();
395  if( !typeFlags.test(Acts::TrackStateFlag::MeasurementFlag) )
396  { return true; }
397 
398  // only fill for state vectors with proper smoothed parameters
399  if( !state.hasSmoothed()) return true;
400 
401  // create svtx state vector with relevant pathlength
402  const float pathlength = state.pathLength() / Acts::UnitConstants::cm;
403  SvtxTrackState_v1 out( pathlength );
404 
405  // get smoothed fitted parameters
406  const Acts::BoundTrackParameters params(state.referenceSurface().getSharedPtr(),
407  state.smoothed(),
408  state.smoothedCovariance());
409 
410  // position
411  const auto global = params.position(geoContext);
412  out.set_x(global.x() / Acts::UnitConstants::cm);
413  out.set_y(global.y() / Acts::UnitConstants::cm);
414  out.set_z(global.z() / Acts::UnitConstants::cm);
415 
416  // momentum
417  const auto momentum = params.momentum();
418  out.set_px(momentum.x());
419  out.set_py(momentum.y());
420  out.set_pz(momentum.z());
421 
423  const auto globalCov = rotateActsCovToSvtxTrack(params, geoContext);
424  for (int i = 0; i < 6; ++i)
425  for (int j = 0; j < 6; ++j)
426  { out.set_error(i, j, globalCov(i,j)); }
427 
428  // print
429  if(m_verbosity > 20)
430  {
431  std::cout << " inserting state with x,y,z ="
432  << " " << global.x() / Acts::UnitConstants::cm
433  << " " << global.y() / Acts::UnitConstants::cm
434  << " " << global.z() / Acts::UnitConstants::cm
435  << " pathlength " << pathlength
436  << " momentum px,py,pz = " << momentum.x() << " " << momentum.y() << " " << momentum.y()
437  << std::endl
438  << "covariance " << globalCov << std::endl;
439  }
440 
441  svtxTrack->insert_state(&out);
442 
443  return true;
444  }
445  );
446 
447  return;
448 }