BrokenLine.h

Go to the documentation of this file.

#ifndef RecoPixelVertexing_PixelTrackFitting_interface_BrokenLine_h
#define RecoPixelVertexing_PixelTrackFitting_interface_BrokenLine_h
#include <alpaka/alpaka.hpp>
#include <Eigen/Eigenvalues>
#include "HeterogeneousCore/AlpakaInterface/interface/config.h"
#include "RecoTracker/PixelTrackFitting/interface/alpaka/FitUtils.h"

namespace ALPAKA_ACCELERATOR_NAMESPACE {
  namespace brokenline {
    using namespace cms::alpakatools;
    using namespace ::riemannFit;


    using karimaki_circle_fit = riemannFit::CircleFit;

    template <int n>
    struct PreparedBrokenLineData {
      int qCharge;                          
      riemannFit::Matrix2xNd<n> radii;      
      riemannFit::VectorNd<n> sTransverse;  
                                            //   starting from the pre-fitted closest approach
      riemannFit::VectorNd<n> sTotal;       
      riemannFit::VectorNd<n> zInSZplane;   
      riemannFit::VectorNd<n> varBeta;      
    };

    template <typename TAcc>
    ALPAKA_FN_ACC ALPAKA_FN_INLINE double multScatt(
        const TAcc& acc, const double& length, const double bField, const double radius, int layer, double slope) {
      // limit R to 20GeV...
      auto pt2 = alpaka::math::min(acc, 20., bField * radius);
      pt2 *= pt2;
      constexpr double inv_X0 = 0.06 / 16.;  
      //if(Layer==1) XXI_0=0.06/16.;
      // else XXI_0=0.06/16.;
      //XX_0*=1;

      constexpr double geometry_factor = 0.7;
      constexpr double fact = geometry_factor * riemannFit::sqr(13.6 / 1000.);
      return fact / (pt2 * (1. + riemannFit::sqr(slope))) * (alpaka::math::abs(acc, length) * inv_X0) *
             riemannFit::sqr(1. + 0.038 * log(alpaka::math::abs(acc, length) * inv_X0));
    }

    template <typename TAcc>
    ALPAKA_FN_ACC ALPAKA_FN_INLINE riemannFit::Matrix2d rotationMatrix(const TAcc& acc, double slope) {
      riemannFit::Matrix2d rot;
      rot(0, 0) = 1. / alpaka::math::sqrt(acc, 1. + riemannFit::sqr(slope));
      rot(0, 1) = slope * rot(0, 0);
      rot(1, 0) = -rot(0, 1);
      rot(1, 1) = rot(0, 0);
      return rot;
    }

    template <typename TAcc>
    ALPAKA_FN_ACC ALPAKA_FN_INLINE void translateKarimaki(
        const TAcc& acc, karimaki_circle_fit& circle, double x0, double y0, riemannFit::Matrix3d& jacobian) {
      // Avoid multiple access to the circle.par vector.
      using scalar = typename std::remove_reference<decltype(circle.par(0))>::type;
      scalar phi = circle.par(0);
      scalar dee = circle.par(1);
      scalar rho = circle.par(2);

      // Avoid repeated trig. computations
      scalar sinPhi = alpaka::math::sin(acc, phi);
      scalar cosPhi = alpaka::math::cos(acc, phi);

      // Intermediate computations for the circle parameters
      scalar deltaPara = x0 * cosPhi + y0 * sinPhi;
      scalar deltaOrth = x0 * sinPhi - y0 * cosPhi + dee;
      scalar tempSmallU = 1 + rho * dee;
      scalar tempC = -rho * y0 + tempSmallU * cosPhi;
      scalar tempB = rho * x0 + tempSmallU * sinPhi;
      scalar tempA = 2. * deltaOrth + rho * (riemannFit::sqr(deltaOrth) + riemannFit::sqr(deltaPara));
      scalar tempU = alpaka::math::sqrt(acc, 1. + rho * tempA);

      // Intermediate computations for the error matrix transform
      scalar xi = 1. / (riemannFit::sqr(tempB) + riemannFit::sqr(tempC));
      scalar tempV = 1. + rho * deltaOrth;
      scalar lambda = (0.5 * tempA) / (riemannFit::sqr(1. + tempU) * tempU);
      scalar mu = 1. / (tempU * (1. + tempU)) + rho * lambda;
      scalar zeta = riemannFit::sqr(deltaOrth) + riemannFit::sqr(deltaPara);
      jacobian << xi * tempSmallU * tempV, -xi * riemannFit::sqr(rho) * deltaOrth, xi * deltaPara,
          2. * mu * tempSmallU * deltaPara, 2. * mu * tempV, mu * zeta - lambda * tempA, 0, 0, 1.;

      // translated circle parameters
      // phi
      circle.par(0) = alpaka::math::atan2(acc, tempB, tempC);
      // d
      circle.par(1) = tempA / (1 + tempU);
      // rho after translation. It is invariant, so noop
      // circle.par(2)= rho;

      // translated error matrix
      circle.cov = jacobian * circle.cov * jacobian.transpose();
    }

    template <typename TAcc, typename M3xN, typename V4, int n>
    ALPAKA_FN_ACC ALPAKA_FN_INLINE void __attribute__((always_inline))
    prepareBrokenLineData(const TAcc& acc,
                          const M3xN& hits,
                          const V4& fast_fit,
                          const double bField,
                          PreparedBrokenLineData<n>& results) {
      riemannFit::Vector2d dVec;
      riemannFit::Vector2d eVec;

      int mId = 1;

      if constexpr (n > 3) {
        riemannFit::Vector2d middle = 0.5 * (hits.block(0, n - 1, 2, 1) + hits.block(0, 0, 2, 1));
        auto d1 = (hits.block(0, n / 2, 2, 1) - middle).squaredNorm();
        auto d2 = (hits.block(0, n / 2 - 1, 2, 1) - middle).squaredNorm();
        mId = d1 < d2 ? n / 2 : n / 2 - 1;
      }

      dVec = hits.block(0, mId, 2, 1) - hits.block(0, 0, 2, 1);
      eVec = hits.block(0, n - 1, 2, 1) - hits.block(0, mId, 2, 1);
      results.qCharge = riemannFit::cross2D(acc, dVec, eVec) > 0 ? -1 : 1;

      const double slope = -results.qCharge / fast_fit(3);

      riemannFit::Matrix2d rotMat = rotationMatrix(acc, slope);

      // calculate radii and s
      results.radii = hits.block(0, 0, 2, n) - fast_fit.head(2) * riemannFit::MatrixXd::Constant(1, n, 1);
      eVec = -fast_fit(2) * fast_fit.head(2) / fast_fit.head(2).norm();
      for (u_int i = 0; i < n; i++) {
        dVec = results.radii.block(0, i, 2, 1);
        results.sTransverse(i) =
            results.qCharge * fast_fit(2) *
            alpaka::math::atan2(
                acc, riemannFit::cross2D(acc, dVec, eVec), dVec.dot(eVec));  // calculates the arc length
      }
      riemannFit::VectorNd<n> zVec = hits.block(2, 0, 1, n).transpose();

      //calculate sTotal and zVec
      riemannFit::Matrix2xNd<n> pointsSZ = riemannFit::Matrix2xNd<n>::Zero();
      for (u_int i = 0; i < n; i++) {
        pointsSZ(0, i) = results.sTransverse(i);
        pointsSZ(1, i) = zVec(i);
        pointsSZ.block(0, i, 2, 1) = rotMat * pointsSZ.block(0, i, 2, 1);
      }
      results.sTotal = pointsSZ.block(0, 0, 1, n).transpose();
      results.zInSZplane = pointsSZ.block(1, 0, 1, n).transpose();

      //calculate varBeta
      results.varBeta(0) = results.varBeta(n - 1) = 0;
      for (u_int i = 1; i < n - 1; i++) {
        results.varBeta(i) =
            multScatt(acc, results.sTotal(i + 1) - results.sTotal(i), bField, fast_fit(2), i + 2, slope) +
            multScatt(acc, results.sTotal(i) - results.sTotal(i - 1), bField, fast_fit(2), i + 1, slope);
      }
    }

    template <typename TAcc, int n>
    ALPAKA_FN_ACC ALPAKA_FN_INLINE riemannFit::MatrixNd<n> matrixC_u(const TAcc& acc,
                                                                     const riemannFit::VectorNd<n>& weights,
                                                                     const riemannFit::VectorNd<n>& sTotal,
                                                                     const riemannFit::VectorNd<n>& varBeta) {
      riemannFit::MatrixNd<n> c_uMat = riemannFit::MatrixNd<n>::Zero();
      for (u_int i = 0; i < n; i++) {
        c_uMat(i, i) = weights(i);
        if (i > 1)
          c_uMat(i, i) += 1. / (varBeta(i - 1) * riemannFit::sqr(sTotal(i) - sTotal(i - 1)));
        if (i > 0 && i < n - 1)
          c_uMat(i, i) +=
              (1. / varBeta(i)) * riemannFit::sqr((sTotal(i + 1) - sTotal(i - 1)) /
                                                  ((sTotal(i + 1) - sTotal(i)) * (sTotal(i) - sTotal(i - 1))));
        if (i < n - 2)
          c_uMat(i, i) += 1. / (varBeta(i + 1) * riemannFit::sqr(sTotal(i + 1) - sTotal(i)));

        if (i > 0 && i < n - 1)
          c_uMat(i, i + 1) =
              1. / (varBeta(i) * (sTotal(i + 1) - sTotal(i))) *
              (-(sTotal(i + 1) - sTotal(i - 1)) / ((sTotal(i + 1) - sTotal(i)) * (sTotal(i) - sTotal(i - 1))));
        if (i < n - 2)
          c_uMat(i, i + 1) +=
              1. / (varBeta(i + 1) * (sTotal(i + 1) - sTotal(i))) *
              (-(sTotal(i + 2) - sTotal(i)) / ((sTotal(i + 2) - sTotal(i + 1)) * (sTotal(i + 1) - sTotal(i))));

        if (i < n - 2)
          c_uMat(i, i + 2) = 1. / (varBeta(i + 1) * (sTotal(i + 2) - sTotal(i + 1)) * (sTotal(i + 1) - sTotal(i)));

        c_uMat(i, i) *= 0.5;
      }
      return c_uMat + c_uMat.transpose();
    }

    template <typename TAcc, typename M3xN, typename V4>
    ALPAKA_FN_ACC ALPAKA_FN_INLINE void fastFit(const TAcc& acc, const M3xN& hits, V4& result) {
      constexpr uint32_t n = M3xN::ColsAtCompileTime;

      int mId = 1;

      if constexpr (n > 3) {
        riemannFit::Vector2d middle = 0.5 * (hits.block(0, n - 1, 2, 1) + hits.block(0, 0, 2, 1));
        auto d1 = (hits.block(0, n / 2, 2, 1) - middle).squaredNorm();
        auto d2 = (hits.block(0, n / 2 - 1, 2, 1) - middle).squaredNorm();
        mId = d1 < d2 ? n / 2 : n / 2 - 1;
      }

      const riemannFit::Vector2d a = hits.block(0, mId, 2, 1) - hits.block(0, 0, 2, 1);
      const riemannFit::Vector2d b = hits.block(0, n - 1, 2, 1) - hits.block(0, mId, 2, 1);
      const riemannFit::Vector2d c = hits.block(0, 0, 2, 1) - hits.block(0, n - 1, 2, 1);

      auto tmp = 0.5 / riemannFit::cross2D(acc, c, a);
      result(0) = hits(0, 0) - (a(1) * c.squaredNorm() + c(1) * a.squaredNorm()) * tmp;
      result(1) = hits(1, 0) + (a(0) * c.squaredNorm() + c(0) * a.squaredNorm()) * tmp;
      // check Wikipedia for these formulas

      result(2) = alpaka::math::sqrt(acc, a.squaredNorm() * b.squaredNorm() * c.squaredNorm()) /
                  (2. * alpaka::math::abs(acc, riemannFit::cross2D(acc, b, a)));
      // Using Math Olympiad's formula R=abc/(4A)

      const riemannFit::Vector2d d = hits.block(0, 0, 2, 1) - result.head(2);
      const riemannFit::Vector2d e = hits.block(0, n - 1, 2, 1) - result.head(2);

      result(3) = result(2) * atan2(riemannFit::cross2D(acc, d, e), d.dot(e)) / (hits(2, n - 1) - hits(2, 0));
      // ds/dz slope between last and first point
    }

    template <typename TAcc, typename M3xN, typename M6xN, typename V4, int n>
    ALPAKA_FN_ACC ALPAKA_FN_INLINE void circleFit(const TAcc& acc,
                                                  const M3xN& hits,
                                                  const M6xN& hits_ge,
                                                  const V4& fast_fit,
                                                  const double bField,
                                                  PreparedBrokenLineData<n>& data,
                                                  karimaki_circle_fit& circle_results) {
      circle_results.qCharge = data.qCharge;
      auto& radii = data.radii;
      const auto& sTransverse = data.sTransverse;
      const auto& sTotal = data.sTotal;
      auto& zInSZplane = data.zInSZplane;
      auto& varBeta = data.varBeta;
      const double slope = -circle_results.qCharge / fast_fit(3);
      varBeta *= 1. + riemannFit::sqr(slope);  // the kink angles are projected!

      for (u_int i = 0; i < n; i++) {
        zInSZplane(i) = radii.block(0, i, 2, 1).norm() - fast_fit(2);
      }

      riemannFit::Matrix2d vMat;           // covariance matrix
      riemannFit::VectorNd<n> weightsVec;  // weights
      riemannFit::Matrix2d rotMat;         // rotation matrix point by point
      for (u_int i = 0; i < n; i++) {
        vMat(0, 0) = hits_ge.col(i)[0];               // x errors
        vMat(0, 1) = vMat(1, 0) = hits_ge.col(i)[1];  // cov_xy
        vMat(1, 1) = hits_ge.col(i)[2];               // y errors
        rotMat = rotationMatrix(acc, -radii(0, i) / radii(1, i));
        weightsVec(i) =
            1. / ((rotMat * vMat * rotMat.transpose())(1, 1));  // compute the orthogonal weight point by point
      }

      riemannFit::VectorNplusONEd<n> r_uVec;
      r_uVec(n) = 0;
      for (u_int i = 0; i < n; i++) {
        r_uVec(i) = weightsVec(i) * zInSZplane(i);
      }

      riemannFit::MatrixNplusONEd<n> c_uMat;
      c_uMat.block(0, 0, n, n) = matrixC_u(acc, weightsVec, sTransverse, varBeta);
      c_uMat(n, n) = 0;
      //add the border to the c_uMat matrix
      for (u_int i = 0; i < n; i++) {
        c_uMat(i, n) = 0;
        if (i > 0 && i < n - 1) {
          c_uMat(i, n) +=
              -(sTransverse(i + 1) - sTransverse(i - 1)) * (sTransverse(i + 1) - sTransverse(i - 1)) /
              (2. * varBeta(i) * (sTransverse(i + 1) - sTransverse(i)) * (sTransverse(i) - sTransverse(i - 1)));
        }
        if (i > 1) {
          c_uMat(i, n) +=
              (sTransverse(i) - sTransverse(i - 2)) / (2. * varBeta(i - 1) * (sTransverse(i) - sTransverse(i - 1)));
        }
        if (i < n - 2) {
          c_uMat(i, n) +=
              (sTransverse(i + 2) - sTransverse(i)) / (2. * varBeta(i + 1) * (sTransverse(i + 1) - sTransverse(i)));
        }
        c_uMat(n, i) = c_uMat(i, n);
        if (i > 0 && i < n - 1)
          c_uMat(n, n) += riemannFit::sqr(sTransverse(i + 1) - sTransverse(i - 1)) / (4. * varBeta(i));
      }

#ifdef CPP_DUMP
      std::cout << "CU5\n" << c_uMat << std::endl;
#endif
      riemannFit::MatrixNplusONEd<n> iMat;
      math::cholesky::invert(c_uMat, iMat);
#ifdef CPP_DUMP
      std::cout << "I5\n" << iMat << std::endl;
#endif
      riemannFit::VectorNplusONEd<n> uVec = iMat * r_uVec;  // obtain the fitted parameters by solving the linear system

      // compute (phi, d_ca, k) in the system in which the midpoint of the first two corrected hits is the origin...

      radii.block(0, 0, 2, 1) /= radii.block(0, 0, 2, 1).norm();
      radii.block(0, 1, 2, 1) /= radii.block(0, 1, 2, 1).norm();

      riemannFit::Vector2d dVec = hits.block(0, 0, 2, 1) + (-zInSZplane(0) + uVec(0)) * radii.block(0, 0, 2, 1);
      riemannFit::Vector2d eVec = hits.block(0, 1, 2, 1) + (-zInSZplane(1) + uVec(1)) * radii.block(0, 1, 2, 1);
      auto eMinusd = eVec - dVec;
      auto eMinusd2 = eMinusd.squaredNorm();
      auto tmp1 = 1. / eMinusd2;
      auto tmp2 = alpaka::math::sqrt(acc, riemannFit::sqr(fast_fit(2)) - 0.25 * eMinusd2);

      circle_results.par << atan2(eMinusd(1), eMinusd(0)), circle_results.qCharge * (tmp2 - fast_fit(2)),
          circle_results.qCharge * (1. / fast_fit(2) + uVec(n));

      tmp2 = 1. / tmp2;

      riemannFit::Matrix3d jacobian;
      jacobian << (radii(1, 0) * eMinusd(0) - eMinusd(1) * radii(0, 0)) * tmp1,
          (radii(1, 1) * eMinusd(0) - eMinusd(1) * radii(0, 1)) * tmp1, 0,
          circle_results.qCharge * (eMinusd(0) * radii(0, 0) + eMinusd(1) * radii(1, 0)) * tmp2,
          circle_results.qCharge * (eMinusd(0) * radii(0, 1) + eMinusd(1) * radii(1, 1)) * tmp2, 0, 0, 0,
          circle_results.qCharge;

      circle_results.cov << iMat(0, 0), iMat(0, 1), iMat(0, n), iMat(1, 0), iMat(1, 1), iMat(1, n), iMat(n, 0),
          iMat(n, 1), iMat(n, n);

      circle_results.cov = jacobian * circle_results.cov * jacobian.transpose();

      //...Translate in the system in which the first corrected hit is the origin, adding the m.s. correction...

      translateKarimaki(acc, circle_results, 0.5 * eMinusd(0), 0.5 * eMinusd(1), jacobian);
      circle_results.cov(0, 0) +=
          (1 + riemannFit::sqr(slope)) * multScatt(acc, sTotal(1) - sTotal(0), bField, fast_fit(2), 2, slope);

      //...And translate back to the original system

      translateKarimaki(acc, circle_results, dVec(0), dVec(1), jacobian);

      // compute chi2
      circle_results.chi2 = 0;
      for (u_int i = 0; i < n; i++) {
        circle_results.chi2 += weightsVec(i) * riemannFit::sqr(zInSZplane(i) - uVec(i));
        if (i > 0 && i < n - 1)
          circle_results.chi2 +=
              riemannFit::sqr(uVec(i - 1) / (sTransverse(i) - sTransverse(i - 1)) -
                              uVec(i) * (sTransverse(i + 1) - sTransverse(i - 1)) /
                                  ((sTransverse(i + 1) - sTransverse(i)) * (sTransverse(i) - sTransverse(i - 1))) +
                              uVec(i + 1) / (sTransverse(i + 1) - sTransverse(i)) +
                              (sTransverse(i + 1) - sTransverse(i - 1)) * uVec(n) / 2) /
              varBeta(i);
      }
    }

    template <typename TAcc, typename V4, typename M6xN, int n>
    ALPAKA_FN_ACC ALPAKA_FN_INLINE void lineFit(const TAcc& acc,
                                                const M6xN& hits_ge,
                                                const V4& fast_fit,
                                                const double bField,
                                                const PreparedBrokenLineData<n>& data,
                                                riemannFit::LineFit& line_results) {
      const auto& radii = data.radii;
      const auto& sTotal = data.sTotal;
      const auto& zInSZplane = data.zInSZplane;
      const auto& varBeta = data.varBeta;

      const double slope = -data.qCharge / fast_fit(3);
      riemannFit::Matrix2d rotMat = rotationMatrix(acc, slope);

      riemannFit::Matrix3d vMat = riemannFit::Matrix3d::Zero();  // covariance matrix XYZ
      riemannFit::Matrix2x3d jacobXYZtosZ =
          riemannFit::Matrix2x3d::Zero();  // jacobian for computation of the error on s (xyz -> sz)
      riemannFit::VectorNd<n> weights = riemannFit::VectorNd<n>::Zero();
      for (u_int i = 0; i < n; i++) {
        vMat(0, 0) = hits_ge.col(i)[0];               // x errors
        vMat(0, 1) = vMat(1, 0) = hits_ge.col(i)[1];  // cov_xy
        vMat(0, 2) = vMat(2, 0) = hits_ge.col(i)[3];  // cov_xz
        vMat(1, 1) = hits_ge.col(i)[2];               // y errors
        vMat(2, 1) = vMat(1, 2) = hits_ge.col(i)[4];  // cov_yz
        vMat(2, 2) = hits_ge.col(i)[5];               // z errors
        auto tmp = 1. / radii.block(0, i, 2, 1).norm();
        jacobXYZtosZ(0, 0) = radii(1, i) * tmp;
        jacobXYZtosZ(0, 1) = -radii(0, i) * tmp;
        jacobXYZtosZ(1, 2) = 1.;
        weights(i) = 1. / ((rotMat * jacobXYZtosZ * vMat * jacobXYZtosZ.transpose() * rotMat.transpose())(
                              1, 1));  // compute the orthogonal weight point by point
      }

      riemannFit::VectorNd<n> r_u;
      for (u_int i = 0; i < n; i++) {
        r_u(i) = weights(i) * zInSZplane(i);
      }
#ifdef CPP_DUMP
      std::cout << "CU4\n" << matrixC_u(w, sTotal, varBeta) << std::endl;
#endif
      riemannFit::MatrixNd<n> iMat;
      math::cholesky::invert(matrixC_u(acc, weights, sTotal, varBeta), iMat);
#ifdef CPP_DUMP
      std::cout << "I4\n" << iMat << std::endl;
#endif

      riemannFit::VectorNd<n> uVec = iMat * r_u;  // obtain the fitted parameters by solving the linear system

      // line parameters in the system in which the first hit is the origin and with axis along SZ
      line_results.par << (uVec(1) - uVec(0)) / (sTotal(1) - sTotal(0)), uVec(0);
      auto idiff = 1. / (sTotal(1) - sTotal(0));
      line_results.cov << (iMat(0, 0) - 2 * iMat(0, 1) + iMat(1, 1)) * riemannFit::sqr(idiff) +
                              multScatt(acc, sTotal(1) - sTotal(0), bField, fast_fit(2), 2, slope),
          (iMat(0, 1) - iMat(0, 0)) * idiff, (iMat(0, 1) - iMat(0, 0)) * idiff, iMat(0, 0);

      // translate to the original SZ system
      riemannFit::Matrix2d jacobian;
      jacobian(0, 0) = 1.;
      jacobian(0, 1) = 0;
      jacobian(1, 0) = -sTotal(0);
      jacobian(1, 1) = 1.;
      line_results.par(1) += -line_results.par(0) * sTotal(0);
      line_results.cov = jacobian * line_results.cov * jacobian.transpose();

      // rotate to the original sz system
      auto tmp = rotMat(0, 0) - line_results.par(0) * rotMat(0, 1);
      jacobian(1, 1) = 1. / tmp;
      jacobian(0, 0) = jacobian(1, 1) * jacobian(1, 1);
      jacobian(0, 1) = 0;
      jacobian(1, 0) = line_results.par(1) * rotMat(0, 1) * jacobian(0, 0);
      line_results.par(1) = line_results.par(1) * jacobian(1, 1);
      line_results.par(0) = (rotMat(0, 1) + line_results.par(0) * rotMat(0, 0)) * jacobian(1, 1);
      line_results.cov = jacobian * line_results.cov * jacobian.transpose();

      // compute chi2
      line_results.chi2 = 0;
      for (u_int i = 0; i < n; i++) {
        line_results.chi2 += weights(i) * riemannFit::sqr(zInSZplane(i) - uVec(i));
        if (i > 0 && i < n - 1)
          line_results.chi2 += riemannFit::sqr(uVec(i - 1) / (sTotal(i) - sTotal(i - 1)) -
                                               uVec(i) * (sTotal(i + 1) - sTotal(i - 1)) /
                                                   ((sTotal(i + 1) - sTotal(i)) * (sTotal(i) - sTotal(i - 1))) +
                                               uVec(i + 1) / (sTotal(i + 1) - sTotal(i))) /
                               varBeta(i);
      }
    }

    template <int n>
    class helixFit {
    public:
      template <typename TAcc>
      ALPAKA_FN_ACC ALPAKA_FN_INLINE void operator()(const TAcc& acc,
                                                     const riemannFit::Matrix3xNd<n>* hits,
                                                     const Eigen::Matrix<float, 6, 4>* hits_ge,
                                                     const double bField,
                                                     riemannFit::HelixFit* helix) const {
        riemannFit::Vector4d fast_fit;
        fastFit(acc, *hits, fast_fit);

        PreparedBrokenLineData<n> data;
        karimaki_circle_fit circle;
        riemannFit::LineFit line;
        riemannFit::Matrix3d jacobian;

        prepareBrokenLineData(acc, *hits, fast_fit, bField, data);
        lineFit(acc, *hits_ge, fast_fit, bField, data, line);
        circleFit(acc, *hits, *hits_ge, fast_fit, bField, data, circle);

        // the circle fit gives k, but here we want p_t, so let's change the parameter and the covariance matrix
        jacobian << 1., 0, 0, 0, 1., 0, 0, 0,
            -alpaka::math::abs(acc, circle.par(2)) * bField / (riemannFit::sqr(circle.par(2)) * circle.par(2));
        circle.par(2) = bField / alpaka::math::abs(acc, circle.par(2));
        circle.cov = jacobian * circle.cov * jacobian.transpose();

        helix->par << circle.par, line.par;
        helix->cov = riemannFit::MatrixXd::Zero(5, 5);
        helix->cov.block(0, 0, 3, 3) = circle.cov;
        helix->cov.block(3, 3, 2, 2) = line.cov;
        helix->qCharge = circle.qCharge;
        helix->chi2_circle = circle.chi2;
        helix->chi2_line = line.chi2;
      }
    };
  }  // namespace brokenline
}  // namespace ALPAKA_ACCELERATOR_NAMESPACE
#endif  // RecoPixelVertexing_PixelTrackFitting_interface_BrokenLine_h