Classes
class	helixFit
	Helix fit by three step: -fast pre-fit (see Fast_fit() for further info); -circle fit of the hits projected in the transverse plane by Broken Line algorithm (see BL_Circle_fit() for further info); -line fit of the hits projected on the (pre-fitted) cilinder surface by Broken Line algorithm (see BL_Line_fit() for further info); Points must be passed ordered (from inner to outer layer). More...

struct	PreparedBrokenLineData
	data needed for the Broken Line fit procedure. More...

Typedefs
using	karimaki_circle_fit = riemannFit::CircleFit
	< Karimäki's parameters: (phi, d, k=1/R) More...

Functions
template<typename TAcc , typename M3xN , typename V4 , int n>
ALPAKA_FN_ACC ALPAKA_FN_INLINE void	__attribute__ ((always_inline)) prepareBrokenLineData(const TAcc &acc
	Computes the data needed for the Broken Line fit procedure that are mainly common for the circle and the line fit. More...

template<typename TAcc , typename M3xN , typename M6xN , typename V4 , int n>
ALPAKA_FN_ACC ALPAKA_FN_INLINE void	circleFit (const TAcc &acc, const M3xN &hits, const M6xN &hits_ge, const V4 &fast_fit, const double bField, PreparedBrokenLineData< n > &data, karimaki_circle_fit &circle_results)
	Performs the Broken Line fit in the curved track case (that is, the fit parameters are the interceptions u and the curvature correction ). More...

if	constexpr (n > 3)

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. More...

	for (u_int i=0;i< n;i++)

template<typename TAcc , typename V4 , typename M6xN , int n>
ALPAKA_FN_ACC ALPAKA_FN_INLINE void	lineFit (const TAcc &acc, const M6xN &hits_ge, const V4 &fast_fit, const double bField, const PreparedBrokenLineData< n > &data, riemannFit::LineFit &line_results)
	Performs the Broken Line fit in the straight track case (that is, the fit parameters are only the interceptions u). More...

template<typename TAcc , int n>
ALPAKA_FN_ACC ALPAKA_FN_INLINE riemannFit::MatrixNd< n >	matrixC_u (const TAcc &acc, const riemannFit::VectorNd< n > &weights, const riemannFit::VectorNd< n > &sTotal, const riemannFit::VectorNd< n > &varBeta)
	Computes the n-by-n band matrix obtained minimizing the Broken Line's cost function w.r.t u. This is the whole matrix in the case of the line fit and the main n-by-n block in the case of the circle fit. More...

template<typename TAcc >
ALPAKA_FN_ACC ALPAKA_FN_INLINE double	multScatt (const TAcc &acc, const double &length, const double bField, const double radius, int layer, double slope)
	Computes the Coulomb multiple scattering variance of the planar angle. More...

template<typename TAcc >
ALPAKA_FN_ACC ALPAKA_FN_INLINE riemannFit::Matrix2d	rotationMatrix (const TAcc &acc, double slope)
	Computes the 2D rotation matrix that transforms the line y=slope*x into the line y=0. More...

template<typename TAcc >
ALPAKA_FN_ACC ALPAKA_FN_INLINE void	translateKarimaki (const TAcc &acc, karimaki_circle_fit &circle, double x0, double y0, riemannFit::Matrix3d &jacobian)
	Changes the Karimäki parameters (and consequently their covariance matrix) under a translation of the coordinate system, such that the old origin has coordinates (x0,y0) in the new coordinate system. The formulas are taken from Karimäki V., 1990, Effective circle fitting for particle trajectories, Nucl. Instr. and Meth. A305 (1991) 187. More...

results	varBeta (0)

Variables
ALPAKA_FN_ACC ALPAKA_FN_INLINE void const M3xN const V4 const double	bField

	dVec = hits.block(0, mId, 2, 1) - hits.block(0, 0, 2, 1)

riemannFit::Vector2d	eVec = hits.block(0, n - 1, 2, 1) - hits.block(0, mId, 2, 1)

ALPAKA_FN_ACC ALPAKA_FN_INLINE void const M3xN const V4 &	fast_fit

ALPAKA_FN_ACC ALPAKA_FN_INLINE void const M3xN &	hits

int	mId = 1

riemannFit::Matrix2xNd< n >	pointsSZ = riemannFit::Matrix2xNd<n>::Zero()

results	qCharge = riemannFit::cross2D(acc, dVec, eVec) > 0 ? -1 : 1

results	radii = hits.block(0, 0, 2, n) - fast_fit.head(2) * riemannFit::MatrixXd::Constant(1, n, 1)

ALPAKA_FN_ACC ALPAKA_FN_INLINE void const M3xN const V4 const double PreparedBrokenLineData< n > &	results

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

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

results	sTotal = pointsSZ.block(0, 0, 1, n).transpose()

results	zInSZplane = pointsSZ.block(1, 0, 1, n).transpose()

riemannFit::VectorNd< n >	zVec = hits.block(2, 0, 1, n).transpose()

Typedef Documentation

◆ karimaki_circle_fit

using ALPAKA_ACCELERATOR_NAMESPACE::brokenline::karimaki_circle_fit = typedef riemannFit::CircleFit

< Karimäki's parameters: (phi, d, k=1/R)

Definition at line 24 of file BrokenLine.h.

Function Documentation

◆ attribute()

template<typename TAcc , typename M3xN , typename V4 , int n>

ALPAKA_FN_ACC ALPAKA_FN_INLINE void ALPAKA_ACCELERATOR_NAMESPACE::brokenline::__attribute__ ( (always_inline) ) const &

Computes the data needed for the Broken Line fit procedure that are mainly common for the circle and the line fit.

Parameters

hits	hits coordinates.
fast_fit	pre-fit result in the form (X0,Y0,R,tan(theta)).
bField	magnetic field in Gev/cm/c.
results	PreparedBrokenLineData to be filled (see description of PreparedBrokenLineData).

◆ circleFit()

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

ALPAKA_FN_ACC ALPAKA_FN_INLINE void ALPAKA_ACCELERATOR_NAMESPACE::brokenline::circleFit	(	const TAcc &	acc,
		const M3xN &	hits,
		const M6xN &	hits_ge,
		const V4 &	fast_fit,
		const double	bField,
		PreparedBrokenLineData< n > &	data,
		karimaki_circle_fit &	circle_results
	)

Performs the Broken Line fit in the curved track case (that is, the fit parameters are the interceptions u and the curvature correction ).

Parameters

hits	hits coordinates.
hits_cov	hits covariance matrix.
fast_fit	pre-fit result in the form (X0,Y0,R,tan(theta)).
bField	magnetic field in Gev/cm/c.
data	PreparedBrokenLineData.
circle_results	struct to be filled with the results in this form: -par parameter of the line in this form: (phi, d, k); -cov covariance matrix of the fitted parameter; -chi2 value of the cost function in the minimum.

The function implements the steps 2 and 3 of the Broken Line fit with the curvature correction.
The step 2 is the least square fit, done by imposing the minimum constraint on the cost function and solving the consequent linear system. It determines the fitted parameters u and and their covariance matrix. The step 3 is the correction of the fast pre-fitted parameters for the innermost part of the track. It is first done in a comfortable coordinate system (the one in which the first hit is the origin) and then the parameters and their covariance matrix are transformed to the original coordinate system.

Definition at line 323 of file BrokenLine.h.

References bField, riemannFit::CircleFit::chi2, gather_cfg::cout, riemannFit::CircleFit::cov, data, dVec, eVec, fast_fit, hits, mps_fire::i, l1tstage2_dqm_sourceclient-live_cfg::invert, matrixC_u(), multScatt(), dqmiodumpmetadata::n, riemannFit::CircleFit::par, riemannFit::CircleFit::qCharge, radii, rotationMatrix(), rotMat, slope, ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::sqr(), mathSSE::sqrt(), sTotal, translateKarimaki(), varBeta(), and zInSZplane.

Referenced by ALPAKA_ACCELERATOR_NAMESPACE::Kernel_BLFit< N, TrackerTraits >::operator()(), and ALPAKA_ACCELERATOR_NAMESPACE::brokenline::helixFit< n >::operator()().

                                                                                      {
     circle_results.qCharge = data.qCharge;
     auto& radii = data.radii;
     const auto& sTransverse = data.sTransverse;
     const auto& sTotal = data.sTotal;
     auto& zInSZplane = data.zInSZplane;
     auto& varBeta = data.varBeta;
     const double slope = -circle_results.qCharge / fast_fit(3);
     varBeta *= 1. + riemannFit::sqr(slope);  // the kink angles are projected!
 
     for (u_int i = 0; i < n; i++) {
       zInSZplane(i) = radii.block(0, i, 2, 1).norm() - fast_fit(2);
     }
 
     riemannFit::Matrix2d vMat;           // covariance matrix
     riemannFit::VectorNd<n> weightsVec;  // weights
     riemannFit::Matrix2d rotMat;         // rotation matrix point by point
     for (u_int i = 0; i < n; i++) {
       vMat(0, 0) = hits_ge.col(i)[0];               // x errors
       vMat(0, 1) = vMat(1, 0) = hits_ge.col(i)[1];  // cov_xy
       vMat(1, 1) = hits_ge.col(i)[2];               // y errors
       rotMat = rotationMatrix(acc, -radii(0, i) / radii(1, i));
       weightsVec(i) =
           1. / ((rotMat * vMat * rotMat.transpose())(1, 1));  // compute the orthogonal weight point by point
     }
 
     riemannFit::VectorNplusONEd<n> r_uVec;
     r_uVec(n) = 0;
     for (u_int i = 0; i < n; i++) {
       r_uVec(i) = weightsVec(i) * zInSZplane(i);
     }
 
     riemannFit::MatrixNplusONEd<n> c_uMat;
     c_uMat.block(0, 0, n, n) = matrixC_u(acc, weightsVec, sTransverse, varBeta);
     c_uMat(n, n) = 0;
     //add the border to the c_uMat matrix
     for (u_int i = 0; i < n; i++) {
       c_uMat(i, n) = 0;
       if (i > 0 && i < n - 1) {
         c_uMat(i, n) +=
             -(sTransverse(i + 1) - sTransverse(i - 1)) * (sTransverse(i + 1) - sTransverse(i - 1)) /
             (2. * varBeta(i) * (sTransverse(i + 1) - sTransverse(i)) * (sTransverse(i) - sTransverse(i - 1)));
       }
       if (i > 1) {
         c_uMat(i, n) +=
             (sTransverse(i) - sTransverse(i - 2)) / (2. * varBeta(i - 1) * (sTransverse(i) - sTransverse(i - 1)));
       }
       if (i < n - 2) {
         c_uMat(i, n) +=
             (sTransverse(i + 2) - sTransverse(i)) / (2. * varBeta(i + 1) * (sTransverse(i + 1) - sTransverse(i)));
       }
       c_uMat(n, i) = c_uMat(i, n);
       if (i > 0 && i < n - 1)
         c_uMat(n, n) += riemannFit::sqr(sTransverse(i + 1) - sTransverse(i - 1)) / (4. * varBeta(i));
     }
 
 #ifdef CPP_DUMP
     std::cout << "CU5\n" << c_uMat << std::endl;
 #endif
     riemannFit::MatrixNplusONEd<n> iMat;
     math::cholesky::invert(c_uMat, iMat);
 #ifdef CPP_DUMP
     std::cout << "I5\n" << iMat << std::endl;
 #endif
     riemannFit::VectorNplusONEd<n> uVec = iMat * r_uVec;  // obtain the fitted parameters by solving the linear system
 
     // compute (phi, d_ca, k) in the system in which the midpoint of the first two corrected hits is the origin...
 
     radii.block(0, 0, 2, 1) /= radii.block(0, 0, 2, 1).norm();
     radii.block(0, 1, 2, 1) /= radii.block(0, 1, 2, 1).norm();
 
     riemannFit::Vector2d dVec = hits.block(0, 0, 2, 1) + (-zInSZplane(0) + uVec(0)) * radii.block(0, 0, 2, 1);
     riemannFit::Vector2d eVec = hits.block(0, 1, 2, 1) + (-zInSZplane(1) + uVec(1)) * radii.block(0, 1, 2, 1);
     auto eMinusd = eVec - dVec;
     auto eMinusd2 = eMinusd.squaredNorm();
     auto tmp1 = 1. / eMinusd2;
     auto tmp2 = alpaka::math::sqrt(acc, riemannFit::sqr(fast_fit(2)) - 0.25 * eMinusd2);
 
     circle_results.par << atan2(eMinusd(1), eMinusd(0)), circle_results.qCharge * (tmp2 - fast_fit(2)),
         circle_results.qCharge * (1. / fast_fit(2) + uVec(n));
 
     tmp2 = 1. / tmp2;
 
     riemannFit::Matrix3d jacobian;
     jacobian << (radii(1, 0) * eMinusd(0) - eMinusd(1) * radii(0, 0)) * tmp1,
         (radii(1, 1) * eMinusd(0) - eMinusd(1) * radii(0, 1)) * tmp1, 0,
         circle_results.qCharge * (eMinusd(0) * radii(0, 0) + eMinusd(1) * radii(1, 0)) * tmp2,
         circle_results.qCharge * (eMinusd(0) * radii(0, 1) + eMinusd(1) * radii(1, 1)) * tmp2, 0, 0, 0,
         circle_results.qCharge;
 
     circle_results.cov << iMat(0, 0), iMat(0, 1), iMat(0, n), iMat(1, 0), iMat(1, 1), iMat(1, n), iMat(n, 0),
         iMat(n, 1), iMat(n, n);
 
     circle_results.cov = jacobian * circle_results.cov * jacobian.transpose();
 
     //...Translate in the system in which the first corrected hit is the origin, adding the m.s. correction...
 
     translateKarimaki(acc, circle_results, 0.5 * eMinusd(0), 0.5 * eMinusd(1), jacobian);
     circle_results.cov(0, 0) +=
         (1 + riemannFit::sqr(slope)) * multScatt(acc, sTotal(1) - sTotal(0), bField, fast_fit(2), 2, slope);
 
     //...And translate back to the original system
 
     translateKarimaki(acc, circle_results, dVec(0), dVec(1), jacobian);
 
     // compute chi2
     circle_results.chi2 = 0;
     for (u_int i = 0; i < n; i++) {
       circle_results.chi2 += weightsVec(i) * riemannFit::sqr(zInSZplane(i) - uVec(i));
       if (i > 0 && i < n - 1)
         circle_results.chi2 +=
             riemannFit::sqr(uVec(i - 1) / (sTransverse(i) - sTransverse(i - 1)) -
                             uVec(i) * (sTransverse(i + 1) - sTransverse(i - 1)) /
                                 ((sTransverse(i + 1) - sTransverse(i)) * (sTransverse(i) - sTransverse(i - 1))) +
                             uVec(i + 1) / (sTransverse(i + 1) - sTransverse(i)) +
                             (sTransverse(i + 1) - sTransverse(i - 1)) * uVec(n) / 2) /
             varBeta(i);
     }
   }

◆ constexpr()

if ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr	(	n	,
		3
	)

Definition at line 164 of file BrokenLine.h.

References d1, hits, mId, and dqmiodumpmetadata::n.

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

◆ fastFit()

template<typename TAcc , typename M3xN , typename V4 >

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

A very fast helix fit.

Parameters

hits	the measured hits.

Returns: (X0,Y0,R,tan(theta)).

Warning: sign of theta is (intentionally, for now) mistaken for negative charges.

Definition at line 266 of file BrokenLine.h.

References a, funct::abs(), b, DummyCfis::c, constexpr(), ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::cross2D(), ztail::d, d1, MillePedeFileConverter_cfg::e, hits, mId, dqmiodumpmetadata::n, mps_fire::result, mathSSE::sqrt(), and createJobs::tmp.

Referenced by ALPAKA_ACCELERATOR_NAMESPACE::Kernel_BLFastFit< N, TrackerTraits >::operator()(), and ALPAKA_ACCELERATOR_NAMESPACE::brokenline::helixFit< n >::operator()().

                                                                                              {
     constexpr uint32_t n = M3xN::ColsAtCompileTime;
 
     int mId = 1;
 
     if constexpr (n > 3) {
       riemannFit::Vector2d middle = 0.5 * (hits.block(0, n - 1, 2, 1) + hits.block(0, 0, 2, 1));
       auto d1 = (hits.block(0, n / 2, 2, 1) - middle).squaredNorm();
       auto d2 = (hits.block(0, n / 2 - 1, 2, 1) - middle).squaredNorm();
       mId = d1 < d2 ? n / 2 : n / 2 - 1;
     }
 
     const riemannFit::Vector2d a = hits.block(0, mId, 2, 1) - hits.block(0, 0, 2, 1);
     const riemannFit::Vector2d b = hits.block(0, n - 1, 2, 1) - hits.block(0, mId, 2, 1);
     const riemannFit::Vector2d c = hits.block(0, 0, 2, 1) - hits.block(0, n - 1, 2, 1);
 
     auto tmp = 0.5 / riemannFit::cross2D(acc, c, a);
     result(0) = hits(0, 0) - (a(1) * c.squaredNorm() + c(1) * a.squaredNorm()) * tmp;
     result(1) = hits(1, 0) + (a(0) * c.squaredNorm() + c(0) * a.squaredNorm()) * tmp;
     // check Wikipedia for these formulas
 
     result(2) = alpaka::math::sqrt(acc, a.squaredNorm() * b.squaredNorm() * c.squaredNorm()) /
                 (2. * alpaka::math::abs(acc, riemannFit::cross2D(acc, b, a)));
     // Using Math Olympiad's formula R=abc/(4A)
 
     const riemannFit::Vector2d d = hits.block(0, 0, 2, 1) - result.head(2);
     const riemannFit::Vector2d e = hits.block(0, n - 1, 2, 1) - result.head(2);
 
     result(3) = result(2) * atan2(riemannFit::cross2D(acc, d, e), d.dot(e)) / (hits(2, n - 1) - hits(2, 0));
     // ds/dz slope between last and first point
   }

◆ for()

ALPAKA_ACCELERATOR_NAMESPACE::brokenline::for ( )

Definition at line 182 of file BrokenLine.h.

References ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::cross2D(), dVec, eVec, fast_fit, mps_fire::i, and results.

                                   {
       dVec = results.radii.block(0, i, 2, 1);
       results.sTransverse(i) =
           results.qCharge * fast_fit(2) *
           alpaka::math::atan2(acc, riemannFit::cross2D(acc, dVec, eVec), dVec.dot(eVec));  // calculates the arc length
     }

◆ lineFit()

template<typename TAcc , typename V4 , typename M6xN , int n>

ALPAKA_FN_ACC ALPAKA_FN_INLINE void ALPAKA_ACCELERATOR_NAMESPACE::brokenline::lineFit	(	const TAcc &	acc,
		const M6xN &	hits_ge,
		const V4 &	fast_fit,
		const double	bField,
		const PreparedBrokenLineData< n > &	data,
		riemannFit::LineFit &	line_results
	)

Performs the Broken Line fit in the straight track case (that is, the fit parameters are only the interceptions u).

Parameters

hits	hits coordinates.
fast_fit	pre-fit result in the form (X0,Y0,R,tan(theta)).
bField	magnetic field in Gev/cm/c.
data	PreparedBrokenLineData.
line_results	struct to be filled with the results in this form: -par parameter of the line in this form: (cot(theta), Zip); -cov covariance matrix of the fitted parameter; -chi2 value of the cost function in the minimum.

The function implements the steps 2 and 3 of the Broken Line fit without the curvature correction.
The step 2 is the least square fit, done by imposing the minimum constraint on the cost function and solving the consequent linear system. It determines the fitted parameters u and their covariance matrix. The step 3 is the correction of the fast pre-fitted parameters for the innermost part of the track. It is first done in a comfortable coordinate system (the one in which the first hit is the origin) and then the parameters and their covariance matrix are transformed to the original coordinate system.

Definition at line 472 of file BrokenLine.h.

References bField, riemannFit::LineFit::chi2, gather_cfg::cout, riemannFit::LineFit::cov, data, fast_fit, mps_fire::i, l1tstage2_dqm_sourceclient-live_cfg::invert, matrixC_u(), multScatt(), dqmiodumpmetadata::n, riemannFit::LineFit::par, radii, rotationMatrix(), rotMat, slope, ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::sqr(), sTotal, createJobs::tmp, varBeta(), w(), hltDeepSecondaryVertexTagInfosPFPuppi_cfi::weights, and zInSZplane.

Referenced by ALPAKA_ACCELERATOR_NAMESPACE::Kernel_BLFit< N, TrackerTraits >::operator()(), and ALPAKA_ACCELERATOR_NAMESPACE::brokenline::helixFit< n >::operator()().

                                                                                {
     const auto& radii = data.radii;
     const auto& sTotal = data.sTotal;
     const auto& zInSZplane = data.zInSZplane;
     const auto& varBeta = data.varBeta;
 
     const double slope = -data.qCharge / fast_fit(3);
     riemannFit::Matrix2d rotMat = rotationMatrix(acc, slope);
 
     riemannFit::Matrix3d vMat = riemannFit::Matrix3d::Zero();  // covariance matrix XYZ
     riemannFit::Matrix2x3d jacobXYZtosZ =
         riemannFit::Matrix2x3d::Zero();  // jacobian for computation of the error on s (xyz -> sz)
     riemannFit::VectorNd<n> weights = riemannFit::VectorNd<n>::Zero();
     for (u_int i = 0; i < n; i++) {
       vMat(0, 0) = hits_ge.col(i)[0];               // x errors
       vMat(0, 1) = vMat(1, 0) = hits_ge.col(i)[1];  // cov_xy
       vMat(0, 2) = vMat(2, 0) = hits_ge.col(i)[3];  // cov_xz
       vMat(1, 1) = hits_ge.col(i)[2];               // y errors
       vMat(2, 1) = vMat(1, 2) = hits_ge.col(i)[4];  // cov_yz
       vMat(2, 2) = hits_ge.col(i)[5];               // z errors
       auto tmp = 1. / radii.block(0, i, 2, 1).norm();
       jacobXYZtosZ(0, 0) = radii(1, i) * tmp;
       jacobXYZtosZ(0, 1) = -radii(0, i) * tmp;
       jacobXYZtosZ(1, 2) = 1.;
       weights(i) = 1. / ((rotMat * jacobXYZtosZ * vMat * jacobXYZtosZ.transpose() * rotMat.transpose())(
                             1, 1));  // compute the orthogonal weight point by point
     }
 
     riemannFit::VectorNd<n> r_u;
     for (u_int i = 0; i < n; i++) {
       r_u(i) = weights(i) * zInSZplane(i);
     }
 #ifdef CPP_DUMP
     std::cout << "CU4\n" << matrixC_u(w, sTotal, varBeta) << std::endl;
 #endif
     riemannFit::MatrixNd<n> iMat;
     math::cholesky::invert(matrixC_u(acc, weights, sTotal, varBeta), iMat);
 #ifdef CPP_DUMP
     std::cout << "I4\n" << iMat << std::endl;
 #endif
 
     riemannFit::VectorNd<n> uVec = iMat * r_u;  // obtain the fitted parameters by solving the linear system
 
     // line parameters in the system in which the first hit is the origin and with axis along SZ
     line_results.par << (uVec(1) - uVec(0)) / (sTotal(1) - sTotal(0)), uVec(0);
     auto idiff = 1. / (sTotal(1) - sTotal(0));
     line_results.cov << (iMat(0, 0) - 2 * iMat(0, 1) + iMat(1, 1)) * riemannFit::sqr(idiff) +
                             multScatt(acc, sTotal(1) - sTotal(0), bField, fast_fit(2), 2, slope),
         (iMat(0, 1) - iMat(0, 0)) * idiff, (iMat(0, 1) - iMat(0, 0)) * idiff, iMat(0, 0);
 
     // translate to the original SZ system
     riemannFit::Matrix2d jacobian;
     jacobian(0, 0) = 1.;
     jacobian(0, 1) = 0;
     jacobian(1, 0) = -sTotal(0);
     jacobian(1, 1) = 1.;
     line_results.par(1) += -line_results.par(0) * sTotal(0);
     line_results.cov = jacobian * line_results.cov * jacobian.transpose();
 
     // rotate to the original sz system
     auto tmp = rotMat(0, 0) - line_results.par(0) * rotMat(0, 1);
     jacobian(1, 1) = 1. / tmp;
     jacobian(0, 0) = jacobian(1, 1) * jacobian(1, 1);
     jacobian(0, 1) = 0;
     jacobian(1, 0) = line_results.par(1) * rotMat(0, 1) * jacobian(0, 0);
     line_results.par(1) = line_results.par(1) * jacobian(1, 1);
     line_results.par(0) = (rotMat(0, 1) + line_results.par(0) * rotMat(0, 0)) * jacobian(1, 1);
     line_results.cov = jacobian * line_results.cov * jacobian.transpose();
 
     // compute chi2
     line_results.chi2 = 0;
     for (u_int i = 0; i < n; i++) {
       line_results.chi2 += weights(i) * riemannFit::sqr(zInSZplane(i) - uVec(i));
       if (i > 0 && i < n - 1)
         line_results.chi2 += riemannFit::sqr(uVec(i - 1) / (sTotal(i) - sTotal(i - 1)) -
                                              uVec(i) * (sTotal(i + 1) - sTotal(i - 1)) /
                                                  ((sTotal(i + 1) - sTotal(i)) * (sTotal(i) - sTotal(i - 1))) +
                                              uVec(i + 1) / (sTotal(i + 1) - sTotal(i))) /
                              varBeta(i);
     }
   }

◆ matrixC_u()

template<typename TAcc , int n>

ALPAKA_FN_ACC ALPAKA_FN_INLINE riemannFit::MatrixNd<n> ALPAKA_ACCELERATOR_NAMESPACE::brokenline::matrixC_u	(	const TAcc &	acc,
		const riemannFit::VectorNd< n > &	weights,
		const riemannFit::VectorNd< n > &	sTotal,
		const riemannFit::VectorNd< n > &	varBeta
	)

Computes the n-by-n band matrix obtained minimizing the Broken Line's cost function w.r.t u. This is the whole matrix in the case of the line fit and the main n-by-n block in the case of the circle fit.

Parameters

weights	weights of the first part of the cost function, the one with the measurements and not the angles ({i=1}^n w*(y_i-u_i)^2).
sTotal	total distance traveled by the particle from the pre-fitted closest approach.
varBeta	kink angles' variance.

Returns: the n-by-n matrix of the linear system

Definition at line 222 of file BrokenLine.h.

References mps_fire::i, dqmiodumpmetadata::n, ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::sqr(), sTotal, varBeta(), and hltDeepSecondaryVertexTagInfosPFPuppi_cfi::weights.

Referenced by circleFit(), and lineFit().

                                                                                                          {
     riemannFit::MatrixNd<n> c_uMat = riemannFit::MatrixNd<n>::Zero();
     for (u_int i = 0; i < n; i++) {
       c_uMat(i, i) = weights(i);
       if (i > 1)
         c_uMat(i, i) += 1. / (varBeta(i - 1) * riemannFit::sqr(sTotal(i) - sTotal(i - 1)));
       if (i > 0 && i < n - 1)
         c_uMat(i, i) +=
             (1. / varBeta(i)) * riemannFit::sqr((sTotal(i + 1) - sTotal(i - 1)) /
                                                 ((sTotal(i + 1) - sTotal(i)) * (sTotal(i) - sTotal(i - 1))));
       if (i < n - 2)
         c_uMat(i, i) += 1. / (varBeta(i + 1) * riemannFit::sqr(sTotal(i + 1) - sTotal(i)));
 
       if (i > 0 && i < n - 1)
         c_uMat(i, i + 1) =
             1. / (varBeta(i) * (sTotal(i + 1) - sTotal(i))) *
             (-(sTotal(i + 1) - sTotal(i - 1)) / ((sTotal(i + 1) - sTotal(i)) * (sTotal(i) - sTotal(i - 1))));
       if (i < n - 2)
         c_uMat(i, i + 1) +=
             1. / (varBeta(i + 1) * (sTotal(i + 1) - sTotal(i))) *
             (-(sTotal(i + 2) - sTotal(i)) / ((sTotal(i + 2) - sTotal(i + 1)) * (sTotal(i + 1) - sTotal(i))));
 
       if (i < n - 2)
         c_uMat(i, i + 2) = 1. / (varBeta(i + 1) * (sTotal(i + 2) - sTotal(i + 1)) * (sTotal(i + 1) - sTotal(i)));
 
       c_uMat(i, i) *= 0.5;
     }
     return c_uMat + c_uMat.transpose();
   }

◆ multScatt()

template<typename TAcc >

ALPAKA_FN_ACC ALPAKA_FN_INLINE double ALPAKA_ACCELERATOR_NAMESPACE::brokenline::multScatt	(	const TAcc &	acc,
		const double &	length,
		const double	bField,
		const double	radius,
		int	layer,
		double	slope
	)

Computes the Coulomb multiple scattering variance of the planar angle.

Parameters

length	length of the track in the material.
bField	magnetic field in Gev/cm/c.
radius	radius of curvature (needed to evaluate p).
layer	denotes which of the four layers of the detector is the endpoint of the multiple scattered track. For example, if Layer=3, then the particle has just gone through the material between the second and the third layer.

Warning: the formula used here assumes beta=1, and so neglects the dependence of theta_0 on the mass of the particle at fixed momentum.

Returns: the variance of the planar angle ((theta_0)^2 /3).

< inverse of radiation length of the material in cm

number between 1/3 (uniform material) and 1 (thin scatterer) to be manually tuned

Definition at line 60 of file BrokenLine.h.

References funct::abs(), bField, constexpr(), fact, dqm-mbProfile::log, SiStripPI::min, HLT_2024v14_cff::pt2, CosmicsPD_Skims::radius, slope, and ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::sqr().

Referenced by circleFit(), and lineFit().

                                                                                                                 {
     // limit R to 20GeV...
     auto pt2 = alpaka::math::min(acc, 20., bField * radius);
     pt2 *= pt2;
     constexpr double inv_X0 = 0.06 / 16.;  
     //if(Layer==1) XXI_0=0.06/16.;
     // else XXI_0=0.06/16.;
     //XX_0*=1;
 
     constexpr double geometry_factor = 0.7;
     constexpr double fact = geometry_factor * riemannFit::sqr(13.6 / 1000.);
     return fact / (pt2 * (1. + riemannFit::sqr(slope))) * (alpaka::math::abs(acc, length) * inv_X0) *
            riemannFit::sqr(1. + 0.038 * log(alpaka::math::abs(acc, length) * inv_X0));
   }

◆ rotationMatrix()

template<typename TAcc >

ALPAKA_FN_ACC ALPAKA_FN_INLINE riemannFit::Matrix2d ALPAKA_ACCELERATOR_NAMESPACE::brokenline::rotationMatrix	(	const TAcc &	acc,
		double	slope
	)

Computes the 2D rotation matrix that transforms the line y=slope*x into the line y=0.

Parameters

slope tangent of the angle of rotation.

Returns: 2D rotation matrix.

Definition at line 85 of file BrokenLine.h.

References makeMuonMisalignmentScenario::rot, slope, ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::sqr(), and mathSSE::sqrt().

Referenced by circleFit(), lineFit(), and reco::BeamSpot::rotatedCovariance3D().

                                                                                                   {
     riemannFit::Matrix2d rot;
     rot(0, 0) = 1. / alpaka::math::sqrt(acc, 1. + riemannFit::sqr(slope));
     rot(0, 1) = slope * rot(0, 0);
     rot(1, 0) = -rot(0, 1);
     rot(1, 1) = rot(0, 0);
     return rot;
   }

◆ translateKarimaki()

template<typename TAcc >

ALPAKA_FN_ACC ALPAKA_FN_INLINE void ALPAKA_ACCELERATOR_NAMESPACE::brokenline::translateKarimaki	(	const TAcc &	acc,
		karimaki_circle_fit &	circle,
		double	x0,
		double	y0,
		riemannFit::Matrix3d &	jacobian
	)

Changes the Karimäki parameters (and consequently their covariance matrix) under a translation of the coordinate system, such that the old origin has coordinates (x0,y0) in the new coordinate system. The formulas are taken from Karimäki V., 1990, Effective circle fitting for particle trajectories, Nucl. Instr. and Meth. A305 (1991) 187.

Parameters

circle	circle fit in the old coordinate system. circle.par(0) is phi, circle.par(1) is d and circle.par(2) is rho.
x0	x coordinate of the translation vector.
y0	y coordinate of the translation vector.
jacobian	passed by reference in order to save stack.

Definition at line 106 of file BrokenLine.h.

References funct::cos(), l1tPhase1JetProducer_cfi::cosPhi, riemannFit::CircleFit::cov, amptDefaultParameters_cff::mu, riemannFit::CircleFit::par, funct::sin(), l1tPhase1JetProducer_cfi::sinPhi, ALPAKA_ACCELERATOR_NAMESPACE::riemannFit::sqr(), mathSSE::sqrt(), and protons_cff::xi.

Referenced by circleFit().

                                                                                                         {
     // Avoid multiple access to the circle.par vector.
     using scalar = typename std::remove_reference<decltype(circle.par(0))>::type;
     scalar phi = circle.par(0);
     scalar dee = circle.par(1);
     scalar rho = circle.par(2);
 
     // Avoid repeated trig. computations
     scalar sinPhi = alpaka::math::sin(acc, phi);
     scalar cosPhi = alpaka::math::cos(acc, phi);
 
     // Intermediate computations for the circle parameters
     scalar deltaPara = x0 * cosPhi + y0 * sinPhi;
     scalar deltaOrth = x0 * sinPhi - y0 * cosPhi + dee;
     scalar tempSmallU = 1 + rho * dee;
     scalar tempC = -rho * y0 + tempSmallU * cosPhi;
     scalar tempB = rho * x0 + tempSmallU * sinPhi;
     scalar tempA = 2. * deltaOrth + rho * (riemannFit::sqr(deltaOrth) + riemannFit::sqr(deltaPara));
     scalar tempU = alpaka::math::sqrt(acc, 1. + rho * tempA);
 
     // Intermediate computations for the error matrix transform
     scalar xi = 1. / (riemannFit::sqr(tempB) + riemannFit::sqr(tempC));
     scalar tempV = 1. + rho * deltaOrth;
     scalar lambda = (0.5 * tempA) / (riemannFit::sqr(1. + tempU) * tempU);
     scalar mu = 1. / (tempU * (1. + tempU)) + rho * lambda;
     scalar zeta = riemannFit::sqr(deltaOrth) + riemannFit::sqr(deltaPara);
     jacobian << xi * tempSmallU * tempV, -xi * riemannFit::sqr(rho) * deltaOrth, xi * deltaPara,
         2. * mu * tempSmallU * deltaPara, 2. * mu * tempV, mu * zeta - lambda * tempA, 0, 0, 1.;
 
     // translated circle parameters
     // phi
     circle.par(0) = alpaka::math::atan2(acc, tempB, tempC);
     // d
     circle.par(1) = tempA / (1 + tempU);
     // rho after translation. It is invariant, so noop
     // circle.par(2)= rho;
 
     // translated error matrix
     circle.cov = jacobian * circle.cov * jacobian.transpose();
   }