CMS 3D CMS Logo
BrokenLine.h
Go to the documentation of this file.
1 #ifndef RecoPixelVertexing_PixelTrackFitting_interface_BrokenLine_h
2 #define RecoPixelVertexing_PixelTrackFitting_interface_BrokenLine_h
3 
4 #include <Eigen/Eigenvalues>
5 
6 #include "RecoPixelVertexing/PixelTrackFitting/interface/FitUtils.h"
7 
8 namespace brokenline {
9 
11 
18  using karimaki_circle_fit = riemannFit::CircleFit;
19 
23  template <int n>
24  struct PreparedBrokenLineData {
25  int qCharge;
26  riemannFit::Matrix2xNd<n> radii;
27  riemannFit::VectorNd<n> sTransverse;
28  // starting from the pre-fitted closest approach
29  riemannFit::VectorNd<n> sTotal;
30  riemannFit::VectorNd<n> zInSZplane;
31  riemannFit::VectorNd<n> varBeta;
32  };
33 
53  __host__ __device__ inline double multScatt(
54  const double& length, const double bField, const double radius, int layer, double slope) {
55  // limit R to 20GeV...
56  auto pt2 = std::min(20., bField * radius);
57  pt2 *= pt2;
58  constexpr double inv_X0 = 0.06 / 16.;
59  //if(Layer==1) XXI_0=0.06/16.;
60  // else XXI_0=0.06/16.;
61  //XX_0*=1;
62 
64  constexpr double geometry_factor = 0.7;
65  constexpr double fact = geometry_factor * riemannFit::sqr(13.6 / 1000.);
66  return fact / (pt2 * (1. + riemannFit::sqr(slope))) * (std::abs(length) * inv_X0) *
67  riemannFit::sqr(1. + 0.038 * log(std::abs(length) * inv_X0));
68  }
69 
77  __host__ __device__ inline riemannFit::Matrix2d rotationMatrix(double slope) {
78  riemannFit::Matrix2d rot;
79  rot(0, 0) = 1. / sqrt(1. + riemannFit::sqr(slope));
80  rot(0, 1) = slope * rot(0, 0);
81  rot(1, 0) = -rot(0, 1);
82  rot(1, 1) = rot(0, 0);
83  return rot;
84  }
85 
97  __host__ __device__ inline void translateKarimaki(karimaki_circle_fit& circle,
98  double x0,
99  double y0,
100  riemannFit::Matrix3d& jacobian) {
101  // Avoid multiple access to the circle.par vector.
102  using scalar = std::remove_reference<decltype(circle.par(0))>::type;
103  scalar phi = circle.par(0);
104  scalar dee = circle.par(1);
105  scalar rho = circle.par(2);
106 
107  // Avoid repeated trig. computations
108  scalar sinPhi = sin(phi);
109  scalar cosPhi = cos(phi);
110 
111  // Intermediate computations for the circle parameters
112  scalar deltaPara = x0 * cosPhi + y0 * sinPhi;
113  scalar deltaOrth = x0 * sinPhi - y0 * cosPhi + dee;
114  scalar tempSmallU = 1 + rho * dee;
115  scalar tempC = -rho * y0 + tempSmallU * cosPhi;
116  scalar tempB = rho * x0 + tempSmallU * sinPhi;
117  scalar tempA = 2. * deltaOrth + rho * (riemannFit::sqr(deltaOrth) + riemannFit::sqr(deltaPara));
118  scalar tempU = sqrt(1. + rho * tempA);
119 
120  // Intermediate computations for the error matrix transform
121  scalar xi = 1. / (riemannFit::sqr(tempB) + riemannFit::sqr(tempC));
122  scalar tempV = 1. + rho * deltaOrth;
123  scalar lambda = (0.5 * tempA) / (riemannFit::sqr(1. + tempU) * tempU);
124  scalar mu = 1. / (tempU * (1. + tempU)) + rho * lambda;
125  scalar zeta = riemannFit::sqr(deltaOrth) + riemannFit::sqr(deltaPara);
126  jacobian << xi * tempSmallU * tempV, -xi * riemannFit::sqr(rho) * deltaOrth, xi * deltaPara,
127  2. * mu * tempSmallU * deltaPara, 2. * mu * tempV, mu * zeta - lambda * tempA, 0, 0, 1.;
128 
129  // translated circle parameters
130  // phi
131  circle.par(0) = atan2(tempB, tempC);
132  // d
133  circle.par(1) = tempA / (1 + tempU);
134  // rho after translation. It is invariant, so noop
135  // circle.par(2)= rho;
136 
137  // translated error matrix
138  circle.cov = jacobian * circle.cov * jacobian.transpose();
139  }
140 
149  template <typename M3xN, typename V4, int n>
150  __host__ __device__ inline void prepareBrokenLineData(const M3xN& hits,
151  const V4& fast_fit,
152  const double bField,
153  PreparedBrokenLineData<n>& results) {
154  riemannFit::Vector2d dVec;
155  riemannFit::Vector2d eVec;
156 
157  dVec = hits.block(0, 1, 2, 1) - hits.block(0, 0, 2, 1);
158  eVec = hits.block(0, n - 1, 2, 1) - hits.block(0, n - 2, 2, 1);
159  results.qCharge = riemannFit::cross2D(dVec, eVec) > 0 ? -1 : 1;
160 
161  const double slope = -results.qCharge / fast_fit(3);
162 
163  riemannFit::Matrix2d rotMat = rotationMatrix(slope);
164 
165  // calculate radii and s
166  results.radii = hits.block(0, 0, 2, n) - fast_fit.head(2) * riemannFit::MatrixXd::Constant(1, n, 1);
167  eVec = -fast_fit(2) * fast_fit.head(2) / fast_fit.head(2).norm();
168  for (u_int i = 0; i < n; i++) {
169  dVec = results.radii.block(0, i, 2, 1);
170  results.sTransverse(i) = results.qCharge * fast_fit(2) *
171  atan2(riemannFit::cross2D(dVec, eVec), dVec.dot(eVec)); // calculates the arc length
172  }
173  riemannFit::VectorNd<n> zVec = hits.block(2, 0, 1, n).transpose();
174 
175  //calculate sTotal and zVec
176  riemannFit::Matrix2xNd<n> pointsSZ = riemannFit::Matrix2xNd<n>::Zero();
177  for (u_int i = 0; i < n; i++) {
178  pointsSZ(0, i) = results.sTransverse(i);
179  pointsSZ(1, i) = zVec(i);
180  pointsSZ.block(0, i, 2, 1) = rotMat * pointsSZ.block(0, i, 2, 1);
181  }
182  results.sTotal = pointsSZ.block(0, 0, 1, n).transpose();
183  results.zInSZplane = pointsSZ.block(1, 0, 1, n).transpose();
184 
185  //calculate varBeta
186  results.varBeta(0) = results.varBeta(n - 1) = 0;
187  for (u_int i = 1; i < n - 1; i++) {
188  results.varBeta(i) = multScatt(results.sTotal(i + 1) - results.sTotal(i), bField, fast_fit(2), i + 2, slope) +
189  multScatt(results.sTotal(i) - results.sTotal(i - 1), bField, fast_fit(2), i + 1, slope);
190  }
191  }
192 
205  template <int n>
206  __host__ __device__ inline riemannFit::MatrixNd<n> matrixC_u(const riemannFit::VectorNd<n>& weights,
207  const riemannFit::VectorNd<n>& sTotal,
208  const riemannFit::VectorNd<n>& varBeta) {
209  riemannFit::MatrixNd<n> c_uMat = riemannFit::MatrixNd<n>::Zero();
210  for (u_int i = 0; i < n; i++) {
211  c_uMat(i, i) = weights(i);
212  if (i > 1)
213  c_uMat(i, i) += 1. / (varBeta(i - 1) * riemannFit::sqr(sTotal(i) - sTotal(i - 1)));
214  if (i > 0 && i < n - 1)
215  c_uMat(i, i) +=
216  (1. / varBeta(i)) * riemannFit::sqr((sTotal(i + 1) - sTotal(i - 1)) /
217  ((sTotal(i + 1) - sTotal(i)) * (sTotal(i) - sTotal(i - 1))));
218  if (i < n - 2)
219  c_uMat(i, i) += 1. / (varBeta(i + 1) * riemannFit::sqr(sTotal(i + 1) - sTotal(i)));
220 
221  if (i > 0 && i < n - 1)
222  c_uMat(i, i + 1) =
223  1. / (varBeta(i) * (sTotal(i + 1) - sTotal(i))) *
224  (-(sTotal(i + 1) - sTotal(i - 1)) / ((sTotal(i + 1) - sTotal(i)) * (sTotal(i) - sTotal(i - 1))));
225  if (i < n - 2)
226  c_uMat(i, i + 1) +=
227  1. / (varBeta(i + 1) * (sTotal(i + 1) - sTotal(i))) *
228  (-(sTotal(i + 2) - sTotal(i)) / ((sTotal(i + 2) - sTotal(i + 1)) * (sTotal(i + 1) - sTotal(i))));
229 
230  if (i < n - 2)
231  c_uMat(i, i + 2) = 1. / (varBeta(i + 1) * (sTotal(i + 2) - sTotal(i + 1)) * (sTotal(i + 1) - sTotal(i)));
232 
233  c_uMat(i, i) *= 0.5;
234  }
235  return c_uMat + c_uMat.transpose();
236  }
237 
248  template <typename M3xN, typename V4>
249  __host__ __device__ inline void fastFit(const M3xN& hits, V4& result) {
250  constexpr uint32_t n = M3xN::ColsAtCompileTime;
251 
252  const riemannFit::Vector2d a = hits.block(0, n / 2, 2, 1) - hits.block(0, 0, 2, 1);
253  const riemannFit::Vector2d b = hits.block(0, n - 1, 2, 1) - hits.block(0, n / 2, 2, 1);
254  const riemannFit::Vector2d c = hits.block(0, 0, 2, 1) - hits.block(0, n - 1, 2, 1);
255 
256  auto tmp = 0.5 / riemannFit::cross2D(c, a);
257  result(0) = hits(0, 0) - (a(1) * c.squaredNorm() + c(1) * a.squaredNorm()) * tmp;
258  result(1) = hits(1, 0) + (a(0) * c.squaredNorm() + c(0) * a.squaredNorm()) * tmp;
259  // check Wikipedia for these formulas
260 
261  result(2) = sqrt(a.squaredNorm() * b.squaredNorm() * c.squaredNorm()) / (2. * std::abs(riemannFit::cross2D(b, a)));
262  // Using Math Olympiad's formula R=abc/(4A)
263 
264  const riemannFit::Vector2d d = hits.block(0, 0, 2, 1) - result.head(2);
265  const riemannFit::Vector2d e = hits.block(0, n - 1, 2, 1) - result.head(2);
266 
267  result(3) = result(2) * atan2(riemannFit::cross2D(d, e), d.dot(e)) / (hits(2, n - 1) - hits(2, 0));
268  // ds/dz slope between last and first point
269  }
270 
295  template <typename M3xN, typename M6xN, typename V4, int n>
296  __host__ __device__ inline void circleFit(const M3xN& hits,
297  const M6xN& hits_ge,
298  const V4& fast_fit,
299  const double bField,
300  PreparedBrokenLineData<n>& data,
301  karimaki_circle_fit& circle_results) {
302  circle_results.qCharge = data.qCharge;
303  auto& radii = data.radii;
304  const auto& sTransverse = data.sTransverse;
305  const auto& sTotal = data.sTotal;
306  auto& zInSZplane = data.zInSZplane;
307  auto& varBeta = data.varBeta;
308  const double slope = -circle_results.qCharge / fast_fit(3);
309  varBeta *= 1. + riemannFit::sqr(slope); // the kink angles are projected!
310 
311  for (u_int i = 0; i < n; i++) {
312  zInSZplane(i) = radii.block(0, i, 2, 1).norm() - fast_fit(2);
313  }
314 
315  riemannFit::Matrix2d vMat; // covariance matrix
316  riemannFit::VectorNd<n> weightsVec; // weights
317  riemannFit::Matrix2d rotMat; // rotation matrix point by point
318  for (u_int i = 0; i < n; i++) {
319  vMat(0, 0) = hits_ge.col(i)[0]; // x errors
320  vMat(0, 1) = vMat(1, 0) = hits_ge.col(i)[1]; // cov_xy
321  vMat(1, 1) = hits_ge.col(i)[2]; // y errors
322  rotMat = rotationMatrix(-radii(0, i) / radii(1, i));
323  weightsVec(i) =
324  1. / ((rotMat * vMat * rotMat.transpose())(1, 1)); // compute the orthogonal weight point by point
325  }
326 
327  riemannFit::VectorNplusONEd<n> r_uVec;
328  r_uVec(n) = 0;
329  for (u_int i = 0; i < n; i++) {
330  r_uVec(i) = weightsVec(i) * zInSZplane(i);
331  }
332 
333  riemannFit::MatrixNplusONEd<n> c_uMat;
334  c_uMat.block(0, 0, n, n) = matrixC_u(weightsVec, sTransverse, varBeta);
335  c_uMat(n, n) = 0;
336  //add the border to the c_uMat matrix
337  for (u_int i = 0; i < n; i++) {
338  c_uMat(i, n) = 0;
339  if (i > 0 && i < n - 1) {
340  c_uMat(i, n) +=
341  -(sTransverse(i + 1) - sTransverse(i - 1)) * (sTransverse(i + 1) - sTransverse(i - 1)) /
342  (2. * varBeta(i) * (sTransverse(i + 1) - sTransverse(i)) * (sTransverse(i) - sTransverse(i - 1)));
343  }
344  if (i > 1) {
345  c_uMat(i, n) +=
346  (sTransverse(i) - sTransverse(i - 2)) / (2. * varBeta(i - 1) * (sTransverse(i) - sTransverse(i - 1)));
347  }
348  if (i < n - 2) {
349  c_uMat(i, n) +=
350  (sTransverse(i + 2) - sTransverse(i)) / (2. * varBeta(i + 1) * (sTransverse(i + 1) - sTransverse(i)));
351  }
352  c_uMat(n, i) = c_uMat(i, n);
353  if (i > 0 && i < n - 1)
354  c_uMat(n, n) += riemannFit::sqr(sTransverse(i + 1) - sTransverse(i - 1)) / (4. * varBeta(i));
355  }
356 
357 #ifdef CPP_DUMP
358  std::cout << "CU5\n" << c_uMat << std::endl;
359 #endif
360  riemannFit::MatrixNplusONEd<n> iMat;
361  math::cholesky::invert(c_uMat, iMat);
362 #ifdef CPP_DUMP
363  std::cout << "I5\n" << iMat << std::endl;
364 #endif
365 
366  riemannFit::VectorNplusONEd<n> uVec = iMat * r_uVec; // obtain the fitted parameters by solving the linear system
367 
368  // compute (phi, d_ca, k) in the system in which the midpoint of the first two corrected hits is the origin...
369 
370  radii.block(0, 0, 2, 1) /= radii.block(0, 0, 2, 1).norm();
371  radii.block(0, 1, 2, 1) /= radii.block(0, 1, 2, 1).norm();
372 
373  riemannFit::Vector2d dVec = hits.block(0, 0, 2, 1) + (-zInSZplane(0) + uVec(0)) * radii.block(0, 0, 2, 1);
374  riemannFit::Vector2d eVec = hits.block(0, 1, 2, 1) + (-zInSZplane(1) + uVec(1)) * radii.block(0, 1, 2, 1);
375 
376  circle_results.par << atan2((eVec - dVec)(1), (eVec - dVec)(0)),
377  -circle_results.qCharge *
378  (fast_fit(2) - sqrt(riemannFit::sqr(fast_fit(2)) - 0.25 * (eVec - dVec).squaredNorm())),
379  circle_results.qCharge * (1. / fast_fit(2) + uVec(n));
380 
381  assert(circle_results.qCharge * circle_results.par(1) <= 0);
382 
383  riemannFit::Vector2d eMinusd = eVec - dVec;
384  double tmp1 = eMinusd.squaredNorm();
385  double tmp2 = sqrt(riemannFit::sqr(2 * fast_fit(2)) - tmp1);
386 
387  riemannFit::Matrix3d jacobian;
388  jacobian << (radii(1, 0) * eMinusd(0) - eMinusd(1) * radii(0, 0)) / tmp1,
389  (radii(1, 1) * eMinusd(0) - eMinusd(1) * radii(0, 1)) / tmp1, 0,
390  (circle_results.qCharge / 2) * (eMinusd(0) * radii(0, 0) + eMinusd(1) * radii(1, 0)) / tmp2,
391  (circle_results.qCharge / 2) * (eMinusd(0) * radii(0, 1) + eMinusd(1) * radii(1, 1)) / tmp2, 0, 0, 0,
392  circle_results.qCharge;
393 
394  circle_results.cov << iMat(0, 0), iMat(0, 1), iMat(0, n), iMat(1, 0), iMat(1, 1), iMat(1, n), iMat(n, 0),
395  iMat(n, 1), iMat(n, n);
396 
397  circle_results.cov = jacobian * circle_results.cov * jacobian.transpose();
398 
399  //...Translate in the system in which the first corrected hit is the origin, adding the m.s. correction...
400 
401  auto eMinusDVec = eVec - dVec;
402  translateKarimaki(circle_results, 0.5 * eMinusDVec(0), 0.5 * eMinusDVec(1), jacobian);
403  circle_results.cov(0, 0) +=
404  (1 + riemannFit::sqr(slope)) * multScatt(sTotal(1) - sTotal(0), bField, fast_fit(2), 2, slope);
405 
406  //...And translate back to the original system
407 
408  translateKarimaki(circle_results, dVec(0), dVec(1), jacobian);
409 
410  // compute chi2
411  circle_results.chi2 = 0;
412  for (u_int i = 0; i < n; i++) {
413  circle_results.chi2 += weightsVec(i) * riemannFit::sqr(zInSZplane(i) - uVec(i));
414  if (i > 0 && i < n - 1)
415  circle_results.chi2 +=
416  riemannFit::sqr(uVec(i - 1) / (sTransverse(i) - sTransverse(i - 1)) -
417  uVec(i) * (sTransverse(i + 1) - sTransverse(i - 1)) /
418  ((sTransverse(i + 1) - sTransverse(i)) * (sTransverse(i) - sTransverse(i - 1))) +
419  uVec(i + 1) / (sTransverse(i + 1) - sTransverse(i)) +
420  (sTransverse(i + 1) - sTransverse(i - 1)) * uVec(n) / 2) /
421  varBeta(i);
422  }
423 
424  // assert(circle_results.chi2>=0);
425  }
426 
449  template <typename V4, typename M6xN, int n>
450  __host__ __device__ inline void lineFit(const M6xN& hits_ge,
451  const V4& fast_fit,
452  const double bField,
453  const PreparedBrokenLineData<n>& data,
454  riemannFit::LineFit& line_results) {
455  const auto& radii = data.radii;
456  const auto& sTotal = data.sTotal;
457  const auto& zInSZplane = data.zInSZplane;
458  const auto& varBeta = data.varBeta;
459 
460  const double slope = -data.qCharge / fast_fit(3);
461  riemannFit::Matrix2d rotMat = rotationMatrix(slope);
462 
463  riemannFit::Matrix3d vMat = riemannFit::Matrix3d::Zero(); // covariance matrix XYZ
464  riemannFit::Matrix2x3d jacobXYZtosZ =
465  riemannFit::Matrix2x3d::Zero(); // jacobian for computation of the error on s (xyz -> sz)
466  riemannFit::VectorNd<n> weights = riemannFit::VectorNd<n>::Zero();
467  for (u_int i = 0; i < n; i++) {
468  vMat(0, 0) = hits_ge.col(i)[0]; // x errors
469  vMat(0, 1) = vMat(1, 0) = hits_ge.col(i)[1]; // cov_xy
470  vMat(0, 2) = vMat(2, 0) = hits_ge.col(i)[3]; // cov_xz
471  vMat(1, 1) = hits_ge.col(i)[2]; // y errors
472  vMat(2, 1) = vMat(1, 2) = hits_ge.col(i)[4]; // cov_yz
473  vMat(2, 2) = hits_ge.col(i)[5]; // z errors
474  auto tmp = 1. / radii.block(0, i, 2, 1).norm();
475  jacobXYZtosZ(0, 0) = radii(1, i) * tmp;
476  jacobXYZtosZ(0, 1) = -radii(0, i) * tmp;
477  jacobXYZtosZ(1, 2) = 1.;
478  weights(i) = 1. / ((rotMat * jacobXYZtosZ * vMat * jacobXYZtosZ.transpose() * rotMat.transpose())(
479  1, 1)); // compute the orthogonal weight point by point
480  }
481 
482  riemannFit::VectorNd<n> r_u;
483  for (u_int i = 0; i < n; i++) {
484  r_u(i) = weights(i) * zInSZplane(i);
485  }
486 #ifdef CPP_DUMP
487  std::cout << "CU4\n" << matrixC_u(w, sTotal, varBeta) << std::endl;
488 #endif
489  riemannFit::MatrixNd<n> iMat;
490  math::cholesky::invert(matrixC_u(weights, sTotal, varBeta), iMat);
491 #ifdef CPP_DUMP
492  std::cout << "I4\n" << iMat << std::endl;
493 #endif
494 
495  riemannFit::VectorNd<n> uVec = iMat * r_u; // obtain the fitted parameters by solving the linear system
496 
497  // line parameters in the system in which the first hit is the origin and with axis along SZ
498  line_results.par << (uVec(1) - uVec(0)) / (sTotal(1) - sTotal(0)), uVec(0);
499  auto idiff = 1. / (sTotal(1) - sTotal(0));
500  line_results.cov << (iMat(0, 0) - 2 * iMat(0, 1) + iMat(1, 1)) * riemannFit::sqr(idiff) +
501  multScatt(sTotal(1) - sTotal(0), bField, fast_fit(2), 2, slope),
502  (iMat(0, 1) - iMat(0, 0)) * idiff, (iMat(0, 1) - iMat(0, 0)) * idiff, iMat(0, 0);
503 
504  // translate to the original SZ system
505  riemannFit::Matrix2d jacobian;
506  jacobian(0, 0) = 1.;
507  jacobian(0, 1) = 0;
508  jacobian(1, 0) = -sTotal(0);
509  jacobian(1, 1) = 1.;
510  line_results.par(1) += -line_results.par(0) * sTotal(0);
511  line_results.cov = jacobian * line_results.cov * jacobian.transpose();
512 
513  // rotate to the original sz system
514  auto tmp = rotMat(0, 0) - line_results.par(0) * rotMat(0, 1);
515  jacobian(1, 1) = 1. / tmp;
516  jacobian(0, 0) = jacobian(1, 1) * jacobian(1, 1);
517  jacobian(0, 1) = 0;
518  jacobian(1, 0) = line_results.par(1) * rotMat(0, 1) * jacobian(0, 0);
519  line_results.par(1) = line_results.par(1) * jacobian(1, 1);
520  line_results.par(0) = (rotMat(0, 1) + line_results.par(0) * rotMat(0, 0)) * jacobian(1, 1);
521  line_results.cov = jacobian * line_results.cov * jacobian.transpose();
522 
523  // compute chi2
524  line_results.chi2 = 0;
525  for (u_int i = 0; i < n; i++) {
526  line_results.chi2 += weights(i) * riemannFit::sqr(zInSZplane(i) - uVec(i));
527  if (i > 0 && i < n - 1)
528  line_results.chi2 += riemannFit::sqr(uVec(i - 1) / (sTotal(i) - sTotal(i - 1)) -
529  uVec(i) * (sTotal(i + 1) - sTotal(i - 1)) /
530  ((sTotal(i + 1) - sTotal(i)) * (sTotal(i) - sTotal(i - 1))) +
531  uVec(i + 1) / (sTotal(i + 1) - sTotal(i))) /
532  varBeta(i);
533  }
534  }
535 
570  template <int n>
571  inline riemannFit::HelixFit helixFit(const riemannFit::Matrix3xNd<n>& hits,
572  const Eigen::Matrix<float, 6, 4>& hits_ge,
573  const double bField) {
574  riemannFit::HelixFit helix;
575  riemannFit::Vector4d fast_fit;
576  fastFit(hits, fast_fit);
577 
578  PreparedBrokenLineData<n> data;
579  karimaki_circle_fit circle;
580  riemannFit::LineFit line;
581  riemannFit::Matrix3d jacobian;
582 
583  prepareBrokenLineData(hits, fast_fit, bField, data);
584  lineFit(hits_ge, fast_fit, bField, data, line);
585  circleFit(hits, hits_ge, fast_fit, bField, data, circle);
586 
587  // the circle fit gives k, but here we want p_t, so let's change the parameter and the covariance matrix
588  jacobian << 1., 0, 0, 0, 1., 0, 0, 0,
589  -std::abs(circle.par(2)) * bField / (riemannFit::sqr(circle.par(2)) * circle.par(2));
590  circle.par(2) = bField / std::abs(circle.par(2));
591  circle.cov = jacobian * circle.cov * jacobian.transpose();
592 
593  helix.par << circle.par, line.par;
594  helix.cov = riemannFit::MatrixXd::Zero(5, 5);
595  helix.cov.block(0, 0, 3, 3) = circle.cov;
596  helix.cov.block(3, 3, 2, 2) = line.cov;
597  helix.qCharge = circle.qCharge;
598  helix.chi2_circle = circle.chi2;
599  helix.chi2_line = line.chi2;
600 
601  return helix;
602  }
603 
604 } // namespace brokenline
605 
606 #endif // RecoPixelVertexing_PixelTrackFitting_interface_BrokenLine_h
riemannFit::Vector4d
Eigen::Vector4d Vector4d
Definition: FitResult.h:15
riemannFit::LineFit::par
Vector2d par
(cotan(theta),Zip)
Definition: FitResult.h:39
riemannFit::sqr
constexpr T sqr(const T a)
raise to square.
Definition: FitUtils.h:67
riemannFit::CircleFit::par
Vector3d par
parameter: (X0,Y0,R)
Definition: FitResult.h:27
mps_fire.i
i
Definition: mps_fire.py:428
funct::Constant
Polynomial< 0 > Constant
Definition: Constant.h:6
dqmiodumpmetadata.n
n
Definition: dqmiodumpmetadata.py:28
riemannFit::LineFit::cov
Matrix2d cov
Definition: FitResult.h:40
riemannFit::LineFit::chi2
double chi2
Definition: FitResult.h:45
hfClusterShapes_cfi.hits
hits
Definition: hfClusterShapes_cfi.py:5
brokenline::prepareBrokenLineData
__host__ __device__ void prepareBrokenLineData(const M3xN &hits, const V4 &fast_fit, const double bField, PreparedBrokenLineData< n > &results)
Computes the data needed for the Broken Line fit procedure that are mainly common for the circle and ...
Definition: BrokenLine.h:150
brokenline::PreparedBrokenLineData::varBeta
riemannFit::VectorNd< n > varBeta
kink angles in the SZ plane
Definition: BrokenLine.h:31
amptDefaultParameters_cff.mu
mu
Definition: amptDefaultParameters_cff.py:16
min
T min(T a, T b)
Definition: MathUtil.h:58
riemannFit::HelixFit::qCharge
int32_t qCharge
particle charge
Definition: FitResult.h:61
gather_cfg.cout
cout
Definition: gather_cfg.py:144
riemannFit::VectorNplusONEd
Eigen::Matrix< double, N+1, 1 > VectorNplusONEd
Definition: FitUtils.h:35
brokenline::PreparedBrokenLineData::sTotal
riemannFit::VectorNd< n > sTotal
total distance traveled (three-dimensional)
Definition: BrokenLine.h:29
cms::cuda::assert
assert(be >=bs)
bookConverter.results
results
Definition: bookConverter.py:144
brokenline::matrixC_u
__host__ __device__ riemannFit::MatrixNd< n > matrixC_u(const riemannFit::VectorNd< n > &weights, const riemannFit::VectorNd< n > &sTotal, const riemannFit::VectorNd< n > &varBeta)
Computes the n-by-n band matrix obtained minimizing the Broken Line's cost function w....
Definition: BrokenLine.h:206
brokenline::circleFit
__host__ __device__ void circleFit(const M3xN &hits, const M6xN &hits_ge, const V4 &fast_fit, const double bField, PreparedBrokenLineData< n > &data, karimaki_circle_fit &circle_results)
Performs the Broken Line fit in the curved track case (that is, the fit parameters are the intercepti...
Definition: BrokenLine.h:296
createJobs.tmp
tmp
align.sh
Definition: createJobs.py:716
riemannFit::Matrix2d
Eigen::Matrix2d Matrix2d
Definition: FitResult.h:17
riemannFit::CircleFit
Definition: FitResult.h:26
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
brokenline
Definition: BrokenLine.h:8
riemannFit::HelixFit::cov
Matrix5d cov
Definition: FitResult.h:50
brokenline::rotationMatrix
__host__ __device__ riemannFit::Matrix2d rotationMatrix(double slope)
Computes the 2D rotation matrix that transforms the line y=slope*x into the line y=0.
Definition: BrokenLine.h:77
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
riemannFit::MatrixNd
Eigen::Matrix< double, N, N > MatrixNd
Definition: FitUtils.h:15
brokenline::lineFit
__host__ __device__ void lineFit(const M6xN &hits_ge, const V4 &fast_fit, const double bField, const PreparedBrokenLineData< n > &data, riemannFit::LineFit &line_results)
Performs the Broken Line fit in the straight track case (that is, the fit parameters are only the int...
Definition: BrokenLine.h:450
riemannFit::CircleFit::qCharge
int32_t qCharge
particle charge
Definition: FitResult.h:34
w
const double w
Definition: UKUtility.cc:23
HLT_FULL_cff.weights
weights
Definition: HLT_FULL_cff.py:99166
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
brokenline::PreparedBrokenLineData::sTransverse
riemannFit::VectorNd< n > sTransverse
total distance traveled in the transverse plane
Definition: BrokenLine.h:27
hitfit::scalar
double scalar(const CLHEP::HepGenMatrix &m)
Return the matrix as a scalar. Raise an assertion if the matris is not .
Definition: matutil.cc:166
riemannFit::Matrix2xNd
Eigen::Matrix< double, 2, N > Matrix2xNd
Definition: FitUtils.h:25
riemannFit::HelixFit::chi2_line
float chi2_line
Definition: FitResult.h:59
b
double b
Definition: hdecay.h:118
phase1PixelTopology::layer
constexpr std::array< uint8_t, layerIndexSize > layer
Definition: phase1PixelTopology.h:99
riemannFit::MatrixNplusONEd
Eigen::Matrix< double, N+1, N+1 > MatrixNplusONEd
Definition: FitUtils.h:17
riemannFit::CircleFit::chi2
float chi2
Definition: FitResult.h:35
a
double a
Definition: hdecay.h:119
riemannFit::Matrix3xNd
Eigen::Matrix< double, 3, N > Matrix3xNd
Definition: FitResult.h:24
type
type
Definition: SiPixelVCal_PayloadInspector.cc:37
riemannFit::HelixFit::chi2_circle
float chi2_circle
Definition: FitResult.h:58
fact
const double fact
Definition: NuclearInteractionFTFSimulator.cc:74
riemannFit::CircleFit::cov
Matrix3d cov
Definition: FitResult.h:28
riemannFit::LineFit
Definition: FitResult.h:38
brokenline::PreparedBrokenLineData::qCharge
int qCharge
particle charge
Definition: BrokenLine.h:25
__device__
#define __device__
Definition: SiPixelGainForHLTonGPU.h:15
HLT_FULL_cff.pt2
pt2
Definition: HLT_FULL_cff.py:9872
riemannFit::HelixFit::par
Vector5d par
(phi,Tip,pt,cotan(theta)),Zip)
Definition: FitResult.h:49
riemannFit::VectorNd
Eigen::Matrix< double, N, 1 > VectorNd
Definition: FitUtils.h:33
brokenline::PreparedBrokenLineData
data needed for the Broken Line fit procedure.
Definition: BrokenLine.h:24
brokenline::PreparedBrokenLineData::zInSZplane
riemannFit::VectorNd< n > zInSZplane
orthogonal coordinate to the pre-fitted line in the sz plane
Definition: BrokenLine.h:30
protons_cff.xi
xi
Definition: protons_cff.py:35
brokenline::translateKarimaki
__host__ __device__ void translateKarimaki(karimaki_circle_fit &circle, double x0, double y0, riemannFit::Matrix3d &jacobian)
Changes the Karimäki parameters (and consequently their covariance matrix) under a translation of the...
Definition: BrokenLine.h:97
Calorimetry_cff.bField
bField
Definition: Calorimetry_cff.py:284
l1tstage2_dqm_sourceclient-live_cfg.invert
invert
Definition: l1tstage2_dqm_sourceclient-live_cfg.py:85
riemannFit::HelixFit
Definition: FitResult.h:48
riemannFit::Vector2d
Eigen::Vector2d Vector2d
Definition: FitResult.h:13
makeMuonMisalignmentScenario.rot
rot
Definition: makeMuonMisalignmentScenario.py:322
CosmicsPD_Skims.radius
radius
Definition: CosmicsPD_Skims.py:135
dqm-mbProfile.log
log
Definition: dqm-mbProfile.py:17
data
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:79
brokenline::multScatt
__host__ __device__ double multScatt(const double &length, const double bField, const double radius, int layer, double slope)
Computes the Coulomb multiple scattering variance of the planar angle.
Definition: BrokenLine.h:53
ztail.d
d
Definition: ztail.py:151
riemannFit::Matrix2x3d
Eigen::Matrix< double, 2, 3 > Matrix2x3d
Definition: FitUtils.h:45
mps_fire.result
result
Definition: mps_fire.py:311
FitUtils.h
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
brokenline::helixFit
riemannFit::HelixFit helixFit(const riemannFit::Matrix3xNd< n > &hits, const Eigen::Matrix< float, 6, 4 > &hits_ge, const double bField)
Helix fit by three step: -fast pre-fit (see Fast_fit() for further info); -circle fit of the hits pr...
Definition: BrokenLine.h:571
brokenline::fastFit
__host__ __device__ void fastFit(const M3xN &hits, V4 &result)
A very fast helix fit.
Definition: BrokenLine.h:249
riemannFit::cross2D
__host__ __device__ double cross2D(const Vector2d &a, const Vector2d &b)
Compute cross product of two 2D vector (assuming z component 0), returning z component of the result.
Definition: FitUtils.h:79
riemannFit::Matrix3d
Eigen::Matrix3d Matrix3d
Definition: FitResult.h:18
c
auto & c
Definition: CAHitNtupletGeneratorKernelsImpl.h:46
slope
static const double slope[3]
Definition: CastorTimeSlew.cc:6
brokenline::PreparedBrokenLineData::radii
riemannFit::Matrix2xNd< n > radii
xy data in the system in which the pre-fitted center is the origin
Definition: BrokenLine.h:26
mps_splice.line
line
Definition: mps_splice.py:76
__host__
#define __host__
Definition: SiPixelGainForHLTonGPU.h:12
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
Generated for CMSSW Reference Manual by doxygen 1.8.16