From 3d13a1784f5c1aef6d49ce981ed855f258736633 Mon Sep 17 00:00:00 2001 From: Andrea Bocci Date: Thu, 8 Nov 2018 09:09:11 +0100 Subject: [PATCH] Riemann fit rework (#190) The Riemann Fit has been reworked so that both barrel and forward cases are naturally supported without branching. The underlying assumption is the uniform material distribution within the Pixel Tracker. The line fit has been reworked and is now using an ordinary least square fit in the S-Z plane. See the motivations and explanations inside the comments in the code. Additional changes: - code clean up - remove unused functions - fix standalone test of RiemannFit on GPU --- .../PixelTrackFitting/interface/RiemannFit.h | 510 +++++++++++------- .../test/PixelTrackRiemannFit.cc | 388 ++++++++----- .../PixelTrackFitting/test/testEigenGPU.cu | 137 +---- 3 files changed, 581 insertions(+), 454 deletions(-) diff --git a/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h b/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h index 6de1c77bbac12..33f8334c8b5a5 100644 --- a/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h +++ b/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h @@ -47,7 +47,7 @@ struct circle_fit |cov(X0,X0)|cov(Y0,X0)|cov( R,X0)| \n |cov(X0,Y0)|cov(Y0,Y0)|cov( R,Y0)| \n |cov(X0, R)|cov(Y0, R)|cov( R, R)| - */ + */ int64_t q; //!< particle charge double chi2 = 0.0; }; @@ -78,7 +78,6 @@ struct helix_fit double chi2_line = 0.0; Vector4d fast_fit; int64_t q; //!< particle charge - // VectorXd time; // TO FIX just for profiling } __attribute__((aligned(16))); template @@ -119,31 +118,40 @@ __host__ __device__ inline double cross2D(const Vector2d& a, const Vector2d& b) return a.x() * b.y() - a.y() * b.x(); } +/*! Compute the Radiation length in the uniform hypothesis + * + * The Pixel detector, barrel and forward, is considered as an omogeneous + * cilinder 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 lenghts are the path travelled by the particle along + * its trajectory, namely the so called S of the helix in 3D space. + * + * \param 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 lenght has + * no reference point (i.e. it has the dca). + * + * \return incremental radiation lengths that correspond to each segment. + */ -__host__ __device__ inline void computeRadLenEff(const Vector4d& fast_fit, - const double B, - double & radlen_eff, - double & theta, - bool & in_forward) { - double X_barrel = 0.015; - double X_forward = 0.05; - theta = atan(fast_fit(3)); - // atan returns values in [-pi/2, pi/2], we need [0, pi] - theta = theta < 0. ? theta + M_PI : theta; - radlen_eff = X_barrel / std::abs(sin(theta)); - in_forward = (theta <= 0.398 or theta >= 2.743); - if (in_forward) - radlen_eff = X_forward / std::abs(cos(theta)); - assert(radlen_eff > 0.); - double p_t = fast_fit(2) * B; - // We have also to correct the radiation lenght in the x-y plane. Since we - // do not know the angle of incidence of the track at this point, we - // arbitrarily set the correction proportional to the inverse of the - // transerse momentum. The cut-off is at 1 Gev, set using a single Muon Pt - // gun and verifying that, at that momentum, not additional correction is, - // in fact, needed. This is an approximation. - if (std::abs(p_t/1.) < 1.) - radlen_eff /= std::abs(p_t/1.); +__host__ __device__ inline +void computeRadLenUniformMaterial(const VectorNd &length_values, + VectorNd & rad_lengths) { + // Radiation length of the pixel detector in the uniform assumption, with + // 0.06 rad_len at 16 cm + const double XX_0 = 16.f/(0.06); +// const double XX_0 = 1000.*16.f/(0.06); + u_int n = length_values.rows(); + rad_lengths(0) = length_values(0)/XX_0; + for (u_int j = 1; j < n; ++j) { + rad_lengths(j) = std::abs(length_values(j)-length_values(j-1))/XX_0; + } } /*! @@ -151,11 +159,29 @@ __host__ __device__ inline void computeRadLenEff(const Vector4d& fast_fit, 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. */ + __host__ __device__ inline MatrixNd Scatter_cov_line(Matrix2Nd& cov_sz, const Vector4d& fast_fit, VectorNd const& s_arcs, VectorNd const& z_values, + const double theta, const double B) { #if RFIT_DEBUG @@ -164,47 +190,24 @@ __host__ __device__ inline MatrixNd Scatter_cov_line(Matrix2Nd& cov_sz, u_int n = s_arcs.rows(); double p_t = fast_fit(2) * B; double p_2 = p_t * p_t * (1. + 1. / (fast_fit(3) * fast_fit(3))); - double radlen_eff = 0.; - double theta = 0.; - bool in_forward = false; - computeRadLenEff(fast_fit, B, radlen_eff, theta, in_forward); - - const double sig2 = .000225 / p_2 * sqr(1 + 0.038 * log(radlen_eff)) * radlen_eff; - for (u_int k = 0; k < n; ++k) - { - for (u_int l = k; l < n; ++l) - { - for (u_int i = 0; i < std::min(k, l); ++i) - { -#if RFIT_DEBUG - printf("Scatter_cov_line - B: %f\n", B); - printf("Scatter_cov_line - radlen_eff: %f, p_t: %f, p2: %f\n", radlen_eff, p_t, p_2); - printf("Scatter_cov_line - sig2:%f, theta: %f\n", sig2, theta); - printf("Scatter_cov_line - Adding to element %d, %d value %f\n", n + k, n + l, (s_arcs(k) - s_arcs(i)) * (s_arcs(l) - s_arcs(i)) * sig2 / sqr(sqr(sin(theta)))); -#endif - if (in_forward) { - cov_sz(k, l) += (z_values(k) - z_values(i)) * (z_values(l) - z_values(i)) * sig2 / sqr(sqr(cos(theta))); - cov_sz(l, k) = cov_sz(k, l); - } else { - cov_sz(n + k, n + l) += (s_arcs(k) - s_arcs(i)) * (s_arcs(l) - s_arcs(i)) * sig2 / sqr(sqr(sin(theta))); - cov_sz(n + l, n + k) = cov_sz(n + k, n + l); - } - } - } - } + VectorNd rad_lengths_S(n); + // 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 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 sig2_S(n); + sig2_S = .000225 / p_2 * (1.f + 0.038 * rad_lengths_S.array().log()).abs2() * rad_lengths_S.array(); #if RFIT_DEBUG Rfit::printIt(&cov_sz, "Scatter_cov_line - cov_sz: "); #endif Matrix2Nd rot = MatrixXd::Zero(2 * n, 2 * n); for (u_int i = 0; i < n; ++i) { - rot(i, i) = cos(theta); - rot(n + i, n + i) = cos(theta); + rot(i, i) = sin(theta); + rot(n + i, n + i) = sin(theta); u_int j = (i + n); - // Signs seem to be wrong for the off-diagonal element, but we are - // inverting x-y in the input vector, since theta is the angle between - // the z axis and the line, and we are putting the s values, which are Y, - // in the first position. A simple sign flip will take care of it. - rot(i, j) = i < j ? sin(theta) : -sin(theta); + rot(i, j) = i < j ? cos(theta) : -cos(theta); } #if RFIT_DEBUG @@ -212,12 +215,23 @@ __host__ __device__ inline MatrixNd Scatter_cov_line(Matrix2Nd& cov_sz, #endif Matrix2Nd tmp = rot*cov_sz*rot.transpose(); - // We are interested only in the errors in the rotated s -axis which, in - // our formalism, are in the upper square matrix. + for (u_int k = 0; k < n; ++k) + { + for (u_int l = k; l < n; ++l) + { + for (u_int 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. #if RFIT_DEBUG Rfit::printIt(&tmp, "Scatter_cov_line - tmp: "); #endif - return tmp.block(0, 0, n, n); + return tmp.block(n, n, n, n); } /*! @@ -233,13 +247,11 @@ __host__ __device__ inline MatrixNd Scatter_cov_line(Matrix2Nd& cov_sz, \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). - \bug currently works only for points in the barrel. \details Only the tangential component is computed (the radial one is negligible). */ -// X in input TO FIX __host__ __device__ inline MatrixNd Scatter_cov_rad(const Matrix2xNd& p2D, const Vector4d& fast_fit, VectorNd const& rad, @@ -248,24 +260,32 @@ __host__ __device__ inline MatrixNd Scatter_cov_rad(const Matrix2xNd& p2D, u_int n = p2D.cols(); double p_t = fast_fit(2) * B; double p_2 = p_t * p_t * (1. + 1. / (fast_fit(3) * fast_fit(3))); - double radlen_eff = 0.; - double theta = 0.; - bool in_forward = false; - computeRadLenEff(fast_fit, B, radlen_eff, theta, in_forward); + double theta = atan(fast_fit(3)); + theta = theta < 0. ? theta + M_PI : theta; + VectorNd s_values(n); + VectorNd rad_lengths(n); + const Vector2d o(fast_fit(0), fast_fit(1)); + // associated Jacobian, used in weights and errors computation + for (u_int i = 0; i < n; ++i) + { // x + Vector2d p = p2D.block(0, i, 2, 1) - o; + const double cross = cross2D(-o, p); + const double dot = (-o).dot(p); + const double atan2_ = atan2(cross, dot); + s_values(i) = std::abs(atan2_ * fast_fit(2)); + } + computeRadLenUniformMaterial(s_values*sqrt(1. + 1./(fast_fit(3)*fast_fit(3))), rad_lengths); MatrixNd scatter_cov_rad = MatrixXd::Zero(n, n); - const double sig2 = .000225 / p_2 * sqr(1 + 0.038 * log(radlen_eff)) * radlen_eff; + VectorNd sig2(n); + sig2 = .000225 / p_2 * (1.f + 0.038 * rad_lengths.array().log()).abs2() * rad_lengths.array(); for (u_int k = 0; k < n; ++k) { for (u_int l = k; l < n; ++l) { for (u_int i = 0; i < std::min(k, l); ++i) { - if (in_forward) { - scatter_cov_rad(k, l) += (rad(k) - rad(i)) * (rad(l) - rad(i)) * sig2 / sqr(cos(theta)); - } else { - scatter_cov_rad(k, l) += (rad(k) - rad(i)) * (rad(l) - rad(i)) * sig2 / sqr(sin(theta)); - } + scatter_cov_rad(k, l) += (rad(k) - rad(i)) * (rad(l) - rad(i)) * sig2(i) / (sqr(sin(theta))); scatter_cov_rad(l, k) = scatter_cov_rad(k, l); } } @@ -409,23 +429,6 @@ __host__ __device__ inline VectorNd Weight_circle(const MatrixNd& cov_rad_inv) return cov_rad_inv.colwise().sum().transpose(); } -/*! - \brief Compute the points' weights' vector for the line fit (ODR). - Results from a pre-fit is needed in order to take the orthogonal (to the - line) component of the errors. - - \param x_err2 squared errors in the x axis. - \param y_err2 squared errors in the y axis. - \param tan_theta tangent of theta (angle between y axis and line). - - \return weight points' weights' vector for the line fit (ODR). -*/ - -__host__ __device__ inline VectorNd Weight_line(const ArrayNd& x_err2, const ArrayNd& y_err2, const double& tan_theta) -{ - return (1. + sqr(tan_theta)) * 1. / (x_err2 + y_err2 * sqr(tan_theta)); -} - /*! \brief Find particle q considering the sign of cross product between particles velocity (estimated by the first 2 hits) and the vector radius @@ -470,40 +473,6 @@ __host__ __device__ inline void par_uvrtopak(circle_fit& circle, const double B, circle.par = par_pak; } -/*! - \brief Compute the error propagation to obtain the square errors in the - x axis for the line fit. If errors have not been computed in the circle fit - than an'approximation is made. - Further information in attached documentation. - - \param V hits' covariance matrix. - \param circle result of the previous circle fit (only the covariance matrix - is needed) TO FIX - \param J Jacobian of the transformation producing x values. - \param error flag for error computation. - - \return x_err2 squared errors in the x axis. -*/ - -__host__ __device__ inline VectorNd X_err2(const Matrix3Nd& V, const circle_fit& circle, const MatrixNx5d& J, - const bool error, u_int n) -{ - VectorNd x_err2(n); - for (u_int i = 0; i < n; ++i) - { - Matrix5d Cov = MatrixXd::Zero(5, 5); - if (error) - Cov.block(0, 0, 3, 3) = circle.cov; - Cov(3, 3) = V(i, i); - Cov(4, 4) = V(i + n, i + n); - Cov(3, 4) = Cov(4, 3) = V(i, i + n); - Eigen::Matrix tmp; - tmp = J.row(i) * Cov * J.row(i).transpose().eval(); - x_err2(i) = tmp(0, 0); - } - return x_err2; -} - /*! \brief Compute the eigenvector associated to the minimum eigenvalue. @@ -1004,7 +973,7 @@ __host__ __device__ inline circle_fit Circle_fit(const Matrix2xNd& hits2D, { const double t = 1. / h; J3 << -v2x2_inv, 0, v(0) * sqr(v2x2_inv) * 2., 0, 0, -v2x2_inv, v(1) * sqr(v2x2_inv) * 2., 0, - 0, 0, -h * sqr(v2x2_inv) * 2. - (2. * c + v(2)) * v2x2_inv * t, -t; + v(0)*v2x2_inv*t, v(1)*v2x2_inv*t, -h * sqr(v2x2_inv) * 2. - (2. * c + v(2)) * v2x2_inv * t, -t; } printIt(&J3, "circle_fit - J3:"); @@ -1059,21 +1028,22 @@ __host__ __device__ inline circle_fit Circle_fit(const Matrix2xNd& hits2D, errors. */ -__host__ __device__ inline line_fit Line_fit(const Matrix3xNd& hits, - const Matrix3Nd& hits_cov, - const circle_fit& circle, - const Vector4d& fast_fit, - const double B, - const bool error = true) +__host__ __device__ inline line_fit Line_fit_odr(const Matrix3xNd& hits, + const Matrix3Nd& hits_cov, + const circle_fit& circle, + const Vector4d& fast_fit, + const double B, + const bool error = true) { u_int n = hits.cols(); + double theta = -circle.q*atan(fast_fit(3)); + theta = theta < 0. ? theta + M_PI : theta; // PROJECTION ON THE CILINDER - Matrix2xNd p2D(2, n); - MatrixNx5d Jx(n, 5); + Matrix2xNd p2D = MatrixXd::Zero(2, n); + Eigen::Matrix Jx; #if RFIT_DEBUG printf("Line_fit - B: %g\n", B); - printIt(&hits, "Line_fit points: "); printIt(&hits_cov, "Line_fit covs: "); #endif @@ -1085,8 +1055,11 @@ __host__ __device__ inline line_fit Line_fit(const Matrix3xNd& hits, const Vector2d o(circle.par(0), circle.par(1)); // associated Jacobian, used in weights and errors computation + Matrix2Nd cov_sz = MatrixXd::Zero(2 * n, 2 * n); for (u_int i = 0; i < n; ++i) { // x + Matrix6d Cov = MatrixXd::Zero(6, 6); + Matrix2d Cov_sz_single = MatrixXd::Zero(2, 2); Vector2d p = hits.block(0, i, 2, 1) - o; const double cross = cross2D(-o, p); const double dot = (-o).dot(p); @@ -1095,9 +1068,9 @@ __host__ __device__ inline line_fit Line_fit(const Matrix3xNd& hits, const double atan2_ = -circle.q * atan2(cross, dot); p2D(0, i) = atan2_ * circle.par(2); - // associated Jacobian, used in weights and errors computation + // associated Jacobian, used in weights and errors- computation const double temp0 = -circle.q * 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 + double d_X0 = 0., d_Y0 = 0., d_R = 0.; // good approximation for big pt and eta if (error) { d_X0 = -temp0 * ((p(1) + o(1)) * dot - (p(0) - o(0)) * cross); @@ -1106,7 +1079,19 @@ __host__ __device__ inline line_fit Line_fit(const Matrix3xNd& hits, } const double d_x = temp0 * (o(1) * dot + o(0) * cross); const double d_y = temp0 * (-o(0) * dot + o(1) * cross); - Jx.row(i) << d_X0, d_Y0, d_R, d_x, d_y; + Jx << d_X0, d_Y0, d_R, d_x, d_y, 0., 0., 0., 0., 0., 0., 1.; +// Jx << d_X0, d_Y0, d_R, p(1)/p.norm(), -p(0)/p.norm(), 0, 0, 0, 0, 0, 0, 1.; + Cov.block(0, 0, 3, 3) = circle.cov; + Cov(3, 3) = hits_cov(i, i); + Cov(4, 4) = hits_cov(i + n, i + n); + Cov(5, 5) = hits_cov(i + 2*n, i + 2*n); + Cov(3, 4) = Cov(4, 3) = hits_cov(i, i + n); + Cov(3, 5) = Cov(5, 3) = hits_cov(i, i + 2*n); + Cov(4, 5) = Cov(5, 4) = hits_cov(i + n, i + 2*n); + Cov_sz_single = Jx * Cov * Jx.transpose(); + cov_sz(i, i) = Cov_sz_single(0, 0); + cov_sz(i + n, i + n) = Cov_sz_single(1, 1); + cov_sz(i, i + n) = cov_sz(i + n, i) = Cov_sz_single(0, 1); } // Math of d_{X0,Y0,R,x,y} all verified by hand @@ -1114,43 +1099,25 @@ __host__ __device__ inline line_fit Line_fit(const Matrix3xNd& hits, p2D.row(1) = hits.row(2); // WEIGHT COMPUTATION - Matrix2Nd cov_sz = MatrixXd::Zero(2 * n, 2 * n); - VectorNd x_err2 = X_err2(hits_cov, circle, Jx, error, n); - VectorNd y_err2 = hits_cov.block(2 * n, 2 * n, n, n).diagonal(); - cov_sz.block(0, 0, n, n) = x_err2.asDiagonal(); - cov_sz.block(n, n, n, n) = y_err2.asDiagonal(); #if RFIT_DEBUG printIt(&cov_sz, "line_fit - cov_sz:"); #endif - MatrixNd cov_with_ms = Scatter_cov_line(cov_sz, fast_fit, p2D.row(0), p2D.row(1), B); + MatrixNd cov_with_ms = Scatter_cov_line(cov_sz, fast_fit, p2D.row(0), p2D.row(1), theta, B); #if RFIT_DEBUG printIt(&cov_with_ms, "line_fit - cov_with_ms: "); #endif - Matrix4d G, G4; - G4 = cov_with_ms.inverse(); + Matrix4d G; + G = cov_with_ms.inverse(); #if RFIT_DEBUG - printIt(&G4, "line_fit - cov_with_ms.inverse():"); + printIt(&G, "line_fit - cov_with_ms.inverse():"); #endif - double renorm = G4.sum(); - G4 *= 1. / renorm; + double renorm = G.sum(); + G *= 1. / renorm; #if RFIT_DEBUG - printIt(&G4, "line_fit - G4:"); + printIt(&G, "line_fit - G4:"); #endif - G = G4; - const VectorNd weight = Weight_circle(G); - - VectorNd err2_inv = cov_with_ms.diagonal(); - err2_inv = err2_inv.cwiseInverse(); -// const VectorNd err2_inv = Weight_line(x_err2, y_err2, fast_fit(3)); -// const VectorNd weight = err2_inv * 1. / err2_inv.sum(); - -#if RFIT_DEBUG - printIt(&x_err2, "Line_fit - x_err2: "); - printIt(&y_err2, "Line_fit - y_err2: "); - printIt(&err2_inv, "Line_fit - err2_inv: "); - printIt(&weight, "Line_fit - weight: "); -#endif + const VectorNd weight = Weight_circle(G); // COST FUNCTION @@ -1162,16 +1129,12 @@ __host__ __device__ inline line_fit Line_fit(const Matrix3xNd& hits, const Matrix2xNd X = p2D.colwise() - r0; Matrix2d A = Matrix2d::Zero(); A = X * G * X.transpose(); -// for (u_int i = 0; i < n; ++i) -// { -// A += err2_inv(i) * (X.col(i) * X.col(i).transpose()); -// } #if RFIT_DEBUG printIt(&A, "Line_fit - A: "); #endif - // minimize + // minimize. v is normalized!! double chi2; Vector2d v = min_eigen2D(A, chi2); #if RFIT_DEBUG @@ -1179,7 +1142,6 @@ __host__ __device__ inline line_fit Line_fit(const Matrix3xNd& hits, printf("Line_fit chi2: %e\n", chi2); #endif - // n *= (chi2>0) ? 1 : -1; //TO FIX // This hack to be able to run on GPU where the automatic assignment to a // double from the vector multiplication is not working. Matrix cm; @@ -1189,8 +1151,8 @@ __host__ __device__ inline line_fit Line_fit(const Matrix3xNd& hits, // COMPUTE LINE PARAMETER line_fit line; line.par << -v(0) / v(1), // cotan(theta)) - -c * sqrt(sqr(v(0)) + sqr(v(1))) * 1. / v(1); // Zip - line.chi2 = abs(chi2); + -c / v(1); // Zip + line.chi2 = abs(chi2*renorm); #if RFIT_DEBUG printIt(&(line.par), "Line_fit - line.par: "); printf("Line_fit - v norm: %e\n", sqrt(v(0)*v(0) + v(1)*v(1))); @@ -1206,19 +1168,21 @@ __host__ __device__ inline line_fit Line_fit(const Matrix3xNd& hits, { // The norm is taken from Chernov, properly adapted to the weights case. double norm = v.transpose() * A * v; - norm /= weight.sum(); +// double norm_empirical = cov_with_ms.diagonal().mean(); #if RFIT_DEBUG + printf("Chi_2: %g\n", chi2); + printf("Norm: %g\n", norm); + printf("weight.sum(): %g\n", weight.sum()); printf("Line_fit - norm: %e\n", norm); #endif - const double sig2 = 1. / (A(0, 0) + A(1, 1)) * norm; + + const double sig2 = norm/(A(0,0) + A(1,1)); C(0, 0) = sig2 * v1_2; C(1, 1) = sig2 * v0_2; - C(0, 1) = C(1, 0) = -sig2 * v(0) * v(1); - const VectorNd weight_2 = (weight).array().square(); - const Vector2d C0(weight_2.dot(x_err2), weight_2.dot(y_err2)); - C.block(0, 2, 2, 1) = C.block(2, 0, 1, 2).transpose() = -C.block(0, 0, 2, 2) * r0; - Matrix tmp = (r0.transpose() * C.block(0, 0, 2, 2) * r0); - C(2, 2) = v0_2 * C0(0) + v1_2 * C0(1) + C0(0) * C(0, 0) + C0(1) * C(1, 1) + tmp(0, 0); + C(1, 0) = C(0, 1) = -sig2 * v(0) * v(1); + C(2, 2) = sig2 * (v(0)*r0(1)-v(1)*r0(0))*(v(0)*r0(1)-v(1)*r0(0)) + (sig2/n)*(A(0,0)+A(1,1)); + C(0, 2) = C(2, 0) = sig2*(v(0)*r0(1)-v(1)*r0(0))*v(1); + C(1, 2) = C(2, 1) = - sig2*(v(0)*r0(1)-v(1)*r0(0))*v(0); } #if RFIT_DEBUG printIt(&C, "line_fit - C:"); @@ -1228,9 +1192,7 @@ __host__ __device__ inline line_fit Line_fit(const Matrix3xNd& hits, { const double t0 = 1. / v(1); const double t1 = sqr(t0); - const double sqrt_ = sqrt(v1_2 + v0_2); - const double t2 = 1. / sqrt_; - J << -t0, v(0) * t1, 0, -c * v(0) * t0 * t2, v0_2 * c * t1 * t2, -sqrt_ * t0; + J << -t0, v(0) * t1, 0., 0., c * t1, -t0; } Matrix JT = J.transpose().eval(); #if RFIT_DEBUG @@ -1245,6 +1207,184 @@ __host__ __device__ inline line_fit Line_fit(const Matrix3xNd& hits, return line; } +/*! \brief 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. + */ + +__host__ __device__ inline line_fit Line_fit(const Matrix3xNd& hits, + const Matrix3Nd& hits_cov, + const circle_fit& circle, + const Vector4d& fast_fit, + const double B, + const bool error = true) { + auto n = hits.cols(); + double theta = -circle.q*atan(fast_fit(3)); + theta = theta < 0. ? theta + M_PI : 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 p2D(2, n); + Eigen::Matrix Jx; + + p2D << MatrixXd::Zero(2, n); + Jx << MatrixXd::Zero(2, 6); + +#if RFIT_DEBUG + printf("Line_fit - B: %g\n", B); + printIt(&hits, "Line_fit points: "); + printIt(&hits_cov, "Line_fit covs: "); +#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 o(circle.par(0), circle.par(1)); + + // associated Jacobian, used in weights and errors computation + Matrix2Nd cov_sz = MatrixXd::Zero(2 * n, 2 * n); + Matrix6d Cov(6,6); + Matrix2d Cov_sz_single(2, 2); + for (u_int i = 0; i < n; ++i) + { + Vector2d p = hits.block(0, i, 2, 1) - o; + const double cross = cross2D(-o, p); + const double dot = (-o).dot(p); + // 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 atan2_ = -circle.q * atan2(cross, dot); +// p2D.coeffRef(1, i) = atan2_ * circle.par(2); + p2D(0, i) = atan2_ * circle.par(2); + + // associated Jacobian, used in weights and errors- computation + const double temp0 = -circle.q * 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 * ((p(1) + o(1)) * dot - (p(0) - o(0)) * cross); + d_Y0 = temp0 * ((p(0) + o(0)) * dot - (o(1) - p(1)) * cross); + d_R = atan2_; + } + const double d_x = temp0 * (o(1) * dot + o(0) * cross); + const double d_y = temp0 * (-o(0) * dot + o(1) * cross); + Jx << d_X0, d_Y0, d_R, d_x, d_y, 0., 0., 0., 0., 0., 0., 1.; + + Cov << MatrixXd::Zero(6, 6); + Cov_sz_single << MatrixXd::Zero(2, 2); + Cov.block(0, 0, 3, 3) = circle.cov; + Cov(3, 3) = hits_cov(i, i); // x errors + Cov(4, 4) = hits_cov(i + n, i + n); // y errors + Cov(5, 5) = hits_cov(i + 2*n, i + 2*n); // z errors + Cov(3, 4) = Cov(4, 3) = hits_cov(i, i + n); // cov_xy + Cov(3, 5) = Cov(5, 3) = hits_cov(i, i + 2*n); // cov_xz + Cov(4, 5) = Cov(5, 4) = hits_cov(i + n, i + 2*n); // cov_yz + Cov_sz_single = Jx * Cov * Jx.transpose(); + cov_sz(i, i) = Cov_sz_single(0, 0); + cov_sz(i + n, i + n) = Cov_sz_single(1, 1); + cov_sz(i, i + n) = cov_sz(i + n, i) = Cov_sz_single(0, 1); + } + // 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 cov_with_ms = Scatter_cov_line(cov_sz, fast_fit, p2D.row(0), p2D.row(1), theta, B); +#if RFIT_DEBUG + printIt(&cov_sz, "line_fit - cov_sz:"); + printIt(&cov_with_ms, "line_fit - cov_with_ms: "); +#endif + + // Prepare the Rotation Matrix to rotate the points + Eigen::Matrix rot = Eigen::Matrix::Zero(); + rot << sin(theta), cos(theta), -cos(theta), sin(theta); + + // Rotate Points with the shape [2, n] + Matrix2xNd p2D_rot = rot*p2D; + +#if 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 A(2,n); + A << MatrixXd::Ones(1, n), p2D_rot.row(0); // rotated s values + +#if RFIT_DEBUG + printIt(&A, "A Matrix:"); +#endif + + // Build A^T V-1 A, where V-1 is the covariance of only the Y components. + MatrixNd Vy_inv = cov_with_ms.inverse(); + Eigen::Matrix Inv_Cov = A*Vy_inv*A.transpose(); + + // Compute the Covariance Matrix of the fit parameters + Eigen::Matrix Cov_params = Inv_Cov.inverse(); + + // Now Compute the Parameters in the form [2,1] + // The first component is q. + // The second component is m. + Eigen::Matrix sol = Cov_params*A*Vy_inv*p2D_rot.row(1).transpose(); + + +#if RFIT_DEBUG + printIt(&sol, "Rotated solutions:"); +#endif + + // We need now to transfer back the results in the original s-z plane + auto common_factor = 1./(sin(theta)-sol(1,0)*cos(theta)); + Matrix J = Matrix::Zero(); + J << 0., common_factor*common_factor, common_factor, sol(0,0)*cos(theta)*common_factor*common_factor; + + double m = common_factor*(sol(1,0)*sin(theta)+cos(theta)); + double q = common_factor*sol(0,0); + auto cov_mq = J * Cov_params * J.transpose(); + + VectorNd res = p2D_rot.row(1).transpose() - A.transpose() * sol; + double chi2 = res.transpose()*Vy_inv*res; + chi2 = chi2 / float(n); + + line_fit line; + line.par << m, q; + line.cov << cov_mq; + line.chi2 = chi2; + +#if RFIT_DEBUG + printf("Common_factor: %g\n", common_factor); + printIt(&J, "Jacobian:"); + printIt(&sol, "Rotated solutions:"); + printIt(&Cov_params, "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; +} + /*! \brief Helix fit by three step: -fast pre-fit (see Fast_fit() for further info); \n diff --git a/RecoPixelVertexing/PixelTrackFitting/test/PixelTrackRiemannFit.cc b/RecoPixelVertexing/PixelTrackFitting/test/PixelTrackRiemannFit.cc index b27ed52473388..f71e3f082ada4 100644 --- a/RecoPixelVertexing/PixelTrackFitting/test/PixelTrackRiemannFit.cc +++ b/RecoPixelVertexing/PixelTrackFitting/test/PixelTrackRiemannFit.cc @@ -6,7 +6,11 @@ #include #include // unique_ptr +#include +#include + #include "RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h" +//#include "RecoPixelVertexing/PixelTrackFitting/interface/BrokenLine.h" using namespace std; using namespace Eigen; @@ -164,161 +168,261 @@ Matrix New_par(const Matrix& gen_par, const int& cha return new_par; } +template +void computePull(std::array & fit, const char * label, + int n_, int iteration, const Vector5d & true_par) { + Matrix score(41, iteration); + + std::string histo_name("Phi Pull"); + histo_name += label; + TH1F phi_pull(histo_name.data(), histo_name.data(), 100, -10., 10.); + histo_name = "dxy Pull "; + histo_name += label; + TH1F dxy_pull(histo_name.data(), histo_name.data(), 100, -10., 10.); + histo_name = "dz Pull "; + histo_name += label; + TH1F dz_pull(histo_name.data(), histo_name.data(), 100, -10., 10.); + histo_name = "Theta Pull "; + histo_name += label; + TH1F theta_pull(histo_name.data(), histo_name.data(), 100, -10., 10.); + histo_name = "Pt Pull "; + histo_name += label; + TH1F pt_pull(histo_name.data(), histo_name.data(), 100, -10., 10.); + histo_name = "Phi Error "; + histo_name += label; + TH1F phi_error(histo_name.data(), histo_name.data(), 100, 0., 0.1); + histo_name = "dxy error "; + histo_name += label; + TH1F dxy_error(histo_name.data(), histo_name.data(), 100, 0., 0.1); + histo_name = "dz error "; + histo_name += label; + TH1F dz_error(histo_name.data(), histo_name.data(), 100, 0., 0.1); + histo_name = "Theta error "; + histo_name += label; + TH1F theta_error(histo_name.data(), histo_name.data(), 100, 0., 0.1); + histo_name = "Pt error "; + histo_name += label; + TH1F pt_error(histo_name.data(), histo_name.data(), 100, 0., 0.1); + for (int x = 0; x < iteration; x++) { + // Compute PULLS information + score(0, x) = (fit[x].par(0) - true_par(0)) / sqrt(fit[x].cov(0, 0)); + score(1, x) = (fit[x].par(1) - true_par(1)) / sqrt(fit[x].cov(1, 1)); + score(2, x) = (fit[x].par(2) - true_par(2)) / sqrt(fit[x].cov(2, 2)); + score(3, x) = (fit[x].par(3) - true_par(3)) / sqrt(fit[x].cov(3, 3)); + score(4, x) = (fit[x].par(4) - true_par(4)) / sqrt(fit[x].cov(4, 4)); + phi_pull.Fill(score(0, x)); + dxy_pull.Fill(score(1, x)); + pt_pull.Fill(score(2, x)); + theta_pull.Fill(score(3, x)); + dz_pull.Fill(score(4, x)); + phi_error.Fill(sqrt(fit[x].cov(0, 0))); + dxy_error.Fill(sqrt(fit[x].cov(1, 1))); + pt_error.Fill(sqrt(fit[x].cov(2, 2))); + theta_error.Fill(sqrt(fit[x].cov(3, 3))); + dz_error.Fill(sqrt(fit[x].cov(4, 4))); + score(5, x) = + (fit[x].par(0) - true_par(0)) * (fit[x].par(1) - true_par(1)) / (fit[x].cov(0, 1)); + score(6, x) = + (fit[x].par(0) - true_par(0)) * (fit[x].par(2) - true_par(2)) / (fit[x].cov(0, 2)); + score(7, x) = + (fit[x].par(1) - true_par(1)) * (fit[x].par(2) - true_par(2)) / (fit[x].cov(1, 2)); + score(8, x) = + (fit[x].par(3) - true_par(3)) * (fit[x].par(4) - true_par(4)) / (fit[x].cov(3, 4)); + score(9, x) = fit[x].chi2_circle; + score(25, x) = fit[x].chi2_line; + score(10, x) = sqrt(fit[x].cov(0, 0)) / fit[x].par(0) * 100; + score(13, x) = sqrt(fit[x].cov(3, 3)) / fit[x].par(3) * 100; + score(14, x) = sqrt(fit[x].cov(4, 4)) / fit[x].par(4) * 100; + score(15, x) = (fit[x].par(0) - true_par(0)) * (fit[x].par(3) - true_par(3)) / + sqrt(fit[x].cov(0, 0)) / sqrt(fit[x].cov(3, 3)); + score(16, x) = (fit[x].par(1) - true_par(1)) * (fit[x].par(3) - true_par(3)) / + sqrt(fit[x].cov(1, 1)) / sqrt(fit[x].cov(3, 3)); + score(17, x) = (fit[x].par(2) - true_par(2)) * (fit[x].par(3) - true_par(3)) / + sqrt(fit[x].cov(2, 2)) / sqrt(fit[x].cov(3, 3)); + score(18, x) = (fit[x].par(0) - true_par(0)) * (fit[x].par(4) - true_par(4)) / + sqrt(fit[x].cov(0, 0)) / sqrt(fit[x].cov(4, 4)); + score(19, x) = (fit[x].par(1) - true_par(1)) * (fit[x].par(4) - true_par(4)) / + sqrt(fit[x].cov(1, 1)) / sqrt(fit[x].cov(4, 4)); + score(20, x) = (fit[x].par(2) - true_par(2)) * (fit[x].par(4) - true_par(4)) / + sqrt(fit[x].cov(2, 2)) / sqrt(fit[x].cov(4, 4)); + score(21, x) = (fit[x].par(0) - true_par(0)) * (fit[x].par(1) - true_par(1)) / + sqrt(fit[x].cov(0, 0)) / sqrt(fit[x].cov(1, 1)); + score(22, x) = (fit[x].par(0) - true_par(0)) * (fit[x].par(2) - true_par(2)) / + sqrt(fit[x].cov(0, 0)) / sqrt(fit[x].cov(2, 2)); + score(23, x) = (fit[x].par(1) - true_par(1)) * (fit[x].par(2) - true_par(2)) / + sqrt(fit[x].cov(1, 1)) / sqrt(fit[x].cov(2, 2)); + score(24, x) = (fit[x].par(3) - true_par(3)) * (fit[x].par(4) - true_par(4)) / + sqrt(fit[x].cov(3, 3)) / sqrt(fit[x].cov(4, 4)); + score(30, x) = fit[x].par(0); + score(31, x) = fit[x].par(1); + score(32, x) = fit[x].par(2); + score(33, x) = fit[x].par(3); + score(34, x) = fit[x].par(4); + score(35, x) = sqrt(fit[x].cov(0,0)); + score(36, x) = sqrt(fit[x].cov(1,1)); + score(37, x) = sqrt(fit[x].cov(2,2)); + score(38, x) = sqrt(fit[x].cov(3,3)); + score(39, x) = sqrt(fit[x].cov(4,4)); + + } + + double phi_ = score.row(0).mean(); + double a_ = score.row(1).mean(); + double pt_ = score.row(2).mean(); + double coT_ = score.row(3).mean(); + double Zip_ = score.row(4).mean(); + std::cout << std::setprecision(5) << std::scientific << label << " AVERAGE FITTED VALUES: \n" + << "phi: " << score.row(30).mean() << " +/- " << score.row(35).mean() << " [+/-] " << sqrt(score.row(35).array().abs2().mean() - score.row(35).mean()*score.row(35).mean()) << std::endl + << "d0: " << score.row(31).mean() << " +/- " << score.row(36).mean() << " [+/-] " << sqrt(score.row(36).array().abs2().mean() - score.row(36).mean()*score.row(36).mean()) << std::endl + << "pt: " << score.row(32).mean() << " +/- " << score.row(37).mean() << " [+/-] " << sqrt(score.row(37).array().abs2().mean() - score.row(37).mean()*score.row(37).mean()) << std::endl + << "coT: " << score.row(33).mean() << " +/- " << score.row(38).mean() << " [+/-] " << sqrt(score.row(38).array().abs2().mean() - score.row(38).mean()*score.row(38).mean()) << std::endl + << "Zip: " << score.row(34).mean() << " +/- " << score.row(39).mean() << " [+/-] " << sqrt(score.row(39).array().abs2().mean() - score.row(39).mean()*score.row(39).mean()) << std::endl; + + Matrix5d correlation; + correlation << 1., score.row(21).mean(), score.row(22).mean(), score.row(15).mean(), + score.row(20).mean(), score.row(21).mean(), 1., score.row(23).mean(), score.row(16).mean(), + score.row(19).mean(), score.row(22).mean(), score.row(23).mean(), 1., score.row(17).mean(), + score.row(20).mean(), score.row(15).mean(), score.row(16).mean(), score.row(17).mean(), 1., + score.row(24).mean(), score.row(18).mean(), score.row(19).mean(), score.row(20).mean(), + score.row(24).mean(), 1.; + + cout << "\n" << label << " PULLS (mean, sigma, relative_error):\n" + << "phi: " << phi_ << " " + << sqrt((score.row(0).array() - phi_).square().sum() / (iteration - 1)) << " " + << abs(score.row(10).mean()) << "%\n" + << "a0 : " << a_ << " " + << sqrt((score.row(1).array() - a_).square().sum() / (iteration - 1)) << " " + << abs(score.row(11).mean()) << "%\n" + << "pt : " << pt_ << " " + << sqrt((score.row(2).array() - pt_).square().sum() / (iteration - 1)) << " " + << abs(score.row(12).mean()) << "%\n" + << "coT: " << coT_ << " " + << sqrt((score.row(3).array() - coT_).square().sum() / (iteration - 1)) << " " + << abs(score.row(13).mean()) << "%\n" + << "Zip: " << Zip_ << " " + << sqrt((score.row(4).array() - Zip_).square().sum() / (iteration - 1)) << " " + << abs(score.row(14).mean()) << "%\n\n" + << "cov(phi,a0)_: " << score.row(5).mean() << "\n" + << "cov(phi,pt)_: " << score.row(6).mean() << "\n" + << "cov(a0,pt)_: " << score.row(7).mean() << "\n" + << "cov(coT,Zip)_: " << score.row(8).mean() << "\n\n" + << "chi2_circle: " << score.row(9).mean() << " vs " << n_ - 3 << "\n" + << "chi2_line: " << score.row(25).mean() << " vs " << n_ - 2 << "\n\n" + << "correlation matrix:\n" + << correlation << "\n\n" + << endl; + + phi_pull.Fit("gaus", "Q"); + dxy_pull.Fit("gaus", "Q"); + dz_pull.Fit("gaus", "Q"); + theta_pull.Fit("gaus", "Q"); + pt_pull.Fit("gaus", "Q"); + phi_pull.Write(); + dxy_pull.Write(); + dz_pull.Write(); + theta_pull.Write(); + pt_pull.Write(); + phi_error.Write(); + dxy_error.Write(); + dz_error.Write(); + theta_error.Write(); + pt_error.Write(); +} + + void test_helix_fit() { int n_; - int iteration; - int debug2 = 0; bool return_err; const double B_field = 3.8 * c_speed / pow(10, 9) / 100; Matrix gen_par; Vector5d true_par; Vector5d err; -// while (1) { - generator.seed(1); - int debug = 0; - debug2 = 0; - std::cout << std::setprecision(6); - cout << "_________________________________________________________________________\n"; - cout << "n x(cm) y(cm) z(cm) phi(grad) R(Gev/c) eta iteration return_err debug" << endl; -// cin >> n_ >> gen_par(0) >> gen_par(1) >> gen_par(2) >> gen_par(3) >> gen_par(4) >> gen_par(5) >> -// iteration >> return_err >> debug2; - n_ = 4; - gen_par(0) = -0.1; // x - gen_par(1) = 0.1; // y - gen_par(2) = -1.; // z - gen_par(3) = 45.; // phi - gen_par(4) = 10.; // R (p_t) - gen_par(5) = 1.; // eta - iteration = 1; - return_err = 1; - debug2 = 1; - - iteration *= 10; - gen_par = New_par(gen_par, 1, B_field); - true_par = True_par(gen_par, 1, B_field); - Matrix3xNd hits; - Matrix3Nd hits_cov; - unique_ptr helix(new helix_fit[iteration]); -// helix_fit* helix = new helix_fit[iteration]; - Matrix score(41, iteration); - - for (int i = 0; i < iteration; i++) { - if (debug2 == 1 && i == (iteration - 1)) { - debug = 1; - } - hits_gen gen; - gen = Hits_gen(n_, gen_par); -// gen.hits = MatrixXd::Zero(3, 4); -// gen.hits_cov = MatrixXd::Zero(3 * 4, 3 * 4); -// gen.hits.col(0) << 1.82917642593, 2.0411875248, 7.18495464325; -// gen.hits.col(1) << 4.47041416168, 4.82704305649, 18.6394691467; -// gen.hits.col(2) << 7.25991010666, 7.74653434753, 30.6931324005; -// gen.hits.col(3) << 8.99161434174, 9.54262828827, 38.1338043213; - helix[i] = Rfit::Helix_fit(gen.hits, gen.hits_cov, B_field, return_err); + generator.seed(1); + std::cout << std::setprecision(6); + cout << "_________________________________________________________________________\n"; + cout << "n x(cm) y(cm) z(cm) phi(grad) R(Gev/c) eta iteration return_err debug" << endl; + cout << "hits: "; + cin >> n_; + cout << "x: "; + cin >> gen_par(0); + cout << "y: "; + cin >> gen_par(1); + cout << "z: "; + cin >> gen_par(2); + cout << "phi: "; + cin >> gen_par(3); + cout << "p_t: "; + cin >> gen_par(4); + cout << "eta: "; + cin >> gen_par(5); + // + /* + n_ = 4; + gen_par(0) = -0.1; // x + gen_par(1) = 0.1; // y + gen_par(2) = -1.; // z + gen_par(3) = 45.; // phi + gen_par(4) = 10.; // R (p_t) + gen_par(5) = 1.; // eta + iteration = 1; + */ + return_err = 1; - if (debug) - cout << std::setprecision(10) - << "phi: " << helix[i].par(0) << " +/- " << sqrt(helix[i].cov(0, 0)) << " vs " - << true_par(0) << endl - << "Tip: " << helix[i].par(1) << " +/- " << sqrt(helix[i].cov(1, 1)) << " vs " - << true_par(1) << endl - << "p_t: " << helix[i].par(2) << " +/- " << sqrt(helix[i].cov(2, 2)) << " vs " - << true_par(2) << endl - << "theta:" << helix[i].par(3) << " +/- " << sqrt(helix[i].cov(3, 3)) << " vs " - << true_par(3) << endl - << "Zip: " << helix[i].par(4) << " +/- " << sqrt(helix[i].cov(4, 4)) << " vs " - << true_par(4) << endl - << "charge:" << helix[i].q << " vs 1" << endl - << "covariance matrix:" << endl - << helix[i].cov << endl - << "Initial hits:\n" << gen.hits << endl - << "Initial Covariance:\n" << gen.hits_cov << endl; - } - - for (int x = 0; x < iteration; x++) { - // Compute PULLS information - score(0, x) = (helix[x].par(0) - true_par(0)) / sqrt(helix[x].cov(0, 0)); - score(1, x) = (helix[x].par(1) - true_par(1)) / sqrt(helix[x].cov(1, 1)); - score(2, x) = (helix[x].par(2) - true_par(2)) / sqrt(helix[x].cov(2, 2)); - score(3, x) = (helix[x].par(3) - true_par(3)) / sqrt(helix[x].cov(3, 3)); - score(4, x) = (helix[x].par(4) - true_par(4)) / sqrt(helix[x].cov(4, 4)); - score(5, x) = - (helix[x].par(0) - true_par(0)) * (helix[x].par(1) - true_par(1)) / (helix[x].cov(0, 1)); - score(6, x) = - (helix[x].par(0) - true_par(0)) * (helix[x].par(2) - true_par(2)) / (helix[x].cov(0, 2)); - score(7, x) = - (helix[x].par(1) - true_par(1)) * (helix[x].par(2) - true_par(2)) / (helix[x].cov(1, 2)); - score(8, x) = - (helix[x].par(3) - true_par(3)) * (helix[x].par(4) - true_par(4)) / (helix[x].cov(3, 4)); - score(9, x) = helix[x].chi2_circle; - score(25, x) = helix[x].chi2_line; - score(10, x) = sqrt(helix[x].cov(0, 0)) / helix[x].par(0) * 100; - score(13, x) = sqrt(helix[x].cov(3, 3)) / helix[x].par(3) * 100; - score(14, x) = sqrt(helix[x].cov(4, 4)) / helix[x].par(4) * 100; - score(15, x) = (helix[x].par(0) - true_par(0)) * (helix[x].par(3) - true_par(3)) / - sqrt(helix[x].cov(0, 0)) / sqrt(helix[x].cov(3, 3)); - score(16, x) = (helix[x].par(1) - true_par(1)) * (helix[x].par(3) - true_par(3)) / - sqrt(helix[x].cov(1, 1)) / sqrt(helix[x].cov(3, 3)); - score(17, x) = (helix[x].par(2) - true_par(2)) * (helix[x].par(3) - true_par(3)) / - sqrt(helix[x].cov(2, 2)) / sqrt(helix[x].cov(3, 3)); - score(18, x) = (helix[x].par(0) - true_par(0)) * (helix[x].par(4) - true_par(4)) / - sqrt(helix[x].cov(0, 0)) / sqrt(helix[x].cov(4, 4)); - score(19, x) = (helix[x].par(1) - true_par(1)) * (helix[x].par(4) - true_par(4)) / - sqrt(helix[x].cov(1, 1)) / sqrt(helix[x].cov(4, 4)); - score(20, x) = (helix[x].par(2) - true_par(2)) * (helix[x].par(4) - true_par(4)) / - sqrt(helix[x].cov(2, 2)) / sqrt(helix[x].cov(4, 4)); - score(21, x) = (helix[x].par(0) - true_par(0)) * (helix[x].par(1) - true_par(1)) / - sqrt(helix[x].cov(0, 0)) / sqrt(helix[x].cov(1, 1)); - score(22, x) = (helix[x].par(0) - true_par(0)) * (helix[x].par(2) - true_par(2)) / - sqrt(helix[x].cov(0, 0)) / sqrt(helix[x].cov(2, 2)); - score(23, x) = (helix[x].par(1) - true_par(1)) * (helix[x].par(2) - true_par(2)) / - sqrt(helix[x].cov(1, 1)) / sqrt(helix[x].cov(2, 2)); - score(24, x) = (helix[x].par(3) - true_par(3)) * (helix[x].par(4) - true_par(4)) / - sqrt(helix[x].cov(3, 3)) / sqrt(helix[x].cov(4, 4)); - } + const int iteration = 5000; + gen_par = New_par(gen_par, 1, B_field); + true_par = True_par(gen_par, 1, B_field); + Matrix3xNd hits; + Matrix3Nd hits_cov; + std::array helixRiemann_fit; +// std::array helixBrokenLine_fit; - double phi_ = score.row(0).mean(); - double a_ = score.row(1).mean(); - double pt_ = score.row(2).mean(); - double coT_ = score.row(3).mean(); - double Zip_ = score.row(4).mean(); - Matrix5d correlation; - correlation << 1., score.row(21).mean(), score.row(22).mean(), score.row(15).mean(), - score.row(20).mean(), score.row(21).mean(), 1., score.row(23).mean(), score.row(16).mean(), - score.row(19).mean(), score.row(22).mean(), score.row(23).mean(), 1., score.row(17).mean(), - score.row(20).mean(), score.row(15).mean(), score.row(16).mean(), score.row(17).mean(), 1., - score.row(24).mean(), score.row(18).mean(), score.row(19).mean(), score.row(20).mean(), - score.row(24).mean(), 1.; + std::cout << "\nTrue parameters: " + << "phi: " << true_par(0) << " " + << "dxy: " << true_par(1) << " " + << "pt: " << true_par(2) << " " + << "CotT: " << true_par(3) << " " + << "Zip: " << true_par(4) << " " + << std::endl; + for (int i = 0; i < iteration; i++) { + hits_gen gen; + gen = Hits_gen(n_, gen_par); + // gen.hits = MatrixXd::Zero(3, 4); + // gen.hits_cov = MatrixXd::Zero(3 * 4, 3 * 4); + // gen.hits.col(0) << 1.82917642593, 2.0411875248, 7.18495464325; + // gen.hits.col(1) << 4.47041416168, 4.82704305649, 18.6394691467; + // gen.hits.col(2) << 7.25991010666, 7.74653434753, 30.6931324005; + // gen.hits.col(3) << 8.99161434174, 9.54262828827, 38.1338043213; + helixRiemann_fit[i] = Rfit::Helix_fit(gen.hits, gen.hits_cov, B_field, return_err); +// helixBrokenLine_fit[i] = BrokenLine::Helix_fit(gen.hits, gen.hits_cov, B_field); - cout << "\nPULLS:\n" - << "phi: " << phi_ << " " - << sqrt((score.row(0).array() - phi_).square().sum() / (iteration - 1)) << " " - << abs(score.row(10).mean()) << "%\n" - << "a0 : " << a_ << " " - << sqrt((score.row(1).array() - a_).square().sum() / (iteration - 1)) << " " - << abs(score.row(11).mean()) << "%\n" - << "pt : " << pt_ << " " - << sqrt((score.row(2).array() - pt_).square().sum() / (iteration - 1)) << " " - << abs(score.row(12).mean()) << "%\n" - << "coT: " << coT_ << " " - << sqrt((score.row(3).array() - coT_).square().sum() / (iteration - 1)) << " " - << abs(score.row(13).mean()) << "%\n" - << "Zip: " << Zip_ << " " - << sqrt((score.row(4).array() - Zip_).square().sum() / (iteration - 1)) << " " - << abs(score.row(14).mean()) << "%\n\n" - << "cov(phi,a0)_: " << score.row(5).mean() << "\n" - << "cov(phi,pt)_: " << score.row(6).mean() << "\n" - << "cov(a0,pt)_: " << score.row(7).mean() << "\n" - << "cov(coT,Zip)_: " << score.row(8).mean() << "\n\n" - << "chi2_circle: " << score.row(9).mean() << " vs " << n_ - 3 << "\n" - << "chi2_line: " << score.row(25).mean() << " vs " << n_ - 2 << "\n\n" - << "correlation matrix:\n" - << correlation << "\n\n" - << endl; -// } + std::cout << std::endl; + /* + if (debug) + cout << std::setprecision(6) + << "phi: " << helixRiemann_fit[i].par(0) << " +/- " << sqrt(helixRiemann_fit[i].cov(0, 0)) << " vs " + << true_par(0) << endl + << "Tip: " << helixRiemann_fit[i].par(1) << " +/- " << sqrt(helixRiemann_fit[i].cov(1, 1)) << " vs " + << true_par(1) << endl + << "p_t: " << helixRiemann_fit[i].par(2) << " +/- " << sqrt(helixRiemann_fit[i].cov(2, 2)) << " vs " + << true_par(2) << endl + << "theta:" << helixRiemann_fit[i].par(3) << " +/- " << sqrt(helixRiemann_fit[i].cov(3, 3)) << " vs " + << true_par(3) << endl + << "Zip: " << helixRiemann_fit[i].par(4) << " +/- " << sqrt(helixRiemann_fit[i].cov(4, 4)) << " vs " + << true_par(4) << endl + << "charge:" << helixRiemann_fit[i].q << " vs 1" << endl + << "covariance matrix:" << endl + << helixRiemann_fit[i].cov << endl + << "Initial hits:\n" << gen.hits << endl + << "Initial Covariance:\n" << gen.hits_cov << endl; + */ + } + computePull(helixRiemann_fit, "Riemann", n_, iteration, true_par); +// computePull(helixBrokenLine_fit, "BrokenLine", n_, iteration, true_par); } int main() { + TFile f("TestFitResults.root", "RECREATE"); test_helix_fit(); + f.Close(); return 0; } diff --git a/RecoPixelVertexing/PixelTrackFitting/test/testEigenGPU.cu b/RecoPixelVertexing/PixelTrackFitting/test/testEigenGPU.cu index 7b1125eebc312..485fac34b00b2 100644 --- a/RecoPixelVertexing/PixelTrackFitting/test/testEigenGPU.cu +++ b/RecoPixelVertexing/PixelTrackFitting/test/testEigenGPU.cu @@ -29,12 +29,10 @@ void kernelFullFit(Rfit::Matrix3xNd * hits, Rfit::Matrix2Nd hits_cov2D_local = (hits_cov->block(0, 0, 2 * n, 2 * n)).eval(); Rfit::printIt(&hits2D_local, "kernelFullFit - hits2D_local: "); Rfit::printIt(&hits_cov2D_local, "kernelFullFit - hits_cov2D_local: "); - /* printf("kernelFullFit - hits address: %p\n", hits); printf("kernelFullFit - hits_cov address: %p\n", hits_cov); printf("kernelFullFit - hits_cov2D address: %p\n", &hits2D_local); printf("kernelFullFit - hits_cov2D_local address: %p\n", &hits_cov2D_local); - */ /* At some point I gave up and locally construct block on the stack, so that the next invocation to Rfit::Circle_fit works properly. Failing to do so implied basically an empty collection of hits and covariances. That could @@ -43,60 +41,20 @@ void kernelFullFit(Rfit::Matrix3xNd * hits, creations of the blocks. To be understood and compared against the myriad of compilation warnings we have. */ + (*circle_fit_resultsGPU) = Rfit::Circle_fit(hits->block(0,0,2,n), hits_cov->block(0, 0, 2 * n, 2 * n), - fast_fit, rad, B, errors); + fast_fit, rad, B, errors); /* - (*circle_fit_resultsGPU) = - Rfit::Circle_fit(hits2D_local, hits_cov2D_local, - fast_fit, rad, B, errors, scattering); - */ - (*line_fit_resultsGPU) = Rfit::Line_fit(*hits, *hits_cov, *circle_fit_resultsGPU, fast_fit, errors); + (*circle_fit_resultsGPU) = + Rfit::Circle_fit(hits2D_local, hits_cov2D_local, + fast_fit, rad, B, errors); + */ + (*line_fit_resultsGPU) = Rfit::Line_fit(*hits, *hits_cov, *circle_fit_resultsGPU, fast_fit, B, errors); return; } -__global__ -void kernelFastFit(Rfit::Matrix3xNd * hits, Vector4d * results) { - (*results) = Rfit::Fast_fit(*hits); -} - -__global__ -void kernelCircleFit(Rfit::Matrix3xNd * hits, - Rfit::Matrix3Nd * hits_cov, Vector4d * fast_fit_input, double B, - Rfit::circle_fit * circle_fit_resultsGPU) { - u_int n = hits->cols(); - Rfit::VectorNd rad = (hits->block(0, 0, 2, n).colwise().norm()); - -#if TEST_DEBUG - printf("fast_fit_input(0): %f\n", (*fast_fit_input)(0)); - printf("fast_fit_input(1): %f\n", (*fast_fit_input)(1)); - printf("fast_fit_input(2): %f\n", (*fast_fit_input)(2)); - printf("fast_fit_input(3): %f\n", (*fast_fit_input)(3)); - printf("rad(0,0): %f\n", rad(0,0)); - printf("rad(1,1): %f\n", rad(1,1)); - printf("rad(2,2): %f\n", rad(2,2)); - printf("hits_cov(0,0): %f\n", (*hits_cov)(0,0)); - printf("hits_cov(1,1): %f\n", (*hits_cov)(1,1)); - printf("hits_cov(2,2): %f\n", (*hits_cov)(2,2)); - printf("hits_cov(11,11): %f\n", (*hits_cov)(11,11)); - printf("B: %f\n", B); -#endif - (*circle_fit_resultsGPU) = - Rfit::Circle_fit(hits->block(0,0,2,n), hits_cov->block(0, 0, 2 * n, 2 * n), - *fast_fit_input, rad, B, false); -} - -__global__ -void kernelLineFit(Rfit::Matrix3xNd * hits, - Rfit::Matrix3Nd * hits_cov, - Rfit::circle_fit * circle_fit, - Vector4d * fast_fit, - Rfit::line_fit * line_fit) -{ - (*line_fit) = Rfit::Line_fit(*hits, *hits_cov, *circle_fit, *fast_fit, true); -} - void fillHitsAndHitsCov(Rfit::Matrix3xNd & hits, Rfit::Matrix3Nd & hits_cov) { hits << 1.98645, 4.72598, 7.65632, 11.3151, 2.18002, 4.88864, 7.75845, 11.3134, @@ -119,79 +77,6 @@ void fillHitsAndHitsCov(Rfit::Matrix3xNd & hits, Rfit::Matrix3Nd & hits_cov) { hits_cov(3,7) = hits_cov(7,3) = -5.28e-06; } -void testFit() { - constexpr double B = 0.0113921; - Rfit::Matrix3xNd hits(3,4); - Rfit::Matrix3Nd hits_cov = MatrixXd::Zero(12,12); - Rfit::Matrix3xNd * hitsGPU = new Rfit::Matrix3xNd(3,4); - Rfit::Matrix3Nd * hits_covGPU = nullptr; - Vector4d * fast_fit_resultsGPU = new Vector4d(); - Vector4d * fast_fit_resultsGPUret = new Vector4d(); - Rfit::circle_fit * circle_fit_resultsGPU = new Rfit::circle_fit(); - Rfit::circle_fit * circle_fit_resultsGPUret = new Rfit::circle_fit(); - - fillHitsAndHitsCov(hits, hits_cov); - - // FAST_FIT_CPU - Vector4d fast_fit_results = Rfit::Fast_fit(hits); -#if TEST_DEBUG - std::cout << "Generated hits:\n" << hits << std::endl; -#endif - std::cout << "Fitted values (FastFit, [X0, Y0, R, tan(theta)]):\n" << fast_fit_results << std::endl; - - // FAST_FIT GPU - cudaMalloc((void**)&hitsGPU, sizeof(Rfit::Matrix3xNd(3,4))); - cudaMalloc((void**)&fast_fit_resultsGPU, sizeof(Vector4d)); - cudaMemcpy(hitsGPU, &hits, sizeof(Rfit::Matrix3xNd(3,4)), cudaMemcpyHostToDevice); - - kernelFastFit<<<1, 1>>>(hitsGPU, fast_fit_resultsGPU); - cudaDeviceSynchronize(); - - cudaMemcpy(fast_fit_resultsGPUret, fast_fit_resultsGPU, sizeof(Vector4d), cudaMemcpyDeviceToHost); - std::cout << "Fitted values (FastFit, [X0, Y0, R, tan(theta)]): GPU\n" << *fast_fit_resultsGPUret << std::endl; - assert(isEqualFuzzy(fast_fit_results, (*fast_fit_resultsGPUret))); - - // CIRCLE_FIT CPU - u_int n = hits.cols(); - Rfit::VectorNd rad = (hits.block(0, 0, 2, n).colwise().norm()); - - Rfit::circle_fit circle_fit_results = Rfit::Circle_fit(hits.block(0, 0, 2, n), - hits_cov.block(0, 0, 2 * n, 2 * n), - fast_fit_results, rad, B, false); - std::cout << "Fitted values (CircleFit):\n" << circle_fit_results.par << std::endl; - - // CIRCLE_FIT GPU - cudaMalloc((void **)&hits_covGPU, sizeof(Rfit::Matrix3Nd(12,12))); - cudaMalloc((void **)&circle_fit_resultsGPU, sizeof(Rfit::circle_fit)); - cudaMemcpy(hits_covGPU, &hits_cov, sizeof(Rfit::Matrix3Nd(12,12)), cudaMemcpyHostToDevice); - - kernelCircleFit<<<1,1>>>(hitsGPU, hits_covGPU, - fast_fit_resultsGPU, B, circle_fit_resultsGPU); - cudaDeviceSynchronize(); - - cudaMemcpy(circle_fit_resultsGPUret, circle_fit_resultsGPU, - sizeof(Rfit::circle_fit), cudaMemcpyDeviceToHost); - std::cout << "Fitted values (CircleFit) GPU:\n" << circle_fit_resultsGPUret->par << std::endl; - assert(isEqualFuzzy(circle_fit_results.par, circle_fit_resultsGPUret->par)); - - // LINE_FIT CPU - Rfit::line_fit line_fit_results = Rfit::Line_fit(hits, hits_cov, circle_fit_results, fast_fit_results, true); - std::cout << "Fitted values (LineFit):\n" << line_fit_results.par << std::endl; - - // LINE_FIT GPU - Rfit::line_fit * line_fit_resultsGPU = nullptr; - Rfit::line_fit * line_fit_resultsGPUret = new Rfit::line_fit(); - - cudaMalloc((void **)&line_fit_resultsGPU, sizeof(Rfit::line_fit)); - - kernelLineFit<<<1,1>>>(hitsGPU, hits_covGPU, circle_fit_resultsGPU, fast_fit_resultsGPU, line_fit_resultsGPU); - cudaDeviceSynchronize(); - - cudaMemcpy(line_fit_resultsGPUret, line_fit_resultsGPU, sizeof(Rfit::line_fit), cudaMemcpyDeviceToHost); - std::cout << "Fitted values (LineFit) GPU:\n" << line_fit_resultsGPUret->par << std::endl; - assert(isEqualFuzzy(line_fit_results.par, line_fit_resultsGPUret->par)); -} - void testFitOneGo(bool errors, double epsilon=1e-6) { constexpr double B = 0.0113921; Rfit::Matrix3xNd hits(3,4); @@ -205,12 +90,12 @@ void testFitOneGo(bool errors, double epsilon=1e-6) { u_int n = hits.cols(); Rfit::VectorNd rad = (hits.block(0, 0, 2, n).colwise().norm()); - Rfit::circle_fit circle_fit_results = Rfit::Circle_fit(hits.block(0, 0, 2, n), + Rfit::circle_fit circle_fit_results = Rfit::Circle_fit(hits.block(0, 0, 2, n), hits_cov.block(0, 0, 2 * n, 2 * n), fast_fit_results, rad, B, errors); // LINE_FIT CPU Rfit::line_fit line_fit_results = Rfit::Line_fit(hits, hits_cov, circle_fit_results, - fast_fit_results, errors); + fast_fit_results, B, errors); // FIT GPU std::cout << "GPU FIT" << std::endl; @@ -253,10 +138,8 @@ void testFitOneGo(bool errors, double epsilon=1e-6) { } int main (int argc, char * argv[]) { -// testFit(); - std::cout << "TEST FIT, NO ERRORS" << std::endl; - testFitOneGo(false); + cudaDeviceSetLimit(cudaLimitStackSize, 32*1024); std::cout << "TEST FIT, ERRORS AND SCATTER" << std::endl; testFitOneGo(true, 1e-5);