Analysis Software
Documentation for sPHENIX simulation 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 "SvtxTrack.h"
3 #include "SvtxTrackState.h"
4 #include "SvtxTrackState_v2.h"
5 
6 #include <trackbase/TrkrCluster.h>
7 #include <trackbase/InttDefs.h>
8 #include <trackbase/MvtxDefs.h>
9 #include <trackbase/ActsSourceLink.h>
10 
11 #include <Acts/EventData/ParticleHypothesis.hpp>
12 
13 namespace
14 {
16  template<class T> inline constexpr T square(const T& x) {return x*x;}
17 
19  template<class T> T radius(const T& x, const T& y)
20  { return std::sqrt(square(x) + square(y));}
21 
23  template<class T>
24  void print_matrix( const std::string &message, const T& matrix)
25  {
26  std::cout << std::endl << message << std::endl;
27  for(int i = 0 ; i < matrix.rows(); ++i)
28  {
29  for(int j = 0; j < matrix.cols(); ++j)
30  { std::cout << matrix(i,j) << ", "; }
31 
32  std::cout << std::endl;
33  }
34  }
35 
36 }
37 
38 //_______________________________________________________________________________
39 Acts::BoundSquareMatrix ActsTransformations::rotateSvtxTrackCovToActs( const SvtxTrack *track ) const
40 { return rotateSvtxTrackCovToActs( track->find_state(0.0)->second ); }
41 
42 //_______________________________________________________________________________
43 Acts::BoundSquareMatrix ActsTransformations::rotateSvtxTrackCovToActs(
44  const SvtxTrackState *state) const
45 {
46  Acts::BoundSquareMatrix svtxCovariance = Acts::BoundSquareMatrix::Zero();
47 
48  for(int i = 0; i < 6; ++i)
49  {
50  for(int j = 0; j < 6; ++j)
51  {
52  svtxCovariance(i,j) = state->get_error(i,j);
54  if(i < 3 && j < 3)
55  svtxCovariance(i,j) *= Acts::UnitConstants::cm2;
56  else if (i < 3)
57  svtxCovariance(i,j) *= Acts::UnitConstants::cm;
58  else if (j < 3)
59  svtxCovariance(i,j) *= Acts::UnitConstants::cm;
60 
61  }
62  }
63 
64  printMatrix("svtx covariance, acts units: ", svtxCovariance);
65 
66  const double px = state->get_px();
67  const double py = state->get_py();
68  const double pz = state->get_pz();
69  const double p2 = square(px) + square(py) + square(pz);
70  const double p = std::sqrt(p2);
71  const double invp = 1./p;
72 
73  const double uPx = px/p;
74  const double uPy = py/p;
75  const double uPz = pz/p;
76 
77  //Acts version
78  const double cosTheta = uPz;
79  const double sinTheta = std::sqrt(square(uPx) + square(uPy));
80  const double invSinTheta = 1. / sinTheta;
81  const double cosPhi = uPx * invSinTheta; // equivalent to x/r
82  const double sinPhi = uPy * invSinTheta; // equivalent to y/r
83 
88 
89  // rotate to (x,y,z,t,upx,upy,upz,q/p)
90  Acts::ActsMatrix<8,6> sphenixRot = Acts::ActsMatrix<8,6>::Zero();
91  sphenixRot(0,0) = 1;
92  sphenixRot(1,1) = 1;
93  sphenixRot(2,2) = 1;
94  sphenixRot(4,3) = invp;
95  sphenixRot(5,4) = invp;
96  sphenixRot(6,5) = invp;
97  sphenixRot(7,3) = uPx/p2;
98  sphenixRot(7,4) = uPy/p2;
99  sphenixRot(7,5) = uPz/p2;
100 
101  const auto rotatedMatrix = sphenixRot*svtxCovariance*sphenixRot.transpose();
102 
103  // global to local transformation
104  /* from https://github.com/acts-project/acts/blob/394cfdb308956de93f90ab3162040bb6d835027d/Core/src/Propagator/detail/CovarianceEngine.cpp#L31 */
105  Acts::FreeToBoundMatrix jacobianGlobalToLocal = Acts::FreeToBoundMatrix::Zero();
106  jacobianGlobalToLocal(Acts::eBoundLoc0, 0) = -sinPhi;
107  jacobianGlobalToLocal(Acts::eBoundLoc0, 1) = cosPhi;
108  jacobianGlobalToLocal(Acts::eBoundLoc1, 0) = -cosPhi * cosTheta;
109  jacobianGlobalToLocal(Acts::eBoundLoc1, 1) = -sinPhi * cosTheta;
110  jacobianGlobalToLocal(Acts::eBoundLoc1, 2) = sinTheta;
111  jacobianGlobalToLocal(Acts::eBoundTime, 3) = 1.;
112  jacobianGlobalToLocal(Acts::eBoundPhi, 4) = -sinPhi * invSinTheta;
113  jacobianGlobalToLocal(Acts::eBoundPhi, 5) = cosPhi * invSinTheta;
114  jacobianGlobalToLocal(Acts::eBoundTheta, 4) = cosPhi * cosTheta;
115  jacobianGlobalToLocal(Acts::eBoundTheta, 5) = sinPhi * cosTheta;
116  jacobianGlobalToLocal(Acts::eBoundTheta, 6) = -sinTheta;
117  jacobianGlobalToLocal(Acts::eBoundQOverP, 7) = 1.;
118 
119  Acts::BoundSquareMatrix actsLocalCov = jacobianGlobalToLocal * rotatedMatrix * jacobianGlobalToLocal.transpose();
120 
121  /*
122  * need to assign the covariance matrix diagonal element corresponding to the time coordinate manually,
123  * since it is lost in the conversion to Svtx coordinates
124  */
125  actsLocalCov( Acts::eBoundTime, Acts::eBoundTime ) = 1.;
126 
127  printMatrix("Rotated to Acts local cov : ",actsLocalCov);
128  return actsLocalCov;
129 
130 }
131 
132 //_______________________________________________________________________________
133 Acts::BoundSquareMatrix ActsTransformations::rotateActsCovToSvtxTrack( const ActsTrackFittingAlgorithm::TrackParameters& params ) const
134 {
135  const auto covarianceMatrix = *params.covariance();
136  printMatrix("Initial Acts covariance: ", covarianceMatrix);
137 
138  const double px = params.momentum().x();
139  const double py = params.momentum().y();
140  const double pz = params.momentum().z();
141  const double p = params.momentum().norm();
142  const double p2 = square(p);
143 
144  const double uPx = px/p;
145  const double uPy = py/p;
146  const double uPz = pz/p;
147 
148  const double cosTheta = uPz;
149  const double sinTheta = std::sqrt(square(uPx) + square(uPy));
150  const double invSinTheta = 1./sinTheta;
151  const double cosPhi = uPx*invSinTheta;
152  const double sinPhi = uPy*invSinTheta;
153 
154  // local to global transformation
155  /* from https://github.com/acts-project/acts/blob/394cfdb308956de93f90ab3162040bb6d835027d/Core/src/Propagator/detail/CovarianceEngine.cpp#L222 */
156  Acts::BoundToFreeMatrix jacobianLocalToGlobal = Acts::BoundToFreeMatrix::Zero();
157  jacobianLocalToGlobal(0, Acts::eBoundLoc0) = -sinPhi;
158  jacobianLocalToGlobal(0, Acts::eBoundLoc1) = -cosPhi * cosTheta;
159  jacobianLocalToGlobal(1, Acts::eBoundLoc0) = cosPhi;
160  jacobianLocalToGlobal(1, Acts::eBoundLoc1) = -sinPhi * cosTheta;
161  jacobianLocalToGlobal(2, Acts::eBoundLoc1) = sinTheta;
162  jacobianLocalToGlobal(3, Acts::eBoundTime) = 1;
163  jacobianLocalToGlobal(4, Acts::eBoundPhi) = -sinTheta * sinPhi;
164  jacobianLocalToGlobal(4, Acts::eBoundTheta) = cosTheta * cosPhi;
165  jacobianLocalToGlobal(5, Acts::eBoundPhi) = sinTheta * cosPhi;
166  jacobianLocalToGlobal(5, Acts::eBoundTheta) = cosTheta * sinPhi;
167  jacobianLocalToGlobal(6, Acts::eBoundTheta) = -sinTheta;
168  jacobianLocalToGlobal(7, Acts::eBoundQOverP) = 1;
169 
170  // Covariance is now an 8x8 matrix in basis (x,y,z,time,Tx,Ty,Tz,q/p)
171  const auto rotatedMatrix = jacobianLocalToGlobal * covarianceMatrix * jacobianLocalToGlobal.transpose();
172 
173  // Now rotate to x,y,z, px,py,pz
174  Acts::ActsMatrix<6,8> sphenixRot = Acts::ActsMatrix<6,8>::Zero();
175  sphenixRot(0,0) = 1;
176  sphenixRot(1,1) = 1;
177  sphenixRot(2,2) = 1;
178  sphenixRot(3,4) = p;
179  sphenixRot(4,5) = p;
180  sphenixRot(5,6) = p;
181  sphenixRot(3,7) = uPx*p2;
182  sphenixRot(4,7) = uPy*p2;
183  sphenixRot(5,7) = uPz*p2;
184 
185  Acts::BoundSquareMatrix globalCov = sphenixRot * rotatedMatrix * sphenixRot.transpose();
186  printMatrix("Global sPHENIX cov : ", globalCov);
187 
189  for(int i = 0; i < 6; ++i)
190  {
191  for(int j = 0; j < 6; ++j)
192  {
193  if(i < 3 && j < 3)
194  globalCov(i,j) /= Acts::UnitConstants::cm2;
195  else if (i < 3)
196  globalCov(i,j) /= Acts::UnitConstants::cm;
197  else if (j < 3)
198  globalCov(i,j) /= Acts::UnitConstants::cm;
199 
200  }
201  }
202 
203  printMatrix("Global sphenix cov after unit conv: " , globalCov);
204 
205  return globalCov;
206 }
207 
208 void ActsTransformations::printMatrix(const std::string &message, const Acts::BoundSquareMatrix& matrix) const
209 { if(m_verbosity > 10) print_matrix( message, matrix ); }
210 
211 void ActsTransformations::calculateDCA(const Acts::BoundTrackParameters param,
212  Acts::Vector3 vertex,
213  Acts::BoundSquareMatrix cov,
214  Acts::GeometryContext& geoCtxt,
215  float &dca3Dxy,
216  float &dca3Dz,
217  float &dca3DxyCov,
218  float &dca3DzCov) const
219 {
220  Acts::Vector3 pos = param.position(geoCtxt);
221  Acts::Vector3 mom = param.momentum();
222 
224  pos -= vertex;
225 
226  Acts::ActsSquareMatrix<3> posCov;
227  for(int i = 0; i < 3; ++i)
228  {
229  for(int j = 0; j < 3; ++j)
230  {
231  posCov(i,j) = cov(i,j);
232  }
233  }
234 
235  Acts::Vector3 r = mom.cross(Acts::Vector3(0.,0.,1.));
236  float phi = atan2(r(1), r(0));
237 
238  Acts::RotationMatrix3 rot;
239  Acts::RotationMatrix3 rot_T;
240  rot(0,0) = cos(phi);
241  rot(0,1) = -sin(phi);
242  rot(0,2) = 0;
243  rot(1,0) = sin(phi);
244  rot(1,1) = cos(phi);
245  rot(1,2) = 0;
246  rot(2,0) = 0;
247  rot(2,1) = 0;
248  rot(2,2) = 1;
249 
250  rot_T = rot.transpose();
251 
252  Acts::Vector3 pos_R = rot * pos;
253  Acts::ActsSquareMatrix<3> rotCov = rot * posCov * rot_T;
254 
255  dca3Dxy = pos_R(0);
256  dca3Dz = pos_R(2);
257  dca3DxyCov = rotCov(0,0);
258  dca3DzCov = rotCov(2,2);
259 
260 }
261 
262 
263 void ActsTransformations::fillSvtxTrackStates(const Acts::ConstVectorMultiTrajectory& traj,
264  const size_t& trackTip,
265  SvtxTrack *svtxTrack,
266  Acts::GeometryContext& geoContext) const
267 {
268  traj.visitBackwards(trackTip, [&](const auto &state)
269  {
270 
272  const auto typeFlags = state.typeFlags();
273  if( !typeFlags.test(Acts::TrackStateFlag::MeasurementFlag) )
274  { return true; }
275 
276  // only fill for state vectors with proper smoothed parameters
277  if( !state.hasSmoothed()) return true;
278 
279  // create svtx state vector with relevant pathlength
280  const float pathlength = state.pathLength() / Acts::UnitConstants::cm;
281  SvtxTrackState_v2 out( pathlength );
282 
283  // get smoothed fitted parameters
284  const Acts::BoundTrackParameters params(
285  state.referenceSurface().getSharedPtr(),
286  state.smoothed(),
287  state.smoothedCovariance(),
288  Acts::ParticleHypothesis::pion());
289 
290  // position
291  const auto global = params.position(geoContext);
292  out.set_x(global.x() / Acts::UnitConstants::cm);
293  out.set_y(global.y() / Acts::UnitConstants::cm);
294  out.set_z(global.z() / Acts::UnitConstants::cm);
295 
296  // momentum
297  const auto momentum = params.momentum();
298  out.set_px(momentum.x());
299  out.set_py(momentum.y());
300  out.set_pz(momentum.z());
301  auto sourceLink = state.getUncalibratedSourceLink().template get<ActsSourceLink>();
302  out.set_cluskey(sourceLink.cluskey());
303 
305  const auto globalCov = rotateActsCovToSvtxTrack(params);
306  for (int i = 0; i < 6; ++i)
307  for (int j = 0; j < 6; ++j)
308  { out.set_error(i, j, globalCov(i,j)); }
309 
310  // print
311  if(m_verbosity > 20)
312  {
313  std::cout << " inserting state with x,y,z ="
314  << " " << global.x() / Acts::UnitConstants::cm
315  << " " << global.y() / Acts::UnitConstants::cm
316  << " " << global.z() / Acts::UnitConstants::cm
317  << " pathlength " << pathlength
318  << " momentum px,py,pz = " << momentum.x() << " " << momentum.y() << " " << momentum.y()
319  << std::endl
320  << "covariance " << globalCov << std::endl;
321  }
322 
323  svtxTrack->insert_state(&out);
324 
325  return true;
326  }
327  );
328 
329  return;
330 }