riemannFit Namespace Reference

Classes
struct	CircleFit
class	helixFit
	Helix fit by three step: -fast pre-fit (see Fast_fit() for further info); -circle fit of hits projected in the transverse plane by Riemann-Chernov algorithm (see Circle_fit() for further info); -line fit of hits projected on cylinder surface by orthogonal distance regression (see Line_fit for further info). Points must be passed ordered (from inner to outer layer). More...
struct	HelixFit
struct	LineFit

Typedefs
template<int N>
using	Array2xNd = Eigen::Array< double, 2, N >
template<int N>
using	ArrayNd = Eigen::Array< double, N, N >
using	Map3x4d = Eigen::Map< Matrix3x4d, 0, Eigen::Stride< 3 *stride, stride > >
template<int N>
using	Map3xNd = Eigen::Map< Matrix3xNd< N >, 0, Eigen::Stride< 3 *stride, stride > >
using	Map4d = Eigen::Map< Vector4d, 0, Eigen::InnerStride< stride > >
using	Map6x4f = Eigen::Map< Matrix6x4f, 0, Eigen::Stride< 6 *stride, stride > >
template<int N>
using	Map6xNf = Eigen::Map< Matrix6xNf< N >, 0, Eigen::Stride< 6 *stride, stride > >
using	Matrix2d = Eigen::Matrix2d
template<int N>
using	Matrix2Nd = Eigen::Matrix< double, 2 *N, 2 *N >
using	Matrix2x3d = Eigen::Matrix< double, 2, 3 >
template<int N>
using	Matrix2xNd = Eigen::Matrix< double, 2, N >
using	Matrix3d = Eigen::Matrix3d
using	Matrix3f = Eigen::Matrix3f
template<int N>
using	Matrix3Nd = Eigen::Matrix< double, 3 *N, 3 *N >
using	Matrix3x4d = Eigen::Matrix< double, 3, 4 >
template<int N>
using	Matrix3xNd = Eigen::Matrix< double, 3, N >
using	Matrix4d = Eigen::Matrix4d
using	Matrix5d = Eigen::Matrix< double, 5, 5 >
using	Matrix6d = Eigen::Matrix< double, 6, 6 >
using	Matrix6x4f = Eigen::Matrix< float, 6, 4 >
template<int N>
using	Matrix6xNf = Eigen::Matrix< float, 6, N >
template<int N>
using	MatrixNd = Eigen::Matrix< double, N, N >
template<int N>
using	MatrixNplusONEd = Eigen::Matrix< double, N+1, N+1 >
template<int N>
using	MatrixNx3d = Eigen::Matrix< double, N, 3 >
template<int N>
using	MatrixNx5d = Eigen::Matrix< double, N, 5 >
using	MatrixXd = Eigen::MatrixXd
template<int N>
using	RowVector2Nd = Eigen::Matrix< double, 1, 2 *N >
template<int N>
using	RowVectorNd = Eigen::Matrix< double, 1, 1, N >
using	Vector2d = Eigen::Vector2d
template<int N>
using	Vector2Nd = Eigen::Matrix< double, 2 *N, 1 >
using	Vector3d = Eigen::Vector3d
using	Vector3f = Eigen::Vector3f
template<int N>
using	Vector3Nd = Eigen::Matrix< double, 3 *N, 1 >
using	Vector4d = Eigen::Vector4d
using	Vector4f = Eigen::Vector4f
using	Vector5d = Eigen::Matrix< double, 5, 1 >
using	Vector6f = Eigen::Matrix< double, 6, 1 >
template<int N>
using	VectorNd = Eigen::Matrix< double, N, 1 >
template<int N>
using	VectorNplusONEd = Eigen::Matrix< double, N+1, 1 >
using	VectorXd = Eigen::VectorXd

Functions
template<typename M2xN >
__host__ __device__ int32_t	charge (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 M2xN , typename V4 , int N>
__host__ __device__ CircleFit	circleFit (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 VNd1 , typename VNd2 >
__host__ __device__ void	computeRadLenUniformMaterial (const VNd1 &length_values, VNd2 &rad_lengths)
template<typename M2xN , int N>
__host__ __device__ VectorNd< N >	cov_carttorad (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 M2xN , typename V4 , int N>
__host__ __device__ VectorNd< N >	cov_carttorad_prefit (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 M2xN , int N>
__host__ __device__ Matrix2Nd< N >	cov_radtocart (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...
__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. More...
template<typename M3xN , typename V4 >
__host__ __device__ void	fastFit (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...
__host__ __device__ void	fromCircleToPerigee (CircleFit &circle)
	Transform circle parameter from (X0,Y0,R) to (phi,Tip,q/R) and consequently covariance matrix. More...
template<int N>
HelixFit	helixFit (const Matrix3xNd< N > &hits, const Eigen::Matrix< float, 6, N > &hits_ge, const double bField, const bool error)
	Helix fit by three step: -fast pre-fit (see Fast_fit() for further info); -circle fit of hits projected in the transverse plane by Riemann-Chernov algorithm (see Circle_fit() for further info); -line fit of hits projected on cylinder surface by orthogonal distance regression (see Line_fit for further info). Points must be passed ordered (from inner to outer layer). More...
template<typename M3xN , typename M6xN , typename V4 >
__host__ __device__ LineFit	lineFit (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 M6xNf , typename M3xNd >
__host__ __device__ void	loadCovariance (M6xNf const &ge, M3xNd &hits_cov)
template<typename M6xNf , typename M2Nd >
__host__ __device__ void	loadCovariance2D (M6xNf const &ge, M2Nd &hits_cov)
__host__ __device__ Vector2d	min_eigen2D (const Matrix2d &aMat, double &chi2)
	2D version of min_eigen3D(). More...
__host__ __device__ Vector3d	min_eigen3D (const Matrix3d &A, double &chi2)
	Compute the eigenvector associated to the minimum eigenvalue. More...
__host__ __device__ Vector3d	min_eigen3D_fast (const Matrix3d &A)
	A faster version of min_eigen3D() where double precision is not needed. More...
__host__ __device__ void	par_uvrtopak (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<class C >
__host__ __device__ void	printIt (C *m, const char *prefix="")
template<auto Start, auto End, auto Inc, class F >
constexpr void	rolling_fits (F &&f)
template<typename M2xN , typename V4 , int N>
__host__ __device__ MatrixNd< N >	scatter_cov_rad (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 V4 , typename VNd1 , typename VNd2 , int N>
__host__ __device__ auto	scatterCovLine (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 VI5 , typename MI5 , typename VO5 , typename MO5 >
void	transformToPerigeePlane (VI5 const &ip, MI5 const &icov, VO5 &op, MO5 &ocov)
template<typename VI5 , typename MI5 , typename VO5 , typename MO5 >
__host__ __device__ void	transformToPerigeePlane (VI5 const &ip, MI5 const &icov, VO5 &op, MO5 &ocov)
template<int N>
__host__ __device__ VectorNd< N >	weightCircle (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...

Variables
constexpr double	epsilon = 1.e-4
	used in numerical derivative (J2 in Circle_fit()) More...
constexpr uint32_t	maxNumberOfConcurrentFits = 32 * 1024
constexpr uint32_t	stride = maxNumberOfConcurrentFits

Typedef Documentation

Array2xNd

template<int N>

using riemannFit::Array2xNd = typedef Eigen::Array<double, 2, N>

Definition at line 31 of file FitUtils.h.

ArrayNd

template<int N>

using riemannFit::ArrayNd = typedef Eigen::Array<double, N, N>

Definition at line 23 of file FitUtils.h.

Map3x4d

typedef Eigen::Map< Matrix3x4d, 0, Eigen::Stride< 3 *stride, stride > > riemannFit::Map3x4d

Definition at line 24 of file HelixFit.h.

Map3xNd

template<int N>

using riemannFit::Map3xNd = typedef Eigen::Map<Matrix3xNd<N>, 0, Eigen::Stride<3 * stride, stride> >

Definition at line 32 of file HelixFit.h.

Map4d

typedef Eigen::Map< Vector4d, 0, Eigen::InnerStride< stride > > riemannFit::Map4d

Definition at line 39 of file HelixFit.h.

Map6x4f

typedef Eigen::Map< Matrix6x4f, 0, Eigen::Stride< 6 *stride, stride > > riemannFit::Map6x4f

Definition at line 26 of file HelixFit.h.

Map6xNf

template<int N>

using riemannFit::Map6xNf = typedef Eigen::Map<Matrix6xNf<N>, 0, Eigen::Stride<6 * stride, stride> >

Definition at line 37 of file HelixFit.h.

Matrix2d

typedef Eigen::Matrix2d riemannFit::Matrix2d

Definition at line 14 of file FitResult.h.

Matrix2Nd

template<int N>

using riemannFit::Matrix2Nd = typedef Eigen::Matrix<double, 2 * N, 2 * N>

Definition at line 25 of file FitUtils.h.

Matrix2x3d

typedef Eigen::Matrix< double, 2, 3 > riemannFit::Matrix2x3d

Definition at line 49 of file FitUtils.h.

Matrix2xNd

template<int N>

using riemannFit::Matrix2xNd = typedef Eigen::Matrix<double, 2, N>

Definition at line 29 of file FitUtils.h.

Matrix3d

typedef Eigen::Matrix3d riemannFit::Matrix3d

Definition at line 15 of file FitResult.h.

Matrix3f

typedef Eigen::Matrix3f riemannFit::Matrix3f

Definition at line 51 of file FitUtils.h.

Matrix3Nd

template<int N>

using riemannFit::Matrix3Nd = typedef Eigen::Matrix<double, 3 * N, 3 * N>

Definition at line 27 of file FitUtils.h.

Matrix3x4d

typedef Eigen::Matrix< double, 3, 4 > riemannFit::Matrix3x4d

Definition at line 23 of file HelixFit.h.

Matrix3xNd

template<int N>

using riemannFit::Matrix3xNd = typedef Eigen::Matrix<double, 3, N>

Definition at line 30 of file HelixFit.h.

Matrix4d

typedef Eigen::Matrix4d riemannFit::Matrix4d

Definition at line 16 of file FitResult.h.

Matrix5d

typedef Eigen::Matrix< double, 5, 5 > riemannFit::Matrix5d

Definition at line 17 of file FitResult.h.

Matrix6d

typedef Eigen::Matrix< double, 6, 6 > riemannFit::Matrix6d

Definition at line 18 of file FitResult.h.

Matrix6x4f

typedef Eigen::Matrix< float, 6, 4 > riemannFit::Matrix6x4f

Definition at line 25 of file HelixFit.h.

Matrix6xNf

template<int N>

using riemannFit::Matrix6xNf = typedef Eigen::Matrix<float, 6, N>

Definition at line 35 of file HelixFit.h.

MatrixNd

template<int N>

using riemannFit::MatrixNd = typedef Eigen::Matrix<double, N, N>

Definition at line 19 of file FitUtils.h.

MatrixNplusONEd

template<int N>

using riemannFit::MatrixNplusONEd = typedef Eigen::Matrix<double, N + 1, N + 1>

Definition at line 21 of file FitUtils.h.

MatrixNx3d

template<int N>

using riemannFit::MatrixNx3d = typedef Eigen::Matrix<double, N, 3>

Definition at line 33 of file FitUtils.h.

MatrixNx5d

template<int N>

using riemannFit::MatrixNx5d = typedef Eigen::Matrix<double, N, 5>

Definition at line 35 of file FitUtils.h.

MatrixXd

typedef Eigen::MatrixXd riemannFit::MatrixXd

Definition at line 17 of file FitUtils.h.

RowVector2Nd

template<int N>

using riemannFit::RowVector2Nd = typedef Eigen::Matrix<double, 1, 2 * N>

Definition at line 47 of file FitUtils.h.

RowVectorNd

template<int N>

using riemannFit::RowVectorNd = typedef Eigen::Matrix<double, 1, 1, N>

Definition at line 45 of file FitUtils.h.

Vector2d

typedef Eigen::Vector2d riemannFit::Vector2d

Definition at line 10 of file FitResult.h.

Vector2Nd

template<int N>

using riemannFit::Vector2Nd = typedef Eigen::Matrix<double, 2 * N, 1>

Definition at line 41 of file FitUtils.h.

Vector3d

typedef Eigen::Vector3d riemannFit::Vector3d

Definition at line 11 of file FitResult.h.

Vector3f

typedef Eigen::Vector3f riemannFit::Vector3f

Definition at line 52 of file FitUtils.h.

Vector3Nd

template<int N>

using riemannFit::Vector3Nd = typedef Eigen::Matrix<double, 3 * N, 1>

Definition at line 43 of file FitUtils.h.

Vector4d

typedef Eigen::Vector4d riemannFit::Vector4d

Definition at line 12 of file FitResult.h.

Vector4f

typedef Eigen::Vector4f riemannFit::Vector4f

Definition at line 53 of file FitUtils.h.

Vector5d

typedef Eigen::Matrix< double, 5, 1 > riemannFit::Vector5d

Definition at line 13 of file FitResult.h.

Vector6f

typedef Eigen::Matrix< double, 6, 1 > riemannFit::Vector6f

Definition at line 54 of file FitUtils.h.

VectorNd

template<int N>

using riemannFit::VectorNd = typedef Eigen::Matrix<double, N, 1>

Definition at line 37 of file FitUtils.h.

VectorNplusONEd

template<int N>

using riemannFit::VectorNplusONEd = typedef Eigen::Matrix<double, N + 1, 1>

Definition at line 39 of file FitUtils.h.

VectorXd

typedef Eigen::VectorXd riemannFit::VectorXd

Definition at line 16 of file FitUtils.h.

Function Documentation

charge()

template<typename M2xN >

__host__ __device__ int32_t riemannFit::charge ( const M2xN &  p2D, const Vector3d &  par_uvr  )

inline

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 291 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 M2xN , typename V4 , int N>

__host__ __device__ CircleFit riemannFit::circleFit ( const M2xN &  hits2D, const Matrix2Nd< N > &  hits_cov2D, const V4 &  fast_fit, const VectorNd< N > &  rad, const double  bField, const bool  error  )

inline

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 458 of file RiemannFit.h.

References a, funct::abs(), b, Calorimetry_cff::bField, charge(), riemannFit::CircleFit::chi2, nano_mu_local_reco_cff::chi2, ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), riemannFit::CircleFit::cov, cov_carttorad_prefit(), cov_radtocart(), dumpMFGeometry_cfg::delta, 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, riemannFit::CircleFit::par, printIt(), riemannFit::CircleFit::qCharge, 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 helixFit().

                                                                   {
#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(hits2D, vMat, fast_fit, rad).asDiagonal();
      MatrixNd<N> scatterCovRadMat = scatter_cov_rad(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(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(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(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(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(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 VNd1 , typename VNd2 >

__host__ __device__ void riemannFit::computeRadLenUniformMaterial ( const VNd1 &  length_values, VNd2 &  rad_lengths  )

inline

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 34 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) = std::abs(length_values(j) - length_values(j - 1)) * xx_0_inv;
    }
  }

cov_carttorad()

template<typename M2xN , int N>

__host__ __device__ VectorNd<N> riemannFit::cov_carttorad ( const M2xN &  p2D, const Matrix2Nd< N > &  cov_cart, const VectorNd< N > &  rad  )

inline

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 211 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 M2xN , typename V4 , int N>

__host__ __device__ VectorNd<N> riemannFit::cov_carttorad_prefit ( const M2xN &  p2D, const Matrix2Nd< N > &  cov_cart, V4 &  fast_fit, const VectorNd< N > &  rad  )

inline

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 242 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(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 M2xN , int N>

__host__ __device__ Matrix2Nd<N> riemannFit::cov_radtocart ( const M2xN &  p2D, const MatrixNd< N > &  cov_rad, const VectorNd< N > &  rad  )

inline

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 174 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()

__host__ __device__ double riemannFit::cross2D ( const Vector2d &  a, const Vector2d &  b  )

inline

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 79 of file FitUtils.h.

References a, and b.

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

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

fastFit()

template<typename M3xN , typename V4 >

__host__ __device__ void riemannFit::fastFit ( const M3xN &  hits, V4 &  result  )

inline

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 378 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_2023v12_cff::flip, hfClusterShapes_cfi::hits, N, dqmiodumpmetadata::n, printIt(), mps_fire::result, sqr(), and mathSSE::sqrt().

Referenced by for(), and helixFit().

                                                                        {
    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(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()

__host__ __device__ void riemannFit::fromCircleToPerigee ( CircleFit &  circle )

inline

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 196 of file FitUtils.h.

References riemannFit::CircleFit::cov, riemannFit::CircleFit::par, riemannFit::CircleFit::qCharge, sqr(), and mathSSE::sqrt().

                                                                         {
    Vector3d par_pak;
    const double temp0 = circle.par.head(2).squaredNorm();
    const double temp1 = sqrt(temp0);
    par_pak << atan2(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;
  }

helixFit()

template<int N>

HelixFit riemannFit::helixFit ( const Matrix3xNd< N > &  hits, const Eigen::Matrix< float, 6, N > &  hits_ge, const double  bField, const bool  error  )

inline

Helix fit by three step: -fast pre-fit (see Fast_fit() for further info);
-circle fit of hits projected in the transverse plane by Riemann-Chernov algorithm (see Circle_fit() for further info);
-line fit of hits projected on cylinder surface by orthogonal distance regression (see Line_fit for further info).
Points must be passed ordered (from inner to outer layer).

Parameters

hits	Matrix3xNd hits coordinates in this form: |x0|x1|x2|...|xn| |y0|y1|y2|...|yn| |z0|z1|z2|...|zn|
hits_cov	Matrix3Nd covariance matrix in this form (()->cov()): |(x0,x0)|(x1,x0)|(x2,x0)|.|(y0,x0)|(y1,x0)|(y2,x0)|.|(z0,x0)|(z1,x0)|(z2,x0)| |(x0,x1)|(x1,x1)|(x2,x1)|.|(y0,x1)|(y1,x1)|(y2,x1)|.|(z0,x1)|(z1,x1)|(z2,x1)| |(x0,x2)|(x1,x2)|(x2,x2)|.|(y0,x2)|(y1,x2)|(y2,x2)|.|(z0,x2)|(z1,x2)|(z2,x2)| . . . . . . . . . . . |(x0,y0)|(x1,y0)|(x2,y0)|.|(y0,y0)|(y1,y0)|(y2,x0)|.|(z0,y0)|(z1,y0)|(z2,y0)| |(x0,y1)|(x1,y1)|(x2,y1)|.|(y0,y1)|(y1,y1)|(y2,x1)|.|(z0,y1)|(z1,y1)|(z2,y1)| |(x0,y2)|(x1,y2)|(x2,y2)|.|(y0,y2)|(y1,y2)|(y2,x2)|.|(z0,y2)|(z1,y2)|(z2,y2)| . . . . . . . . . . . |(x0,z0)|(x1,z0)|(x2,z0)|.|(y0,z0)|(y1,z0)|(y2,z0)|.|(z0,z0)|(z1,z0)|(z2,z0)| |(x0,z1)|(x1,z1)|(x2,z1)|.|(y0,z1)|(y1,z1)|(y2,z1)|.|(z0,z1)|(z1,z1)|(z2,z1)| |(x0,z2)|(x1,z2)|(x2,z2)|.|(y0,z2)|(y1,z2)|(y2,z2)|.|(z0,z2)|(z1,z2)|(z2,z2)|
bField	magnetic field in the center of the detector in Gev/cm/c unit, in order to perform pt calculation.
error	flag for error computation.
scattering	flag for multiple scattering treatment. (see Circle_fit() documentation for further info).

Warning

see Circle_fit(), Line_fit() and Fast_fit() warnings.

Definition at line 978 of file RiemannFit.h.

References Calorimetry_cff::bField, riemannFit::CircleFit::chi2, riemannFit::HelixFit::chi2_circle, riemannFit::HelixFit::chi2_line, circleFit(), ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), riemannFit::CircleFit::cov, riemannFit::HelixFit::cov, relativeConstraints::error, ALPAKA_ACCELERATOR_NAMESPACE::brokenline::fast_fit, fastFit(), hfClusterShapes_cfi::hits, mps_splice::line, lineFit(), loadCovariance2D(), N, dqmiodumpmetadata::n, riemannFit::CircleFit::par, riemannFit::HelixFit::par, par_uvrtopak(), riemannFit::CircleFit::qCharge, riemannFit::HelixFit::qCharge, and parallelization::uint.

Referenced by PixelNtupletsFitter::run().

                                             {
    constexpr uint n = N;
    VectorNd<4> rad = (hits.block(0, 0, 2, n).colwise().norm());

    // Fast_fit gives back (X0, Y0, R, theta) w/o errors, using only 3 points.
    Vector4d fast_fit;
    fastFit(hits, fast_fit);
    riemannFit::Matrix2Nd<N> hits_cov = MatrixXd::Zero(2 * n, 2 * n);
    riemannFit::loadCovariance2D(hits_ge, hits_cov);
    CircleFit circle = circleFit(hits.block(0, 0, 2, n), hits_cov, fast_fit, rad, bField, error);
    LineFit line = lineFit(hits, hits_ge, circle, fast_fit, bField, error);

    par_uvrtopak(circle, bField, error);

    HelixFit helix;
    helix.par << circle.par, line.par;
    if (error) {
      helix.cov = 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;

    return helix;
  }

lineFit()

template<typename M3xN , typename M6xN , typename V4 >

__host__ __device__ LineFit riemannFit::lineFit ( const M3xN &  hits, const M6xN &  hits_ge, const CircleFit &  circle, const V4 &  fast_fit, const double  bField, const bool  error  )

inline

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 788 of file RiemannFit.h.

References Calorimetry_cff::bField, nano_mu_local_reco_cff::chi2, ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), funct::cos(), riemannFit::CircleFit::cov, 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, riemannFit::CircleFit::par, printIt(), riemannFit::CircleFit::qCharge, makeMuonMisalignmentScenario::rot, scatterCovLine(), funct::sin(), mkfit::Const::sol, sqr(), theta(), createJobs::tmp, and parallelization::uint.

Referenced by helixFit().

                                                               {
    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(-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(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 M6xNf , typename M3xNd >

__host__ __device__ void riemannFit::loadCovariance	(	M6xNf const &	ge,
		M3xNd &	hits_cov
	)

Definition at line 120 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 M6xNf , typename M2Nd >

__host__ __device__ void riemannFit::loadCovariance2D	(	M6xNf const &	ge,
		M2Nd &	hits_cov
	)

load error in CMSSW format to our formalism

Definition at line 88 of file FitUtils.h.

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

Referenced by helixFit().

                                                                             {
    // 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()

__host__ __device__ Vector2d riemannFit::min_eigen2D ( const Matrix2d &  aMat, double &  chi2  )

inline

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 356 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()

__host__ __device__ Vector3d riemannFit::min_eigen3D ( const Matrix3d &  A, double &  chi2  )

inline

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 313 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()

__host__ __device__ Vector3d riemannFit::min_eigen3D_fast ( const Matrix3d &  A )

inline

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 338 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()

__host__ __device__ void riemannFit::par_uvrtopak ( CircleFit &  circle, const double  B, const bool  error  )

inline

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 173 of file FitUtils.h.

References B, riemannFit::CircleFit::cov, relativeConstraints::error, riemannFit::CircleFit::par, riemannFit::CircleFit::qCharge, sqr(), and mathSSE::sqrt().

Referenced by helixFit().

                                                                                                    {
    Vector3d par_pak;
    const double temp0 = circle.par.head(2).squaredNorm();
    const double temp1 = sqrt(temp0);
    par_pak << atan2(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<class C >

__host__ __device__ void riemannFit::printIt	(	C *	m,
		const char *	prefix = ""
	)

Definition at line 53 of file FitUtils.h.

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

Referenced by circleFit(), cov_radtocart(), fastFit(), lineFit(), scatter_cov_rad(), and scatterCovLine().

                                                                  {
#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
  }

rolling_fits()

template<auto Start, auto End, auto Inc, class F >

constexpr void riemannFit::rolling_fits

(

F &&

f

)

Definition at line 42 of file HelixFit.h.

References ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), PixelRegions::End, and f.

                                     {
    if constexpr (Start < End) {
      f(std::integral_constant<decltype(Start), Start>());
      rolling_fits<Start + Inc, End, Inc>(f);
    }
  }

scatter_cov_rad()

template<typename M2xN , typename V4 , int N>

__host__ __device__ MatrixNd<N> riemannFit::scatter_cov_rad ( const M2xN &  p2D, const V4 &  fast_fit, VectorNd< N > const &  rad, double  B  )

inline

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 126 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 = std::min(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(-oVec, pVec);
      const double dot = (-oVec).dot(pVec);
      const double tempAtan2 = atan2(cross, dot);
      s_values(i) = std::abs(tempAtan2 * fast_fit(2));
    }
    computeRadLenUniformMaterial(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 < std::min(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 V4 , typename VNd1 , typename VNd2 , int N>

__host__ __device__ auto riemannFit::scatterCovLine ( 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  )

inline

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 65 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 = std::min(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(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 < std::min(k, l); ++i) {
          tmp(k + n, l + n) += std::abs(s_values(k) - s_values(i)) * std::abs(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 riemannFit::sqr

(

const T

a

)

raise to square.

Definition at line 67 of file FitUtils.h.

References a.

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

                             {
    return a * a;
  }

transformToPerigeePlane() [1/2]

template<typename VI5 , typename MI5 , typename VO5 , typename MO5 >

void riemannFit::transformToPerigeePlane ( VI5 const &  ip, MI5 const &  icov, VO5 &  op, MO5 &  ocov  )

inline

Definition at line 60 of file FitUtils.h.

References findAndChange::op, and mathSSE::sqrt().

Referenced by L2TauNNProducerAlpaka::impactParameter(), L2TauNNProducer::impactParameter(), SeedProducerFromSoAT< TrackerTraits >::produce(), PixelTrackProducerFromSoAAlpaka< TrackerTraits >::produce(), and PixelTrackProducerFromSoAT< TrackerTraits >::produce().

                                                                                          {
    auto sinTheta2 = 1. / (1. + ip(3) * ip(3));
    auto sinTheta = std::sqrt(sinTheta2);
    auto cosTheta = ip(3) * sinTheta;

    op(0) = sinTheta * ip(2);
    op(1) = 0.;
    op(2) = -ip(3);
    op(3) = ip(1);
    op(4) = -ip(4);

    Matrix5d jMat = Matrix5d::Zero();

    jMat(0, 2) = sinTheta;
    jMat(0, 3) = -sinTheta2 * cosTheta * ip(2);
    jMat(1, 0) = 1.;
    jMat(2, 3) = -1.;
    jMat(3, 1) = 1.;
    jMat(4, 4) = -1;

    ocov = jMat * icov * jMat.transpose();
  }

transformToPerigeePlane() [2/2]

template<typename VI5 , typename MI5 , typename VO5 , typename MO5 >

__host__ __device__ void riemannFit::transformToPerigeePlane ( VI5 const &  ip, MI5 const &  icov, VO5 &  op, MO5 &  ocov  )

inline

Definition at line 218 of file FitUtils.h.

References findAndChange::op, and mathSSE::sqrt().

                                                                                                              {
    auto sinTheta2 = 1. / (1. + ip(3) * ip(3));
    auto sinTheta = std::sqrt(sinTheta2);
    auto cosTheta = ip(3) * sinTheta;

    op(0) = sinTheta * ip(2);
    op(1) = 0.;
    op(2) = -ip(3);
    op(3) = ip(1);
    op(4) = -ip(4);

    Matrix5d jMat = Matrix5d::Zero();

    jMat(0, 2) = sinTheta;
    jMat(0, 3) = -sinTheta2 * cosTheta * ip(2);
    jMat(1, 0) = 1.;
    jMat(2, 3) = -1.;
    jMat(3, 1) = 1.;
    jMat(4, 4) = -1;

    ocov = jMat * icov * jMat.transpose();
  }

weightCircle()

template<int N>

__host__ __device__ VectorNd<N> riemannFit::weightCircle ( const MatrixNd< N > &  cov_rad_inv )

inline

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 278 of file RiemannFit.h.

Referenced by circleFit().

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

Variable Documentation

epsilon

constexpr double riemannFit::epsilon = 1.e-4

used in numerical derivative (J2 in Circle_fit())

Definition at line 14 of file FitUtils.h.

Referenced by circleFit().

maxNumberOfConcurrentFits

constexpr uint32_t riemannFit::maxNumberOfConcurrentFits = 32 * 1024

Definition at line 21 of file HelixFit.h.

Referenced by ALPAKA_ACCELERATOR_NAMESPACE::Kernel_BLFastFit< N, TrackerTraits >::operator()(), ALPAKA_ACCELERATOR_NAMESPACE::Kernel_FastFit< N, TrackerTraits >::operator()(), ALPAKA_ACCELERATOR_NAMESPACE::Kernel_CircleFit< N, TrackerTraits >::operator()(), ALPAKA_ACCELERATOR_NAMESPACE::Kernel_LineFit< N, TrackerTraits >::operator()(), and ALPAKA_ACCELERATOR_NAMESPACE::Kernel_BLFit< N, TrackerTraits >::operator()().

stride

constexpr uint32_t riemannFit::stride = maxNumberOfConcurrentFits

Definition at line 22 of file HelixFit.h.

Referenced by caHitNtupletGeneratorKernels::__attribute__(), HcalDigisProducerGPU::acquire(), cms::rocmtest::add_vectors_d(), cms::rocmtest::add_vectors_f(), PhysicsPerformanceDBWriterFromFile_WPandPayload::beginJob(), PhysicsPerformanceDBWriterFromFile_WPandPayload_IOV::beginJob(), ALPAKA_ACCELERATOR_NAMESPACE::CAHitNtupletGeneratorKernels< TTTraits >::buildDoublets(), cms::soa::SoAParametersImpl< SoAColumnType::eigen, T >::checkAlignment(), calo::multifit::compute_decomposition_unrolled(), cms::alpakatools::uniform_groups_along< TAcc, Dim, typename >::const_iterator::const_iterator(), cms::alpakatools::independent_groups_along< TAcc, Dim, typename >::const_iterator::const_iterator(), npstat::ArrayND< Numeric >::contract(), npstat::ArrayND< Numeric >::contractLoop(), npstat::ArrayND< Numeric >::convertToLastDimCdfLoop(), npstat::ArrayND< Numeric >::functorFillLoop(), ALPAKA_ACCELERATOR_NAMESPACE::hcalFastCluster_exotic(), ALPAKA_ACCELERATOR_NAMESPACE::hcalFastCluster_multiSeedIterative(), ALPAKA_ACCELERATOR_NAMESPACE::hcalFastCluster_multiSeedParallel(), npstat::ArrayND< Numeric >::jointSliceLoop(), ALPAKA_ACCELERATOR_NAMESPACE::CAHitNtupletGeneratorKernels< TTTraits >::launchKernels(), npstat::ArrayND< Numeric >::linearFillLoop(), npstat::ArrayND< Numeric >::makeUnit(), npstat::ArrayND< Numeric >::processSubrangeLoop(), npstat::ArrayND< Numeric >::projectInnerLoop(), npstat::ArrayND< Numeric >::projectInnerLoop2(), npstat::ArrayND< Numeric >::projectLoop(), npstat::ArrayND< Numeric >::projectLoop2(), npstat::ArrayND< Numeric >::scaleBySliceInnerLoop(), npstat::ArrayND< Numeric >::scaleBySliceLoop(), l1tVertexFinder::VertexFinder::strided_iota(), npstat::ArrayND< Numeric >::sumBelowLoop(), and ALPAKA_ACCELERATOR_NAMESPACE::testESAlgoAsync().