ALPAKA_ACCELERATOR_NAMESPACE::riemannFit Namespace Reference

Functions
template<typename TAcc , typename M2xN >
ALPAKA_FN_ACC ALPAKA_FN_INLINE int32_t	charge (const TAcc &acc, const M2xN &p2D, const Vector3d &par_uvr)
	Find particle q considering the sign of cross product between particles velocity (estimated by the first 2 hits) and the vector radius between the first hit and the center of the fitted circle. More...
template<typename TAcc , typename M2xN , typename V4 , int N>
ALPAKA_FN_ACC ALPAKA_FN_INLINE CircleFit	circleFit (const TAcc &acc, const M2xN &hits2D, const Matrix2Nd< N > &hits_cov2D, const V4 &fast_fit, const VectorNd< N > &rad, const double bField, const bool error)
	Fit a generic number of 2D points with a circle using Riemann-Chernov algorithm. Covariance matrix of fitted parameter is optionally computed. Multiple scattering (currently only in barrel layer) is optionally handled. More...
template<typename TAcc , typename VNd1 , typename VNd2 >
ALPAKA_FN_ACC ALPAKA_FN_INLINE void	computeRadLenUniformMaterial (const TAcc &acc, const VNd1 &length_values, VNd2 &rad_lengths)
template<typename TAcc , typename M2xN , int N>
ALPAKA_FN_ACC ALPAKA_FN_INLINE VectorNd< N >	cov_carttorad (const TAcc &acc, const M2xN &p2D, const Matrix2Nd< N > &cov_cart, const VectorNd< N > &rad)
	Transform covariance matrix from Cartesian coordinates (only transverse plane component) to radial coordinates (both radial and tangential component but only diagonal terms, correlation between different point are not managed). More...
template<typename TAcc , typename M2xN , typename V4 , int N>
ALPAKA_FN_ACC ALPAKA_FN_INLINE VectorNd< N >	cov_carttorad_prefit (const TAcc &acc, const M2xN &p2D, const Matrix2Nd< N > &cov_cart, V4 &fast_fit, const VectorNd< N > &rad)
	Transform covariance matrix from Cartesian coordinates (only transverse plane component) to coordinates system orthogonal to the pre-fitted circle in each point. Further information in attached documentation. More...
template<typename TAcc , typename M2xN , int N>
ALPAKA_FN_ACC ALPAKA_FN_INLINE Matrix2Nd< N >	cov_radtocart (const TAcc &acc, const M2xN &p2D, const MatrixNd< N > &cov_rad, const VectorNd< N > &rad)
	Transform covariance matrix from radial (only tangential component) to Cartesian coordinates (only transverse plane component). More...
template<typename TAcc >
ALPAKA_FN_ACC ALPAKA_FN_INLINE double	cross2D (const TAcc &acc, const Vector2d &a, const Vector2d &b)
	Compute cross product of two 2D vector (assuming z component 0), returning z component of the result. More...
template<typename TAcc , typename M3xN , typename V4 >
ALPAKA_FN_ACC ALPAKA_FN_INLINE void	fastFit (const TAcc &acc, const M3xN &hits, V4 &result)
	A very fast helix fit: it fits a circle by three points (first, middle and last point) and a line by two points (first and last). More...
template<typename TAcc >
ALPAKA_FN_ACC ALPAKA_FN_INLINE void	fromCircleToPerigee (const TAcc &acc, CircleFit &circle)
	Transform circle parameter from (X0,Y0,R) to (phi,Tip,q/R) and consequently covariance matrix. More...
template<typename TAcc , typename M3xN , typename M6xN , typename V4 >
ALPAKA_FN_ACC ALPAKA_FN_INLINE LineFit	lineFit (const TAcc &acc, const M3xN &hits, const M6xN &hits_ge, const CircleFit &circle, const V4 &fast_fit, const double bField, const bool error)
	Perform an ordinary least square fit in the s-z plane to compute the parameters cotTheta and Zip. More...
template<typename TAcc , typename M6xNf , typename M3xNd >
ALPAKA_FN_ACC void	loadCovariance (const TAcc &acc, M6xNf const &ge, M3xNd &hits_cov)
template<typename TAcc , typename M6xNf , typename M2Nd >
ALPAKA_FN_ACC void	loadCovariance2D (const TAcc &acc, M6xNf const &ge, M2Nd &hits_cov)
template<typename TAcc >
ALPAKA_FN_ACC ALPAKA_FN_INLINE Vector2d	min_eigen2D (const TAcc &acc, const Matrix2d &aMat, double &chi2)
	2D version of min_eigen3D(). More...
template<typename TAcc >
ALPAKA_FN_ACC ALPAKA_FN_INLINE Vector3d	min_eigen3D (const TAcc &acc, const Matrix3d &A, double &chi2)
	Compute the eigenvector associated to the minimum eigenvalue. More...
template<typename TAcc >
ALPAKA_FN_ACC ALPAKA_FN_INLINE Vector3d	min_eigen3D_fast (const TAcc &acc, const Matrix3d &A)
	A faster version of min_eigen3D() where double precision is not needed. More...
template<typename TAcc >
ALPAKA_FN_ACC ALPAKA_FN_INLINE void	par_uvrtopak (const TAcc &acc, CircleFit &circle, const double B, const bool error)
	Transform circle parameter from (X0,Y0,R) to (phi,Tip,p_t) and consequently covariance matrix. More...
template<typename TAcc , class C >
ALPAKA_FN_ACC void	printIt (const TAcc &acc, C *m, const char *prefix="")
template<typename TAcc , typename M2xN , typename V4 , int N>
ALPAKA_FN_ACC ALPAKA_FN_INLINE MatrixNd< N >	scatter_cov_rad (const TAcc &acc, const M2xN &p2D, const V4 &fast_fit, VectorNd< N > const &rad, double B)
	Compute the covariance matrix (in radial coordinates) of points in the transverse plane due to multiple Coulomb scattering. More...
template<typename TAcc , typename V4 , typename VNd1 , typename VNd2 , int N>
ALPAKA_FN_ACC ALPAKA_FN_INLINE auto	scatterCovLine (const TAcc &acc, Matrix2d const *cov_sz, const V4 &fast_fit, VNd1 const &s_arcs, VNd2 const &z_values, const double theta, const double bField, MatrixNd< N > &ret)
	Compute the covariance matrix along cartesian S-Z of points due to multiple Coulomb scattering to be used in the line_fit, for the barrel and forward cases. The input covariance matrix is in the variables s-z, original and unrotated. The multiple scattering component is computed in the usual linear approximation, using the 3D path which is computed as the squared root of the squared sum of the s and z components passed in. Internally a rotation by theta is performed and the covariance matrix returned is the one in the direction orthogonal to the rotated S3D axis, i.e. along the rotated Z axis. The choice of the rotation is not arbitrary, but derived from the fact that putting the horizontal axis along the S3D direction allows the usage of the ordinary least squared fitting techiques with the trivial parametrization y = mx + q, avoiding the patological case with m = +/- inf, that would correspond to the case at eta = 0. More...
template<typename T >
constexpr T	sqr (const T a)
	raise to square. More...
template<typename TAcc , int N>
ALPAKA_FN_ACC ALPAKA_FN_INLINE VectorNd< N >	weightCircle (const TAcc &acc, const MatrixNd< N > &cov_rad_inv)
	Compute the points' weights' vector for the circle fit when multiple scattering is managed. Further information in attached documentation. More...

Function Documentation

charge()

template<typename TAcc , typename M2xN >

ALPAKA_FN_ACC ALPAKA_FN_INLINE int32_t ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::charge	(	const TAcc &	acc,
		const M2xN &	p2D,
		const Vector3d &	par_uvr
	)

Find particle q considering the sign of cross product between particles velocity (estimated by the first 2 hits) and the vector radius between the first hit and the center of the fitted circle.

Parameters

p2D	2D points in transverse plane.
par_uvr	result of the circle fit in this form: (X0,Y0,R).

Returns

q int 1 or -1.

Definition at line 297 of file RiemannFit.h.

Referenced by circleFit().

                                                                                                           {
    return ((p2D(0, 1) - p2D(0, 0)) * (par_uvr.y() - p2D(1, 0)) - (p2D(1, 1) - p2D(1, 0)) * (par_uvr.x() - p2D(0, 0)) >
            0)
               ? -1
               : 1;
  }

circleFit()

template<typename TAcc , typename M2xN , typename V4 , int N>

ALPAKA_FN_ACC ALPAKA_FN_INLINE CircleFit ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::circleFit	(	const TAcc &	acc,
		const M2xN &	hits2D,
		const Matrix2Nd< N > &	hits_cov2D,
		const V4 &	fast_fit,
		const VectorNd< N > &	rad,
		const double	bField,
		const bool	error
	)

Fit a generic number of 2D points with a circle using Riemann-Chernov algorithm. Covariance matrix of fitted parameter is optionally computed. Multiple scattering (currently only in barrel layer) is optionally handled.

Parameters

hits2D	2D points to be fitted.
hits_cov2D	covariance matrix of 2D points.
fast_fit	pre-fit result in this form: (X0,Y0,R,tan(theta)). (tan(theta) is not used).
bField	magnetic field
error	flag for error computation.
scattering	flag for multiple scattering

Returns

circle circle_fit: -par parameter of the fitted circle in this form (X0,Y0,R);
-cov covariance matrix of the fitted parameter (not initialized if error = false);
-q charge of the particle;
-chi2.

Warning

hits must be passed ordered from inner to outer layer (double hits on the same layer must be ordered too) so that multiple scattering is treated properly.

Multiple scattering for barrel is still not tested.

Multiple scattering for endcap hits is not handled (yet). Do not fit endcap hits with scattering = true !

Definition at line 465 of file RiemannFit.h.

References a, funct::abs(), b, Calorimetry_cff::bField, charge(), nano_mu_local_reco_cff::chi2, ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), cov_carttorad_prefit(), cov_radtocart(), dumpMFGeometry_cfg::delta, geometryDiff::epsilon, relativeConstraints::error, ALPAKA_ACCELERATOR_NAMESPACE::brokenline::fast_fit, mps_fire::i, l1tstage2_dqm_sourceclient-live_cfg::invert, dqmiolumiharvest::j, dqmdumpme::k, MainPageGenerator::l, CaloTowersParam_cfi::mc, min_eigen3D(), min_eigen3D_fast(), N, dqmiodumpmetadata::n, printIt(), scatter_cov_rad(), Validation_hcalonly_cfi::sign, sqr(), mathSSE::sqrt(), submitPVValidationJobs::t, FrontierCondition_GT_autoExpress_cfi::t0, RandomServiceHelper::t1, createJobs::tmp, parallelization::uint, mps_merge::weight, and weightCircle().

Referenced by ALPAKA_ACCELERATOR_NAMESPACE::Kernel_CircleFit< N, TrackerTraits >::operator()(), and riemannFit::helixFit< N >::operator()().

                                                                       {
#ifdef RFIT_DEBUG
    printf("circle_fit - enter\n");
#endif
    // INITIALIZATION
    Matrix2Nd<N> vMat = hits_cov2D;
    constexpr uint n = N;
    printIt(&hits2D, "circle_fit - hits2D:");
    printIt(&hits_cov2D, "circle_fit - hits_cov2D:");

#ifdef RFIT_DEBUG
    printf("circle_fit - WEIGHT COMPUTATION\n");
#endif
    // WEIGHT COMPUTATION
    VectorNd<N> weight;
    MatrixNd<N> gMat;
    double renorm;
    {
      MatrixNd<N> cov_rad = cov_carttorad_prefit(acc, hits2D, vMat, fast_fit, rad).asDiagonal();
      MatrixNd<N> scatterCovRadMat = scatter_cov_rad(acc, hits2D, fast_fit, rad, bField);
      printIt(&scatterCovRadMat, "circle_fit - scatter_cov_rad:");
      printIt(&hits2D, "circle_fit - hits2D bis:");
#ifdef RFIT_DEBUG
      printf("Address of hits2D: a) %p\n", &hits2D);
#endif
      vMat += cov_radtocart(acc, hits2D, scatterCovRadMat, rad);
      printIt(&vMat, "circle_fit - V:");
      cov_rad += scatterCovRadMat;
      printIt(&cov_rad, "circle_fit - cov_rad:");
      math::cholesky::invert(cov_rad, gMat);
      // gMat = cov_rad.inverse();
      renorm = gMat.sum();
      gMat *= 1. / renorm;
      weight = weightCircle(acc, gMat);
    }
    printIt(&weight, "circle_fit - weight:");

    // SPACE TRANSFORMATION
#ifdef RFIT_DEBUG
    printf("circle_fit - SPACE TRANSFORMATION\n");
#endif

    // center
#ifdef RFIT_DEBUG
    printf("Address of hits2D: b) %p\n", &hits2D);
#endif
    const Vector2d hCentroid = hits2D.rowwise().mean();  // centroid
    printIt(&hCentroid, "circle_fit - h_:");
    Matrix3xNd<N> p3D;
    p3D.block(0, 0, 2, n) = hits2D.colwise() - hCentroid;
    printIt(&p3D, "circle_fit - p3D: a)");
    Vector2Nd<N> mc;  // centered hits, used in error computation
    mc << p3D.row(0).transpose(), p3D.row(1).transpose();
    printIt(&mc, "circle_fit - mc(centered hits):");

    // scale
    const double tempQ = mc.squaredNorm();
    const double tempS = sqrt(n * 1. / tempQ);  // scaling factor
    p3D.block(0, 0, 2, n) *= tempS;

    // project on paraboloid
    p3D.row(2) = p3D.block(0, 0, 2, n).colwise().squaredNorm();
    printIt(&p3D, "circle_fit - p3D: b)");

#ifdef RFIT_DEBUG
    printf("circle_fit - COST FUNCTION\n");
#endif
    // COST FUNCTION

    // compute
    Vector3d r0;
    r0.noalias() = p3D * weight;  // center of gravity
    const Matrix3xNd<N> xMat = p3D.colwise() - r0;
    Matrix3d aMat = xMat * gMat * xMat.transpose();
    printIt(&aMat, "circle_fit - A:");

#ifdef RFIT_DEBUG
    printf("circle_fit - MINIMIZE\n");
#endif
    // minimize
    double chi2;
    Vector3d vVec = min_eigen3D(acc, aMat, chi2);
#ifdef RFIT_DEBUG
    printf("circle_fit - AFTER MIN_EIGEN\n");
#endif
    printIt(&vVec, "v BEFORE INVERSION");
    vVec *= (vVec(2) > 0) ? 1 : -1;  // TO FIX dovrebbe essere N(3)>0
    printIt(&vVec, "v AFTER INVERSION");
    // This hack to be able to run on GPU where the automatic assignment to a
    // double from the vector multiplication is not working.
#ifdef RFIT_DEBUG
    printf("circle_fit - AFTER MIN_EIGEN 1\n");
#endif
    Eigen::Matrix<double, 1, 1> cm;
#ifdef RFIT_DEBUG
    printf("circle_fit - AFTER MIN_EIGEN 2\n");
#endif
    cm = -vVec.transpose() * r0;
#ifdef RFIT_DEBUG
    printf("circle_fit - AFTER MIN_EIGEN 3\n");
#endif
    const double tempC = cm(0, 0);

#ifdef RFIT_DEBUG
    printf("circle_fit - COMPUTE CIRCLE PARAMETER\n");
#endif
    // COMPUTE CIRCLE PARAMETER

    // auxiliary quantities
    const double tempH = sqrt(1. - sqr(vVec(2)) - 4. * tempC * vVec(2));
    const double v2x2_inv = 1. / (2. * vVec(2));
    const double s_inv = 1. / tempS;
    Vector3d par_uvr;  // used in error propagation
    par_uvr << -vVec(0) * v2x2_inv, -vVec(1) * v2x2_inv, tempH * v2x2_inv;

    CircleFit circle;
    circle.par << par_uvr(0) * s_inv + hCentroid(0), par_uvr(1) * s_inv + hCentroid(1), par_uvr(2) * s_inv;
    circle.qCharge = charge(acc, hits2D, circle.par);
    circle.chi2 = abs(chi2) * renorm / sqr(2 * vVec(2) * par_uvr(2) * tempS);
    printIt(&circle.par, "circle_fit - CIRCLE PARAMETERS:");
    printIt(&circle.cov, "circle_fit - CIRCLE COVARIANCE:");
#ifdef RFIT_DEBUG
    printf("circle_fit - CIRCLE CHARGE: %d\n", circle.qCharge);
#endif

#ifdef RFIT_DEBUG
    printf("circle_fit - ERROR PROPAGATION\n");
#endif
    // ERROR PROPAGATION
    if (error) {
#ifdef RFIT_DEBUG
      printf("circle_fit - ERROR PRPAGATION ACTIVATED\n");
#endif
      ArrayNd<N> vcsMat[2][2];  // cov matrix of center & scaled points
      MatrixNd<N> cMat[3][3];   // cov matrix of 3D transformed points
#ifdef RFIT_DEBUG
      printf("circle_fit - ERROR PRPAGATION ACTIVATED 2\n");
#endif
      {
        Eigen::Matrix<double, 1, 1> cm;
        Eigen::Matrix<double, 1, 1> cm2;
        cm = mc.transpose() * vMat * mc;
        const double tempC2 = cm(0, 0);
        Matrix2Nd<N> tempVcsMat;
        tempVcsMat.template triangularView<Eigen::Upper>() =
            (sqr(tempS) * vMat + sqr(sqr(tempS)) * 1. / (4. * tempQ * n) *
                                     (2. * vMat.squaredNorm() + 4. * tempC2) *  // mc.transpose() * V * mc) *
                                     (mc * mc.transpose()));

        printIt(&tempVcsMat, "circle_fit - Vcs:");
        cMat[0][0] = tempVcsMat.block(0, 0, n, n).template selfadjointView<Eigen::Upper>();
        vcsMat[0][1] = tempVcsMat.block(0, n, n, n);
        cMat[1][1] = tempVcsMat.block(n, n, n, n).template selfadjointView<Eigen::Upper>();
        vcsMat[1][0] = vcsMat[0][1].transpose();
        printIt(&tempVcsMat, "circle_fit - Vcs:");
      }

      {
        const ArrayNd<N> t0 = (VectorXd::Constant(n, 1.) * p3D.row(0));
        const ArrayNd<N> t1 = (VectorXd::Constant(n, 1.) * p3D.row(1));
        const ArrayNd<N> t00 = p3D.row(0).transpose() * p3D.row(0);
        const ArrayNd<N> t01 = p3D.row(0).transpose() * p3D.row(1);
        const ArrayNd<N> t11 = p3D.row(1).transpose() * p3D.row(1);
        const ArrayNd<N> t10 = t01.transpose();
        vcsMat[0][0] = cMat[0][0];
        cMat[0][1] = vcsMat[0][1];
        cMat[0][2] = 2. * (vcsMat[0][0] * t0 + vcsMat[0][1] * t1);
        vcsMat[1][1] = cMat[1][1];
        cMat[1][2] = 2. * (vcsMat[1][0] * t0 + vcsMat[1][1] * t1);
        MatrixNd<N> tmp;
        tmp.template triangularView<Eigen::Upper>() =
            (2. * (vcsMat[0][0] * vcsMat[0][0] + vcsMat[0][0] * vcsMat[0][1] + vcsMat[1][1] * vcsMat[1][0] +
                   vcsMat[1][1] * vcsMat[1][1]) +
             4. * (vcsMat[0][0] * t00 + vcsMat[0][1] * t01 + vcsMat[1][0] * t10 + vcsMat[1][1] * t11))
                .matrix();
        cMat[2][2] = tmp.template selfadjointView<Eigen::Upper>();
      }
      printIt(&cMat[0][0], "circle_fit - C[0][0]:");

      Matrix3d c0Mat;  // cov matrix of center of gravity (r0.x,r0.y,r0.z)
      for (uint i = 0; i < 3; ++i) {
        for (uint j = i; j < 3; ++j) {
          Eigen::Matrix<double, 1, 1> tmp;
          tmp = weight.transpose() * cMat[i][j] * weight;
          // Workaround to get things working in GPU
          const double tempC = tmp(0, 0);
          c0Mat(i, j) = tempC;  //weight.transpose() * C[i][j] * weight;
          c0Mat(j, i) = c0Mat(i, j);
        }
      }
      printIt(&c0Mat, "circle_fit - C0:");

      const MatrixNd<N> wMat = weight * weight.transpose();
      const MatrixNd<N> hMat = MatrixNd<N>::Identity().rowwise() - weight.transpose();
      const MatrixNx3d<N> s_v = hMat * p3D.transpose();
      printIt(&wMat, "circle_fit - W:");
      printIt(&hMat, "circle_fit - H:");
      printIt(&s_v, "circle_fit - s_v:");

      MatrixNd<N> dMat[3][3];  // cov(s_v)
      dMat[0][0] = (hMat * cMat[0][0] * hMat.transpose()).cwiseProduct(wMat);
      dMat[0][1] = (hMat * cMat[0][1] * hMat.transpose()).cwiseProduct(wMat);
      dMat[0][2] = (hMat * cMat[0][2] * hMat.transpose()).cwiseProduct(wMat);
      dMat[1][1] = (hMat * cMat[1][1] * hMat.transpose()).cwiseProduct(wMat);
      dMat[1][2] = (hMat * cMat[1][2] * hMat.transpose()).cwiseProduct(wMat);
      dMat[2][2] = (hMat * cMat[2][2] * hMat.transpose()).cwiseProduct(wMat);
      dMat[1][0] = dMat[0][1].transpose();
      dMat[2][0] = dMat[0][2].transpose();
      dMat[2][1] = dMat[1][2].transpose();
      printIt(&dMat[0][0], "circle_fit - D_[0][0]:");

      constexpr uint nu[6][2] = {{0, 0}, {0, 1}, {0, 2}, {1, 1}, {1, 2}, {2, 2}};

      Matrix6d eMat;  // cov matrix of the 6 independent elements of A
      for (uint a = 0; a < 6; ++a) {
        const uint i = nu[a][0], j = nu[a][1];
        for (uint b = a; b < 6; ++b) {
          const uint k = nu[b][0], l = nu[b][1];
          VectorNd<N> t0(n);
          VectorNd<N> t1(n);
          if (l == k) {
            t0 = 2. * dMat[j][l] * s_v.col(l);
            if (i == j)
              t1 = t0;
            else
              t1 = 2. * dMat[i][l] * s_v.col(l);
          } else {
            t0 = dMat[j][l] * s_v.col(k) + dMat[j][k] * s_v.col(l);
            if (i == j)
              t1 = t0;
            else
              t1 = dMat[i][l] * s_v.col(k) + dMat[i][k] * s_v.col(l);
          }

          if (i == j) {
            Eigen::Matrix<double, 1, 1> cm;
            cm = s_v.col(i).transpose() * (t0 + t1);
            // Workaround to get things working in GPU
            const double tempC = cm(0, 0);
            eMat(a, b) = 0. + tempC;
          } else {
            Eigen::Matrix<double, 1, 1> cm;
            cm = (s_v.col(i).transpose() * t0) + (s_v.col(j).transpose() * t1);
            // Workaround to get things working in GPU
            const double tempC = cm(0, 0);
            eMat(a, b) = 0. + tempC;  //(s_v.col(i).transpose() * t0) + (s_v.col(j).transpose() * t1);
          }
          if (b != a)
            eMat(b, a) = eMat(a, b);
        }
      }
      printIt(&eMat, "circle_fit - E:");

      Eigen::Matrix<double, 3, 6> j2Mat;  // Jacobian of min_eigen() (numerically computed)
      for (uint a = 0; a < 6; ++a) {
        const uint i = nu[a][0], j = nu[a][1];
        Matrix3d delta = Matrix3d::Zero();
        delta(i, j) = delta(j, i) = abs(aMat(i, j) * epsilon);
        j2Mat.col(a) = min_eigen3D_fast(acc, aMat + delta);
        const int sign = (j2Mat.col(a)(2) > 0) ? 1 : -1;
        j2Mat.col(a) = (j2Mat.col(a) * sign - vVec) / delta(i, j);
      }
      printIt(&j2Mat, "circle_fit - J2:");

      Matrix4d cvcMat;  // joint cov matrix of (v0,v1,v2,c)
      {
        Matrix3d t0 = j2Mat * eMat * j2Mat.transpose();
        Vector3d t1 = -t0 * r0;
        cvcMat.block(0, 0, 3, 3) = t0;
        cvcMat.block(0, 3, 3, 1) = t1;
        cvcMat.block(3, 0, 1, 3) = t1.transpose();
        Eigen::Matrix<double, 1, 1> cm1;
        Eigen::Matrix<double, 1, 1> cm3;
        cm1 = (vVec.transpose() * c0Mat * vVec);
        //      cm2 = (c0Mat.cwiseProduct(t0)).sum();
        cm3 = (r0.transpose() * t0 * r0);
        // Workaround to get things working in GPU
        const double tempC = cm1(0, 0) + (c0Mat.cwiseProduct(t0)).sum() + cm3(0, 0);
        cvcMat(3, 3) = tempC;
        // (v.transpose() * c0Mat * v) + (c0Mat.cwiseProduct(t0)).sum() + (r0.transpose() * t0 * r0);
      }
      printIt(&cvcMat, "circle_fit - Cvc:");

      Eigen::Matrix<double, 3, 4> j3Mat;  // Jacobian (v0,v1,v2,c)->(X0,Y0,R)
      {
        const double t = 1. / tempH;
        j3Mat << -v2x2_inv, 0, vVec(0) * sqr(v2x2_inv) * 2., 0, 0, -v2x2_inv, vVec(1) * sqr(v2x2_inv) * 2., 0,
            vVec(0) * v2x2_inv * t, vVec(1) * v2x2_inv * t,
            -tempH * sqr(v2x2_inv) * 2. - (2. * tempC + vVec(2)) * v2x2_inv * t, -t;
      }
      printIt(&j3Mat, "circle_fit - J3:");

      const RowVector2Nd<N> Jq = mc.transpose() * tempS * 1. / n;  // var(q)
      printIt(&Jq, "circle_fit - Jq:");

      Matrix3d cov_uvr = j3Mat * cvcMat * j3Mat.transpose() * sqr(s_inv)  // cov(X0,Y0,R)
                         + (par_uvr * par_uvr.transpose()) * (Jq * vMat * Jq.transpose());

      circle.cov = cov_uvr;
    }

    printIt(&circle.cov, "Circle cov:");
#ifdef RFIT_DEBUG
    printf("circle_fit - exit\n");
#endif
    return circle;
  }

computeRadLenUniformMaterial()

template<typename TAcc , typename VNd1 , typename VNd2 >

ALPAKA_FN_ACC ALPAKA_FN_INLINE void ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::computeRadLenUniformMaterial	(	const TAcc &	acc,
		const VNd1 &	length_values,
		VNd2 &	rad_lengths
	)

Compute the Radiation length in the uniform hypothesis

The Pixel detector, barrel and forward, is considered as an homogeneous cylinder of material, whose radiation lengths has been derived from the TDR plot that shows that 16cm correspond to 0.06 radiation lengths. Therefore one radiation length corresponds to 16cm/0.06 =~ 267 cm. All radiation lengths are computed using this unique number, in both regions, barrel and endcap.

NB: no angle corrections nor projections are computed inside this routine. It is therefore the responsibility of the caller to supply the proper lengths in input. These lengths are the path traveled by the particle along its trajectory, namely the so called S of the helix in 3D space.

Parameters

length_values

vector of incremental distances that will be translated into radiation length equivalent. Each radiation length i is computed incrementally with respect to the previous length i-1. The first length has no reference point (i.e. it has the dca).

Returns

incremental radiation lengths that correspond to each segment.

Definition at line 38 of file RiemannFit.h.

References funct::abs(), ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), dqmiolumiharvest::j, dqmiodumpmetadata::n, and parallelization::uint.

Referenced by scatter_cov_rad(), and scatterCovLine().

                                                                                      {
    // Radiation length of the pixel detector in the uniform assumption, with
    // 0.06 rad_len at 16 cm
    constexpr double xx_0_inv = 0.06 / 16.;
    uint n = length_values.rows();
    rad_lengths(0) = length_values(0) * xx_0_inv;
    for (uint j = 1; j < n; ++j) {
      rad_lengths(j) = alpaka::math::abs(acc, length_values(j) - length_values(j - 1)) * xx_0_inv;
    }
  }

cov_carttorad()

template<typename TAcc , typename M2xN , int N>

ALPAKA_FN_ACC ALPAKA_FN_INLINE VectorNd<N> ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::cov_carttorad	(	const TAcc &	acc,
		const M2xN &	p2D,
		const Matrix2Nd< N > &	cov_cart,
		const VectorNd< N > &	rad
	)

Transform covariance matrix from Cartesian coordinates (only transverse plane component) to radial coordinates (both radial and tangential component but only diagonal terms, correlation between different point are not managed).

Parameters

p2D	2D points in transverse plane.
cov_cart	covariance matrix in Cartesian coordinates.

Returns

cov_rad covariance matrix in raidal coordinate.

Warning

correlation between different point are not computed.

< in case you have (0,0) to avoid dividing by 0 radius

Definition at line 218 of file RiemannFit.h.

References ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), MillePedeFileConverter_cfg::e, mps_fire::i, N, dqmiodumpmetadata::n, sqr(), and parallelization::uint.

                                                                                   {
    constexpr uint n = N;
    VectorNd<N> cov_rad;
    const VectorNd<N> rad_inv2 = rad.cwiseInverse().array().square();
    for (uint i = 0; i < n; ++i) {
      if (rad(i) < 1.e-4)
        cov_rad(i) = cov_cart(i, i);
      else {
        cov_rad(i) = rad_inv2(i) * (cov_cart(i, i) * sqr(p2D(1, i)) + cov_cart(i + n, i + n) * sqr(p2D(0, i)) -
                                    2. * cov_cart(i, i + n) * p2D(0, i) * p2D(1, i));
      }
    }
    return cov_rad;
  }

cov_carttorad_prefit()

template<typename TAcc , typename M2xN , typename V4 , int N>

ALPAKA_FN_ACC ALPAKA_FN_INLINE VectorNd<N> ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::cov_carttorad_prefit	(	const TAcc &	acc,
		const M2xN &	p2D,
		const Matrix2Nd< N > &	cov_cart,
		V4 &	fast_fit,
		const VectorNd< N > &	rad
	)

Transform covariance matrix from Cartesian coordinates (only transverse plane component) to coordinates system orthogonal to the pre-fitted circle in each point. Further information in attached documentation.

Parameters

p2D	2D points in transverse plane.
cov_cart	covariance matrix in Cartesian coordinates.
fast_fit	fast_fit Vector4d result of the previous pre-fit structured in this form:(X0, Y0, R, tan(theta))).

Returns

cov_rad covariance matrix in the pre-fitted circle's orthogonal system.

< in case you have (0,0) to avoid dividing by 0 radius

Definition at line 250 of file RiemannFit.h.

References a, b, ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), cross2D(), MillePedeFileConverter_cfg::e, ALPAKA_ACCELERATOR_NAMESPACE::brokenline::fast_fit, mps_fire::i, N, dqmiodumpmetadata::n, sqr(), parallelization::uint, and testProducerWithPsetDescEmpty_cfi::y2.

Referenced by circleFit().

                                                                                                            {
    constexpr uint n = N;
    VectorNd<N> cov_rad;
    for (uint i = 0; i < n; ++i) {
      if (rad(i) < 1.e-4)
        cov_rad(i) = cov_cart(i, i);  // TO FIX
      else {
        Vector2d a = p2D.col(i);
        Vector2d b = p2D.col(i) - fast_fit.head(2);
        const double x2 = a.dot(b);
        const double y2 = cross2D(acc, a, b);
        const double tan_c = -y2 / x2;
        const double tan_c2 = sqr(tan_c);
        cov_rad(i) =
            1. / (1. + tan_c2) * (cov_cart(i, i) + cov_cart(i + n, i + n) * tan_c2 + 2 * cov_cart(i, i + n) * tan_c);
      }
    }
    return cov_rad;
  }

cov_radtocart()

template<typename TAcc , typename M2xN , int N>

ALPAKA_FN_ACC ALPAKA_FN_INLINE Matrix2Nd<N> ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::cov_radtocart	(	const TAcc &	acc,
		const M2xN &	p2D,
		const MatrixNd< N > &	cov_rad,
		const VectorNd< N > &	rad
	)

Transform covariance matrix from radial (only tangential component) to Cartesian coordinates (only transverse plane component).

Parameters

p2D	2D points in the transverse plane.
cov_rad	covariance matrix in radial coordinate.

Returns

cov_cart covariance matrix in Cartesian coordinates.

Definition at line 180 of file RiemannFit.h.

References ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), mps_fire::i, dqmiolumiharvest::j, N, dqmiodumpmetadata::n, printIt(), and parallelization::uint.

Referenced by circleFit().

                                                                                    {
#ifdef RFIT_DEBUG
    printf("Address of p2D: %p\n", &p2D);
#endif
    printIt(&p2D, "cov_radtocart - p2D:");
    constexpr uint n = N;
    Matrix2Nd<N> cov_cart = Matrix2Nd<N>::Zero();
    VectorNd<N> rad_inv = rad.cwiseInverse();
    printIt(&rad_inv, "cov_radtocart - rad_inv:");
    for (uint i = 0; i < n; ++i) {
      for (uint j = i; j < n; ++j) {
        cov_cart(i, j) = cov_rad(i, j) * p2D(1, i) * rad_inv(i) * p2D(1, j) * rad_inv(j);
        cov_cart(i + n, j + n) = cov_rad(i, j) * p2D(0, i) * rad_inv(i) * p2D(0, j) * rad_inv(j);
        cov_cart(i, j + n) = -cov_rad(i, j) * p2D(1, i) * rad_inv(i) * p2D(0, j) * rad_inv(j);
        cov_cart(i + n, j) = -cov_rad(i, j) * p2D(0, i) * rad_inv(i) * p2D(1, j) * rad_inv(j);
        cov_cart(j, i) = cov_cart(i, j);
        cov_cart(j + n, i + n) = cov_cart(i + n, j + n);
        cov_cart(j + n, i) = cov_cart(i, j + n);
        cov_cart(j, i + n) = cov_cart(i + n, j);
      }
    }
    return cov_cart;
  }

cross2D()

template<typename TAcc >

ALPAKA_FN_ACC ALPAKA_FN_INLINE double ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::cross2D	(	const TAcc &	acc,
		const Vector2d &	a,
		const Vector2d &	b
	)

Compute cross product of two 2D vector (assuming z component 0), returning z component of the result.

Parameters

a	first 2D vector in the product.
b	second 2D vector in the product.

Returns

z component of the cross product.

Definition at line 117 of file FitUtils.h.

References a, and b.

Referenced by cov_carttorad_prefit(), ALPAKA_ACCELERATOR_NAMESPACE::brokenline::fastFit(), fastFit(), ALPAKA_ACCELERATOR_NAMESPACE::brokenline::for(), lineFit(), and scatter_cov_rad().

                                                                                                         {
      return a.x() * b.y() - a.y() * b.x();
    }

fastFit()

template<typename TAcc , typename M3xN , typename V4 >

ALPAKA_FN_ACC ALPAKA_FN_INLINE void ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::fastFit	(	const TAcc &	acc,
		const M3xN &	hits,
		V4 &	result
	)

A very fast helix fit: it fits a circle by three points (first, middle and last point) and a line by two points (first and last).

Parameters

hits	points to be fitted

Returns

result in this form: (X0,Y0,R,tan(theta)).

Warning

points must be passed ordered (from internal layer to external) in order to maximize accuracy and do not mistake tan(theta) sign.

This fast fit is used as pre-fit which is needed for:

• weights estimation and chi2 computation in line fit (fundamental);
• weights estimation and chi2 computation in circle fit (useful);
• computation of error due to multiple scattering.

< in case b.x is 0 (2 hits with same x)

Definition at line 385 of file RiemannFit.h.

References funct::abs(), b2, nano_mu_digi_cff::bx, ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), cross2D(), ALPAKA_ACCELERATOR_NAMESPACE::brokenline::dVec, PVValHelper::dz, ALPAKA_ACCELERATOR_NAMESPACE::brokenline::eVec, HLT_2024v14_cff::flip, hfClusterShapes_cfi::hits, N, dqmiodumpmetadata::n, printIt(), mps_fire::result, sqr(), and mathSSE::sqrt().

Referenced by ALPAKA_ACCELERATOR_NAMESPACE::Kernel_FastFit< N, TrackerTraits >::operator()(), and riemannFit::helixFit< N >::operator()().

                                                                                             {
    constexpr uint32_t N = M3xN::ColsAtCompileTime;
    constexpr auto n = N;  // get the number of hits
    printIt(&hits, "Fast_fit - hits: ");

    // CIRCLE FIT
    // Make segments between middle-to-first(b) and last-to-first(c) hits
    const Vector2d bVec = hits.block(0, n / 2, 2, 1) - hits.block(0, 0, 2, 1);
    const Vector2d cVec = hits.block(0, n - 1, 2, 1) - hits.block(0, 0, 2, 1);
    printIt(&bVec, "Fast_fit - b: ");
    printIt(&cVec, "Fast_fit - c: ");
    // Compute their lengths
    auto b2 = bVec.squaredNorm();
    auto c2 = cVec.squaredNorm();
    // The algebra has been verified (MR). The usual approach has been followed:
    // * use an orthogonal reference frame passing from the first point.
    // * build the segments (chords)
    // * build orthogonal lines through mid points
    // * make a system and solve for X0 and Y0.
    // * add the initial point
    bool flip = abs(bVec.x()) < abs(bVec.y());
    auto bx = flip ? bVec.y() : bVec.x();
    auto by = flip ? bVec.x() : bVec.y();
    auto cx = flip ? cVec.y() : cVec.x();
    auto cy = flip ? cVec.x() : cVec.y();
    auto div = 2. * (cx * by - bx * cy);
    // if aligned TO FIX
    auto y0 = (cx * b2 - bx * c2) / div;
    auto x0 = (0.5 * b2 - y0 * by) / bx;
    result(0) = hits(0, 0) + (flip ? y0 : x0);
    result(1) = hits(1, 0) + (flip ? x0 : y0);
    result(2) = sqrt(sqr(x0) + sqr(y0));
    printIt(&result, "Fast_fit - result: ");

    // LINE FIT
    const Vector2d dVec = hits.block(0, 0, 2, 1) - result.head(2);
    const Vector2d eVec = hits.block(0, n - 1, 2, 1) - result.head(2);
    printIt(&eVec, "Fast_fit - e: ");
    printIt(&dVec, "Fast_fit - d: ");
    // Compute the arc-length between first and last point: L = R * theta = R * atan (tan (Theta) )
    auto dr = result(2) * atan2(cross2D(acc, dVec, eVec), dVec.dot(eVec));
    // Simple difference in Z between last and first hit
    auto dz = hits(2, n - 1) - hits(2, 0);

    result(3) = (dr / dz);

#ifdef RFIT_DEBUG
    printf("Fast_fit: [%f, %f, %f, %f]\n", result(0), result(1), result(2), result(3));
#endif
  }

fromCircleToPerigee()

template<typename TAcc >

ALPAKA_FN_ACC ALPAKA_FN_INLINE void ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::fromCircleToPerigee	(	const TAcc &	acc,
		CircleFit &	circle
	)

Transform circle parameter from (X0,Y0,R) to (phi,Tip,q/R) and consequently covariance matrix.

Parameters

circle_uvr

parameter (X0,Y0,R), covariance matrix to be transformed and particle charge.

Definition at line 239 of file FitUtils.h.

References sqr(), and mathSSE::sqrt().

Referenced by ALPAKA_ACCELERATOR_NAMESPACE::Kernel_LineFit< N, TrackerTraits >::operator()().

                                                                                                {
      Vector3d par_pak;
      const double temp0 = circle.par.head(2).squaredNorm();
      const double temp1 = alpaka::math::sqrt(acc, temp0);
      par_pak << alpaka::math::atan2(acc, circle.qCharge * circle.par(0), -circle.qCharge * circle.par(1)),
          circle.qCharge * (temp1 - circle.par(2)), circle.qCharge / circle.par(2);

      const double temp2 = sqr(circle.par(0)) * 1. / temp0;
      const double temp3 = 1. / temp1 * circle.qCharge;
      Matrix3d j4Mat;
      j4Mat << -circle.par(1) * temp2 * 1. / sqr(circle.par(0)), temp2 * 1. / circle.par(0), 0., circle.par(0) * temp3,
          circle.par(1) * temp3, -circle.qCharge, 0., 0., -circle.qCharge / (circle.par(2) * circle.par(2));
      circle.cov = j4Mat * circle.cov * j4Mat.transpose();

      circle.par = par_pak;
    }

lineFit()

template<typename TAcc , typename M3xN , typename M6xN , typename V4 >

ALPAKA_FN_ACC ALPAKA_FN_INLINE LineFit ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::lineFit	(	const TAcc &	acc,
		const M3xN &	hits,
		const M6xN &	hits_ge,
		const CircleFit &	circle,
		const V4 &	fast_fit,
		const double	bField,
		const bool	error
	)

Perform an ordinary least square fit in the s-z plane to compute the parameters cotTheta and Zip.

The fit is performed in the rotated S3D-Z' plane, following the formalism of Frodesen, Chapter 10, p. 259.

The system has been rotated to both try to use the combined errors in s-z along Z', as errors in the Y direction and to avoid the patological case of degenerate lines with angular coefficient m = +/- inf.

The rotation is using the information on the theta angle computed in the fast fit. The rotation is such that the S3D axis will be the X-direction, while the rotated Z-axis will be the Y-direction. This pretty much follows what is done in the same fit in the Broken Line approach.

Definition at line 796 of file RiemannFit.h.

References Calorimetry_cff::bField, nano_mu_local_reco_cff::chi2, ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), funct::cos(), cross(), cross2D(), alignCSCRings::d_x, alignCSCRings::d_y, dot(), relativeConstraints::error, ALPAKA_ACCELERATOR_NAMESPACE::brokenline::fast_fit, hfClusterShapes_cfi::hits, mps_fire::i, l1tstage2_dqm_sourceclient-live_cfg::invert, mps_splice::line, M_PI, N, dqmiodumpmetadata::n, printIt(), makeMuonMisalignmentScenario::rot, scatterCovLine(), funct::sin(), mkfit::Const::sol, sqr(), theta(), createJobs::tmp, and parallelization::uint.

Referenced by ALPAKA_ACCELERATOR_NAMESPACE::Kernel_LineFit< N, TrackerTraits >::operator()(), and riemannFit::helixFit< N >::operator()().

                                                                   {
    constexpr uint32_t N = M3xN::ColsAtCompileTime;
    constexpr auto n = N;
    double theta = -circle.qCharge * atan(fast_fit(3));
    theta = theta < 0. ? theta + M_PI : theta;

    // Prepare the Rotation Matrix to rotate the points
    Eigen::Matrix<double, 2, 2> rot;
    rot << sin(theta), cos(theta), -cos(theta), sin(theta);

    // PROJECTION ON THE CILINDER
    //
    // p2D will be:
    // [s1, s2, s3, ..., sn]
    // [z1, z2, z3, ..., zn]
    // s values will be ordinary x-values
    // z values will be ordinary y-values

    Matrix2xNd<N> p2D = Matrix2xNd<N>::Zero();
    Eigen::Matrix<double, 2, 6> jxMat;

#ifdef RFIT_DEBUG
    printf("Line_fit - B: %g\n", bField);
    printIt(&hits, "Line_fit points: ");
    printIt(&hits_ge, "Line_fit covs: ");
    printIt(&rot, "Line_fit rot: ");
#endif
    // x & associated Jacobian
    // cfr https://indico.cern.ch/event/663159/contributions/2707659/attachments/1517175/2368189/Riemann_fit.pdf
    // Slide 11
    // a ==> -o i.e. the origin of the circle in XY plane, negative
    // b ==> p i.e. distances of the points wrt the origin of the circle.
    const Vector2d oVec(circle.par(0), circle.par(1));

    // associated Jacobian, used in weights and errors computation
    Matrix6d covMat = Matrix6d::Zero();
    Matrix2d cov_sz[N];
    for (uint i = 0; i < n; ++i) {
      Vector2d pVec = hits.block(0, i, 2, 1) - oVec;
      const double cross = cross2D(acc, -oVec, pVec);
      const double dot = (-oVec).dot(pVec);
      // atan2(cross, dot) give back the angle in the transverse plane so tha the
      // final equation reads: x_i = -q*R*theta (theta = angle returned by atan2)
      const double tempQAtan2 = -circle.qCharge * atan2(cross, dot);
      //    p2D.coeffRef(1, i) = atan2_ * circle.par(2);
      p2D(0, i) = tempQAtan2 * circle.par(2);

      // associated Jacobian, used in weights and errors- computation
      const double temp0 = -circle.qCharge * circle.par(2) * 1. / (sqr(dot) + sqr(cross));
      double d_X0 = 0., d_Y0 = 0., d_R = 0.;  // good approximation for big pt and eta
      if (error) {
        d_X0 = -temp0 * ((pVec(1) + oVec(1)) * dot - (pVec(0) - oVec(0)) * cross);
        d_Y0 = temp0 * ((pVec(0) + oVec(0)) * dot - (oVec(1) - pVec(1)) * cross);
        d_R = tempQAtan2;
      }
      const double d_x = temp0 * (oVec(1) * dot + oVec(0) * cross);
      const double d_y = temp0 * (-oVec(0) * dot + oVec(1) * cross);
      jxMat << d_X0, d_Y0, d_R, d_x, d_y, 0., 0., 0., 0., 0., 0., 1.;

      covMat.block(0, 0, 3, 3) = circle.cov;
      covMat(3, 3) = hits_ge.col(i)[0];                 // x errors
      covMat(4, 4) = hits_ge.col(i)[2];                 // y errors
      covMat(5, 5) = hits_ge.col(i)[5];                 // z errors
      covMat(3, 4) = covMat(4, 3) = hits_ge.col(i)[1];  // cov_xy
      covMat(3, 5) = covMat(5, 3) = hits_ge.col(i)[3];  // cov_xz
      covMat(4, 5) = covMat(5, 4) = hits_ge.col(i)[4];  // cov_yz
      Matrix2d tmp = jxMat * covMat * jxMat.transpose();
      cov_sz[i].noalias() = rot * tmp * rot.transpose();
    }
    // Math of d_{X0,Y0,R,x,y} all verified by hand
    p2D.row(1) = hits.row(2);

    // The following matrix will contain errors orthogonal to the rotated S
    // component only, with the Multiple Scattering properly treated!!
    MatrixNd<N> cov_with_ms;
    scatterCovLine(acc, cov_sz, fast_fit, p2D.row(0), p2D.row(1), theta, bField, cov_with_ms);
#ifdef RFIT_DEBUG
    printIt(cov_sz, "line_fit - cov_sz:");
    printIt(&cov_with_ms, "line_fit - cov_with_ms: ");
#endif

    // Rotate Points with the shape [2, n]
    Matrix2xNd<N> p2D_rot = rot * p2D;

#ifdef RFIT_DEBUG
    printf("Fast fit Tan(theta): %g\n", fast_fit(3));
    printf("Rotation angle: %g\n", theta);
    printIt(&rot, "Rotation Matrix:");
    printIt(&p2D, "Original Hits(s,z):");
    printIt(&p2D_rot, "Rotated hits(S3D, Z'):");
    printIt(&rot, "Rotation Matrix:");
#endif

    // Build the A Matrix
    Matrix2xNd<N> aMat;
    aMat << MatrixXd::Ones(1, n), p2D_rot.row(0);  // rotated s values

#ifdef RFIT_DEBUG
    printIt(&aMat, "A Matrix:");
#endif

    // Build A^T V-1 A, where V-1 is the covariance of only the Y components.
    MatrixNd<N> vyInvMat;
    math::cholesky::invert(cov_with_ms, vyInvMat);
    // MatrixNd<N> vyInvMat = cov_with_ms.inverse();
    Eigen::Matrix<double, 2, 2> covParamsMat = aMat * vyInvMat * aMat.transpose();
    // Compute the Covariance Matrix of the fit parameters
    math::cholesky::invert(covParamsMat, covParamsMat);

    // Now Compute the Parameters in the form [2,1]
    // The first component is q.
    // The second component is m.
    Eigen::Matrix<double, 2, 1> sol = covParamsMat * aMat * vyInvMat * p2D_rot.row(1).transpose();

#ifdef RFIT_DEBUG
    printIt(&sol, "Rotated solutions:");
#endif

    // We need now to transfer back the results in the original s-z plane
    const auto sinTheta = sin(theta);
    const auto cosTheta = cos(theta);
    auto common_factor = 1. / (sinTheta - sol(1, 0) * cosTheta);
    Eigen::Matrix<double, 2, 2> jMat;
    jMat << 0., common_factor * common_factor, common_factor, sol(0, 0) * cosTheta * common_factor * common_factor;

    double tempM = common_factor * (sol(1, 0) * sinTheta + cosTheta);
    double tempQ = common_factor * sol(0, 0);
    auto cov_mq = jMat * covParamsMat * jMat.transpose();

    VectorNd<N> res = p2D_rot.row(1).transpose() - aMat.transpose() * sol;
    double chi2 = res.transpose() * vyInvMat * res;

    LineFit line;
    line.par << tempM, tempQ;
    line.cov << cov_mq;
    line.chi2 = chi2;

#ifdef RFIT_DEBUG
    printf("Common_factor: %g\n", common_factor);
    printIt(&jMat, "Jacobian:");
    printIt(&sol, "Rotated solutions:");
    printIt(&covParamsMat, "Cov_params:");
    printIt(&cov_mq, "Rotated Covariance Matrix:");
    printIt(&(line.par), "Real Parameters:");
    printIt(&(line.cov), "Real Covariance Matrix:");
    printf("Chi2: %g\n", chi2);
#endif

    return line;
  }

loadCovariance()

template<typename TAcc , typename M6xNf , typename M3xNd >

ALPAKA_FN_ACC void ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::loadCovariance	(	const TAcc &	acc,
		M6xNf const &	ge,
		M3xNd &	hits_cov
	)

Definition at line 158 of file FitUtils.h.

References ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), mps_fire::i, dqmiolumiharvest::j, and MainPageGenerator::l.

                                                                                         {
      // Index numerology:
      // i: index of the hits/point (0,..,3)
      // j: index of space component (x,y,z)
      // l: index of space components (x,y,z)
      // ge is always in sync with the index i and is formatted as:
      // ge[] ==> [xx, xy, yy, xz, yz, zz]
      // in (j,l) notation, we have:
      // ge[] ==> [(0,0), (0,1), (1,1), (0,2), (1,2), (2,2)]
      // so the index ge_idx corresponds to the matrix elements:
      // | 0  1  3 |
      // | 1  2  4 |
      // | 3  4  5 |
      constexpr uint32_t hits_in_fit = M6xNf::ColsAtCompileTime;
      for (uint32_t i = 0; i < hits_in_fit; ++i) {
        {
          constexpr uint32_t ge_idx = 0, j = 0, l = 0;
          hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) = ge.col(i)[ge_idx];
        }
        {
          constexpr uint32_t ge_idx = 2, j = 1, l = 1;
          hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) = ge.col(i)[ge_idx];
        }
        {
          constexpr uint32_t ge_idx = 5, j = 2, l = 2;
          hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) = ge.col(i)[ge_idx];
        }
        {
          constexpr uint32_t ge_idx = 1, j = 1, l = 0;
          hits_cov(i + l * hits_in_fit, i + j * hits_in_fit) = hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) =
              ge.col(i)[ge_idx];
        }
        {
          constexpr uint32_t ge_idx = 3, j = 2, l = 0;
          hits_cov(i + l * hits_in_fit, i + j * hits_in_fit) = hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) =
              ge.col(i)[ge_idx];
        }
        {
          constexpr uint32_t ge_idx = 4, j = 2, l = 1;
          hits_cov(i + l * hits_in_fit, i + j * hits_in_fit) = hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) =
              ge.col(i)[ge_idx];
        }
      }
    }

loadCovariance2D()

template<typename TAcc , typename M6xNf , typename M2Nd >

ALPAKA_FN_ACC void ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::loadCovariance2D	(	const TAcc &	acc,
		M6xNf const &	ge,
		M2Nd &	hits_cov
	)

load error in CMSSW format to our formalism

Definition at line 126 of file FitUtils.h.

References ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), mps_fire::i, dqmiolumiharvest::j, and MainPageGenerator::l.

Referenced by ALPAKA_ACCELERATOR_NAMESPACE::Kernel_CircleFit< N, TrackerTraits >::operator()(), and riemannFit::helixFit< N >::operator()().

                                                                                          {
      // Index numerology:
      // i: index of the hits/point (0,..,3)
      // j: index of space component (x,y,z)
      // l: index of space components (x,y,z)
      // ge is always in sync with the index i and is formatted as:
      // ge[] ==> [xx, xy, yy, xz, yz, zz]
      // in (j,l) notation, we have:
      // ge[] ==> [(0,0), (0,1), (1,1), (0,2), (1,2), (2,2)]
      // so the index ge_idx corresponds to the matrix elements:
      // | 0  1  3 |
      // | 1  2  4 |
      // | 3  4  5 |
      constexpr uint32_t hits_in_fit = M6xNf::ColsAtCompileTime;
      for (uint32_t i = 0; i < hits_in_fit; ++i) {
        {
          constexpr uint32_t ge_idx = 0, j = 0, l = 0;
          hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) = ge.col(i)[ge_idx];
        }
        {
          constexpr uint32_t ge_idx = 2, j = 1, l = 1;
          hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) = ge.col(i)[ge_idx];
        }
        {
          constexpr uint32_t ge_idx = 1, j = 1, l = 0;
          hits_cov(i + l * hits_in_fit, i + j * hits_in_fit) = hits_cov(i + j * hits_in_fit, i + l * hits_in_fit) =
              ge.col(i)[ge_idx];
        }
      }
    }

min_eigen2D()

template<typename TAcc >

ALPAKA_FN_ACC ALPAKA_FN_INLINE Vector2d ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::min_eigen2D	(	const TAcc &	acc,
		const Matrix2d &	aMat,
		double &	chi2
	)

2D version of min_eigen3D().

Parameters

aMat	the Matrix you want to know eigenvector and eigenvalue.
chi2	the double were the chi2-related quantity will be stored

Returns

the eigenvector associated to the minimum eigenvalue. The computedDirect() method of SelfAdjointEigenSolver for 2x2 Matrix do not use special math function (just sqrt) therefore it doesn't speed up significantly in single precision.

Definition at line 363 of file RiemannFit.h.

References nano_mu_local_reco_cff::chi2.

                                                                                                           {
    Eigen::SelfAdjointEigenSolver<Matrix2d> solver(2);
    solver.computeDirect(aMat);
    int min_index;
    chi2 = solver.eigenvalues().minCoeff(&min_index);
    return solver.eigenvectors().col(min_index);
  }

min_eigen3D()

template<typename TAcc >

ALPAKA_FN_ACC ALPAKA_FN_INLINE Vector3d ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::min_eigen3D	(	const TAcc &	acc,
		const Matrix3d &	A,
		double &	chi2
	)

Compute the eigenvector associated to the minimum eigenvalue.

Parameters

A	the Matrix you want to know eigenvector and eigenvalue.
chi2	the double were the chi2-related quantity will be stored.

Returns

the eigenvector associated to the minimum eigenvalue.

Warning

double precision is needed for a correct assessment of chi2.

The minimus eigenvalue is related to chi2. We exploit the fact that the matrix is symmetrical and small (2x2 for line fit and 3x3 for circle fit), so the SelfAdjointEigenSolver from Eigen library is used, with the computedDirect method (available only for 2x2 and 3x3 Matrix) wich computes eigendecomposition of given matrix using a fast closed-form algorithm. For this optimization the matrix type must be known at compiling time.

Definition at line 319 of file RiemannFit.h.

References A, and nano_mu_local_reco_cff::chi2.

Referenced by circleFit().

                                                                                                        {
#ifdef RFIT_DEBUG
    printf("min_eigen3D - enter\n");
#endif
    Eigen::SelfAdjointEigenSolver<Matrix3d> solver(3);
    solver.computeDirect(A);
    int min_index;
    chi2 = solver.eigenvalues().minCoeff(&min_index);
#ifdef RFIT_DEBUG
    printf("min_eigen3D - exit\n");
#endif
    return solver.eigenvectors().col(min_index);
  }

min_eigen3D_fast()

template<typename TAcc >

ALPAKA_FN_ACC ALPAKA_FN_INLINE Vector3d ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::min_eigen3D_fast	(	const TAcc &	acc,
		const Matrix3d &	A
	)

A faster version of min_eigen3D() where double precision is not needed.

Parameters

A	the Matrix you want to know eigenvector and eigenvalue.
chi2	the double were the chi2-related quantity will be stored

Returns

the eigenvector associated to the minimum eigenvalue. The computedDirect() method of SelfAdjointEigenSolver for 3x3 Matrix indeed, use trigonometry function (it solves a third degree equation) which speed up in single precision.

Definition at line 345 of file RiemannFit.h.

References A.

Referenced by circleFit().

                                                                                               {
    Eigen::SelfAdjointEigenSolver<Matrix3f> solver(3);
    solver.computeDirect(A.cast<float>());
    int min_index;
    solver.eigenvalues().minCoeff(&min_index);
    return solver.eigenvectors().col(min_index).cast<double>();
  }

par_uvrtopak()

template<typename TAcc >

ALPAKA_FN_ACC ALPAKA_FN_INLINE void ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::par_uvrtopak	(	const TAcc &	acc,
		CircleFit &	circle,
		const double	B,
		const bool	error
	)

Transform circle parameter from (X0,Y0,R) to (phi,Tip,p_t) and consequently covariance matrix.

Parameters

circle_uvr	parameter (X0,Y0,R), covariance matrix to be transformed and particle charge.
B	magnetic field in Gev/cm/c unit.
error	flag for errors computation.

Definition at line 212 of file FitUtils.h.

References B, relativeConstraints::error, sqr(), and mathSSE::sqrt().

Referenced by riemannFit::helixFit< N >::operator()().

                                                                       {
      Vector3d par_pak;
      const double temp0 = circle.par.head(2).squaredNorm();
      const double temp1 = alpaka::math::sqrt(acc, temp0);
      par_pak << alpaka::math::atan2(acc, circle.qCharge * circle.par(0), -circle.qCharge * circle.par(1)),
          circle.qCharge * (temp1 - circle.par(2)), circle.par(2) * B;
      if (error) {
        const double temp2 = sqr(circle.par(0)) * 1. / temp0;
        const double temp3 = 1. / temp1 * circle.qCharge;
        Matrix3d j4Mat;
        j4Mat << -circle.par(1) * temp2 * 1. / sqr(circle.par(0)), temp2 * 1. / circle.par(0), 0.,
            circle.par(0) * temp3, circle.par(1) * temp3, -circle.qCharge, 0., 0., B;
        circle.cov = j4Mat * circle.cov * j4Mat.transpose();
      }
      circle.par = par_pak;
    }

printIt()

template<typename TAcc , class C >

ALPAKA_FN_ACC void ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::printIt	(	const TAcc &	acc,
		C *	m,
		const char *	prefix = ""
	)

Definition at line 90 of file FitUtils.h.

References DummyCfis::c, visualization-live-secondInstance_cfg::m, hcallasereventfilter2012_cfi::prefix, and parallelization::uint.

Referenced by circleFit(), cov_radtocart(), PrintSensitive::dumpTouch(), PrintMTDSens::dumpTouch(), fastFit(), lineFit(), scatter_cov_rad(), scatterCovLine(), and fwlite::Scanner< Collection >::setPrintFullEventId().

                                                                               {
#ifdef RFIT_DEBUG
      for (uint r = 0; r < m->rows(); ++r) {
        for (uint c = 0; c < m->cols(); ++c) {
          printf("%s Matrix(%d,%d) = %g\n", prefix, r, c, (*m)(r, c));
        }
      }
#endif
    }

scatter_cov_rad()

template<typename TAcc , typename M2xN , typename V4 , int N>

ALPAKA_FN_ACC ALPAKA_FN_INLINE MatrixNd<N> ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::scatter_cov_rad	(	const TAcc &	acc,
		const M2xN &	p2D,
		const V4 &	fast_fit,
		VectorNd< N > const &	rad,
		double	B
	)

Compute the covariance matrix (in radial coordinates) of points in the transverse plane due to multiple Coulomb scattering.

Parameters

p2D	2D points in the transverse plane.
fast_fit	fast_fit Vector4d result of the previous pre-fit structured in this form:(X0, Y0, R, Tan(Theta))).
B	magnetic field use to compute p

Returns

scatter_cov_rad errors due to multiple scattering.

Warning

input points must be ordered radially from the detector center (from inner layer to outer ones; points on the same layer must ordered too).

Only the tangential component is computed (the radial one is negligible).

Definition at line 134 of file RiemannFit.h.

References funct::abs(), B, computeRadLenUniformMaterial(), ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), cross(), cross2D(), dot(), ALPAKA_ACCELERATOR_NAMESPACE::brokenline::fast_fit, mps_fire::i, dqmdumpme::k, MainPageGenerator::l, M_PI, SiStripPI::min, N, dqmiodumpmetadata::n, printIt(), funct::sin(), sqr(), mathSSE::sqrt(), theta(), and parallelization::uint.

Referenced by circleFit().

                                                                                              {
    constexpr uint n = N;
    double p_t = alpaka::math::min(acc, 20., fast_fit(2) * B);  // limit pt to avoid too small error!!!
    double p_2 = p_t * p_t * (1. + 1. / sqr(fast_fit(3)));
    double theta = atan(fast_fit(3));
    theta = theta < 0. ? theta + M_PI : theta;
    VectorNd<N> s_values;
    VectorNd<N> rad_lengths;
    const Vector2d oVec(fast_fit(0), fast_fit(1));

    // associated Jacobian, used in weights and errors computation
    for (uint i = 0; i < n; ++i) {  // x
      Vector2d pVec = p2D.block(0, i, 2, 1) - oVec;
      const double cross = cross2D(acc, -oVec, pVec);
      const double dot = (-oVec).dot(pVec);
      const double tempAtan2 = atan2(cross, dot);
      s_values(i) = alpaka::math::abs(acc, tempAtan2 * fast_fit(2));
    }
    computeRadLenUniformMaterial(acc, s_values * sqrt(1. + 1. / sqr(fast_fit(3))), rad_lengths);
    MatrixNd<N> scatter_cov_rad = MatrixNd<N>::Zero();
    VectorNd<N> sig2 = (1. + 0.038 * rad_lengths.array().log()).abs2() * rad_lengths.array();
    sig2 *= 0.000225 / (p_2 * sqr(sin(theta)));
    for (uint k = 0; k < n; ++k) {
      for (uint l = k; l < n; ++l) {
        for (uint i = 0; i < uint(alpaka::math::min(acc, k, l)); ++i) {
          scatter_cov_rad(k, l) += (rad(k) - rad(i)) * (rad(l) - rad(i)) * sig2(i);
        }
        scatter_cov_rad(l, k) = scatter_cov_rad(k, l);
      }
    }
#ifdef RFIT_DEBUG
    riemannFit::printIt(&scatter_cov_rad, "Scatter_cov_rad - scatter_cov_rad: ");
#endif
    return scatter_cov_rad;
  }

scatterCovLine()

template<typename TAcc , typename V4 , typename VNd1 , typename VNd2 , int N>

ALPAKA_FN_ACC ALPAKA_FN_INLINE auto ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::scatterCovLine	(	const TAcc &	acc,
		Matrix2d const *	cov_sz,
		const V4 &	fast_fit,
		VNd1 const &	s_arcs,
		VNd2 const &	z_values,
		const double	theta,
		const double	bField,
		MatrixNd< N > &	ret
	)

Compute the covariance matrix along cartesian S-Z of points due to multiple Coulomb scattering to be used in the line_fit, for the barrel and forward cases. The input covariance matrix is in the variables s-z, original and unrotated. The multiple scattering component is computed in the usual linear approximation, using the 3D path which is computed as the squared root of the squared sum of the s and z components passed in. Internally a rotation by theta is performed and the covariance matrix returned is the one in the direction orthogonal to the rotated S3D axis, i.e. along the rotated Z axis. The choice of the rotation is not arbitrary, but derived from the fact that putting the horizontal axis along the S3D direction allows the usage of the ordinary least squared fitting techiques with the trivial parametrization y = mx + q, avoiding the patological case with m = +/- inf, that would correspond to the case at eta = 0.

Definition at line 71 of file RiemannFit.h.

References funct::abs(), Calorimetry_cff::bField, computeRadLenUniformMaterial(), ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), ALPAKA_ACCELERATOR_NAMESPACE::brokenline::fast_fit, mps_fire::i, dqmdumpme::k, MainPageGenerator::l, SiStripPI::min, N, dqmiodumpmetadata::n, printIt(), runTheMatrix::ret, sqr(), createJobs::tmp, and parallelization::uint.

Referenced by lineFit().

                                                                       {
#ifdef RFIT_DEBUG
    riemannFit::printIt(&s_arcs, "Scatter_cov_line - s_arcs: ");
#endif
    constexpr uint n = N;
    double p_t = alpaka::math::min(acc, 20., fast_fit(2) * bField);  // limit pt to avoid too small error!!!
    double p_2 = p_t * p_t * (1. + 1. / sqr(fast_fit(3)));
    VectorNd<N> rad_lengths_S;
    // See documentation at http://eigen.tuxfamily.org/dox/group__TutorialArrayClass.html
    // Basically, to perform cwise operations on Matrices and Vectors, you need
    // to transform them into Array-like objects.
    VectorNd<N> s_values = s_arcs.array() * s_arcs.array() + z_values.array() * z_values.array();
    s_values = s_values.array().sqrt();
    computeRadLenUniformMaterial(acc, s_values, rad_lengths_S);
    VectorNd<N> sig2_S;
    sig2_S = .000225 / p_2 * (1. + 0.038 * rad_lengths_S.array().log()).abs2() * rad_lengths_S.array();
#ifdef RFIT_DEBUG
    riemannFit::printIt(cov_sz, "Scatter_cov_line - cov_sz: ");
#endif
    Matrix2Nd<N> tmp = Matrix2Nd<N>::Zero();
    for (uint k = 0; k < n; ++k) {
      tmp(k, k) = cov_sz[k](0, 0);
      tmp(k + n, k + n) = cov_sz[k](1, 1);
      tmp(k, k + n) = tmp(k + n, k) = cov_sz[k](0, 1);
    }
    for (uint k = 0; k < n; ++k) {
      for (uint l = k; l < n; ++l) {
        for (uint i = 0; i < uint(alpaka::math::min(acc, k, l)); ++i) {
          tmp(k + n, l + n) += alpaka::math::abs(acc, s_values(k) - s_values(i)) *
                               alpaka::math::abs(acc, s_values(l) - s_values(i)) * sig2_S(i);
        }
        tmp(l + n, k + n) = tmp(k + n, l + n);
      }
    }
    // We are interested only in the errors orthogonal to the rotated s-axis
    // which, in our formalism, are in the lower square matrix.
#ifdef RFIT_DEBUG
    riemannFit::printIt(&tmp, "Scatter_cov_line - tmp: ");
#endif
    ret = tmp.block(n, n, n, n);
  }

sqr()

template<typename T >

constexpr T ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::sqr

(

const T

a

)

raise to square.

Definition at line 104 of file FitUtils.h.

References a.

Referenced by ALPAKA_ACCELERATOR_NAMESPACE::brokenline::circleFit(), circleFit(), cov_carttorad(), cov_carttorad_prefit(), fastFit(), fromCircleToPerigee(), ALPAKA_ACCELERATOR_NAMESPACE::brokenline::lineFit(), lineFit(), ALPAKA_ACCELERATOR_NAMESPACE::brokenline::matrixC_u(), ALPAKA_ACCELERATOR_NAMESPACE::brokenline::multScatt(), ALPAKA_ACCELERATOR_NAMESPACE::brokenline::helixFit< n >::operator()(), par_uvrtopak(), ALPAKA_ACCELERATOR_NAMESPACE::brokenline::rotationMatrix(), scatter_cov_rad(), scatterCovLine(), and ALPAKA_ACCELERATOR_NAMESPACE::brokenline::translateKarimaki().

                               {
      return a * a;
    }

weightCircle()

template<typename TAcc , int N>

ALPAKA_FN_ACC ALPAKA_FN_INLINE VectorNd<N> ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::weightCircle	(	const TAcc &	acc,
		const MatrixNd< N > &	cov_rad_inv
	)

Compute the points' weights' vector for the circle fit when multiple scattering is managed. Further information in attached documentation.

Parameters

cov_rad_inv

covariance matrix inverse in radial coordinated (or, beter, pre-fitted circle's orthogonal system).

Returns

weight VectorNd points' weights' vector.

Definition at line 284 of file RiemannFit.h.

Referenced by circleFit().

                                                                                                           {
    return cov_rad_inv.colwise().sum().transpose();
  }