From a0f04b54a13ac087876457e5721e1788ab7c5fb0 Mon Sep 17 00:00:00 2001 From: Andrea Bocci Date: Wed, 14 Nov 2018 23:56:26 +0100 Subject: [PATCH] Synchronise with CMSSW_10_4_0_pre2 --- .../PixelTrackFitting/interface/RiemannFit.h | 510 +++++++----------- .../test/PixelTrackRiemannFit.cc | 388 +++++-------- .../PixelTrackFitting/test/testEigenGPU.cu | 137 ++++- .../PixelTriplets/plugins/BuildFile.xml | 3 +- .../python/plotting/trackingPlots.py | 16 +- 5 files changed, 469 insertions(+), 585 deletions(-) diff --git a/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h b/RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h index 33f8334c8b5a5..6de1c77bbac12 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,6 +78,7 @@ 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 @@ -118,40 +119,31 @@ __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 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; - } +__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.); } /*! @@ -159,29 +151,11 @@ void computeRadLenUniformMaterial(const VectorNd &length_values, 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 @@ -190,24 +164,47 @@ __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))); - 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(); + 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); + } + } + } + } #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) = sin(theta); - rot(n + i, n + i) = sin(theta); + rot(i, i) = cos(theta); + rot(n + i, n + i) = cos(theta); u_int j = (i + n); - rot(i, j) = i < j ? cos(theta) : -cos(theta); + // 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); } #if RFIT_DEBUG @@ -215,23 +212,12 @@ __host__ __device__ inline MatrixNd Scatter_cov_line(Matrix2Nd& cov_sz, #endif Matrix2Nd tmp = rot*cov_sz*rot.transpose(); - 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. + // We are interested only in the errors in the rotated s -axis which, in + // our formalism, are in the upper square matrix. #if RFIT_DEBUG Rfit::printIt(&tmp, "Scatter_cov_line - tmp: "); #endif - return tmp.block(n, n, n, n); + return tmp.block(0, 0, n, n); } /*! @@ -247,11 +233,13 @@ __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, @@ -260,32 +248,24 @@ __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 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)); + double radlen_eff = 0.; + double theta = 0.; + bool in_forward = false; + computeRadLenEff(fast_fit, B, radlen_eff, theta, in_forward); - // 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); - VectorNd sig2(n); - sig2 = .000225 / p_2 * (1.f + 0.038 * rad_lengths.array().log()).abs2() * rad_lengths.array(); + 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) { - scatter_cov_rad(k, l) += (rad(k) - rad(i)) * (rad(l) - rad(i)) * sig2(i) / (sqr(sin(theta))); + 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(l, k) = scatter_cov_rad(k, l); } } @@ -429,6 +409,23 @@ __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 @@ -473,6 +470,40 @@ __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. @@ -973,7 +1004,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, - v(0)*v2x2_inv*t, v(1)*v2x2_inv*t, -h * sqr(v2x2_inv) * 2. - (2. * c + v(2)) * v2x2_inv * t, -t; + 0, 0, -h * sqr(v2x2_inv) * 2. - (2. * c + v(2)) * v2x2_inv * t, -t; } printIt(&J3, "circle_fit - J3:"); @@ -1028,22 +1059,21 @@ __host__ __device__ inline circle_fit Circle_fit(const Matrix2xNd& hits2D, errors. */ -__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) +__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) { 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 = MatrixXd::Zero(2, n); - Eigen::Matrix Jx; + Matrix2xNd p2D(2, n); + MatrixNx5d Jx(n, 5); #if RFIT_DEBUG printf("Line_fit - B: %g\n", B); + printIt(&hits, "Line_fit points: "); printIt(&hits_cov, "Line_fit covs: "); #endif @@ -1055,11 +1085,8 @@ __host__ __device__ inline line_fit Line_fit_odr(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); @@ -1068,9 +1095,9 @@ __host__ __device__ inline line_fit Line_fit_odr(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); @@ -1079,19 +1106,7 @@ __host__ __device__ inline line_fit Line_fit_odr(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 << 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); + Jx.row(i) << d_X0, d_Y0, d_R, d_x, d_y; } // Math of d_{X0,Y0,R,x,y} all verified by hand @@ -1099,26 +1114,44 @@ __host__ __device__ inline line_fit Line_fit_odr(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), theta, B); + MatrixNd cov_with_ms = Scatter_cov_line(cov_sz, fast_fit, p2D.row(0), p2D.row(1), B); #if RFIT_DEBUG printIt(&cov_with_ms, "line_fit - cov_with_ms: "); #endif - Matrix4d G; - G = cov_with_ms.inverse(); + Matrix4d G, G4; + G4 = cov_with_ms.inverse(); #if RFIT_DEBUG - printIt(&G, "line_fit - cov_with_ms.inverse():"); + printIt(&G4, "line_fit - cov_with_ms.inverse():"); #endif - double renorm = G.sum(); - G *= 1. / renorm; + double renorm = G4.sum(); + G4 *= 1. / renorm; #if RFIT_DEBUG - printIt(&G, "line_fit - G4:"); + printIt(&G4, "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 + // COST FUNCTION // compute @@ -1129,12 +1162,16 @@ __host__ __device__ inline line_fit Line_fit_odr(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. v is normalized!! + // minimize double chi2; Vector2d v = min_eigen2D(A, chi2); #if RFIT_DEBUG @@ -1142,6 +1179,7 @@ __host__ __device__ inline line_fit Line_fit_odr(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; @@ -1151,8 +1189,8 @@ __host__ __device__ inline line_fit Line_fit_odr(const Matrix3xNd& hits, // COMPUTE LINE PARAMETER line_fit line; line.par << -v(0) / v(1), // cotan(theta)) - -c / v(1); // Zip - line.chi2 = abs(chi2*renorm); + -c * sqrt(sqr(v(0)) + sqr(v(1))) * 1. / v(1); // Zip + line.chi2 = abs(chi2); #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))); @@ -1168,21 +1206,19 @@ __host__ __device__ inline line_fit Line_fit_odr(const Matrix3xNd& hits, { // The norm is taken from Chernov, properly adapted to the weights case. double norm = v.transpose() * A * v; -// double norm_empirical = cov_with_ms.diagonal().mean(); + norm /= weight.sum(); #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 = norm/(A(0,0) + A(1,1)); + const double sig2 = 1. / (A(0, 0) + A(1, 1)) * norm; C(0, 0) = sig2 * v1_2; C(1, 1) = sig2 * v0_2; - 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); + 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); } #if RFIT_DEBUG printIt(&C, "line_fit - C:"); @@ -1192,7 +1228,9 @@ __host__ __device__ inline line_fit Line_fit_odr(const Matrix3xNd& hits, { const double t0 = 1. / v(1); const double t1 = sqr(t0); - J << -t0, v(0) * t1, 0., 0., c * t1, -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; } Matrix JT = J.transpose().eval(); #if RFIT_DEBUG @@ -1207,184 +1245,6 @@ __host__ __device__ inline line_fit Line_fit_odr(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 f71e3f082ada4..b27ed52473388 100644 --- a/RecoPixelVertexing/PixelTrackFitting/test/PixelTrackRiemannFit.cc +++ b/RecoPixelVertexing/PixelTrackFitting/test/PixelTrackRiemannFit.cc @@ -6,11 +6,7 @@ #include #include // unique_ptr -#include -#include - #include "RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h" -//#include "RecoPixelVertexing/PixelTrackFitting/interface/BrokenLine.h" using namespace std; using namespace Eigen; @@ -168,261 +164,161 @@ 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; - 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; +// 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; - 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; + 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); - 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); + 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); - 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); + 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)); + } + + 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.; + + 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; +// } } 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 485fac34b00b2..7b1125eebc312 100644 --- a/RecoPixelVertexing/PixelTrackFitting/test/testEigenGPU.cu +++ b/RecoPixelVertexing/PixelTrackFitting/test/testEigenGPU.cu @@ -29,10 +29,12 @@ 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 @@ -41,20 +43,60 @@ 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); - */ - (*line_fit_resultsGPU) = Rfit::Line_fit(*hits, *hits_cov, *circle_fit_resultsGPU, fast_fit, 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); 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, @@ -77,25 +119,98 @@ void fillHitsAndHitsCov(Rfit::Matrix3xNd & hits, Rfit::Matrix3Nd & hits_cov) { hits_cov(3,7) = hits_cov(7,3) = -5.28e-06; } -void testFitOneGo(bool errors, double epsilon=1e-6) { +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); + Rfit::Matrix3Nd hits_cov = MatrixXd::Zero(12,12); + + fillHitsAndHitsCov(hits, hits_cov); + + // FAST_FIT_CPU + Vector4d fast_fit_results = Rfit::Fast_fit(hits); + // 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, errors); // LINE_FIT CPU Rfit::line_fit line_fit_results = Rfit::Line_fit(hits, hits_cov, circle_fit_results, - fast_fit_results, B, errors); + fast_fit_results, errors); // FIT GPU std::cout << "GPU FIT" << std::endl; @@ -138,8 +253,10 @@ 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); diff --git a/RecoPixelVertexing/PixelTriplets/plugins/BuildFile.xml b/RecoPixelVertexing/PixelTriplets/plugins/BuildFile.xml index 8956caca42899..77cb6c4da68a4 100644 --- a/RecoPixelVertexing/PixelTriplets/plugins/BuildFile.xml +++ b/RecoPixelVertexing/PixelTriplets/plugins/BuildFile.xml @@ -1,5 +1,6 @@ + @@ -10,7 +11,7 @@ + - diff --git a/Validation/RecoTrack/python/plotting/trackingPlots.py b/Validation/RecoTrack/python/plotting/trackingPlots.py index ac81473c843cb..f14650a8d92ff 100644 --- a/Validation/RecoTrack/python/plotting/trackingPlots.py +++ b/Validation/RecoTrack/python/plotting/trackingPlots.py @@ -1223,6 +1223,12 @@ def _trackingFolders(lastDirName="Track"): + _makeMVAPlots(3) \ + _makeMVAPlots(3, hp=True) # add more if needed +_buildingExtendedPlots = [ + _pulls, + _resolutionsEta, + _resolutionsPt, + _tuning, +] _extendedPlots = [ _extDistPtEtaPhi, _extDistDxyDzBS, @@ -1279,7 +1285,7 @@ def _trackingFolders(lastDirName="Track"): ] plotter = Plotter() plotterExt = Plotter() -def _appendTrackingPlots(lastDirName, name, algoPlots, onlyForPileup=False, onlyForElectron=False, onlyForConversion=False, onlyForBHadron=False, seeding=False, rawSummary=False, highPuritySummary=True): +def _appendTrackingPlots(lastDirName, name, algoPlots, onlyForPileup=False, onlyForElectron=False, onlyForConversion=False, onlyForBHadron=False, seeding=False, building=False, rawSummary=False, highPuritySummary=True): folders = _trackingFolders(lastDirName) # to keep backward compatibility, this set of plots has empty name limiters = dict(onlyForPileup=onlyForPileup, onlyForElectron=onlyForElectron, onlyForConversion=onlyForConversion, onlyForBHadron=onlyForBHadron) @@ -1290,7 +1296,11 @@ def _appendTrackingPlots(lastDirName, name, algoPlots, onlyForPileup=False, only _trackingSubFoldersFallbackSLHC_Phase1PU140, _trackingSubFoldersFallbackFromPV, _trackingSubFoldersFallbackConversion]) plotter.append(name, folders, TrackingPlotFolder(*algoPlots, **commonForTPF), **common) - plotterExt.append(name, folders, TrackingPlotFolder(*_extendedPlots, **commonForTPF), **common) + extendedPlots = [] + if building: + extendedPlots.extend(_buildingExtendedPlots) + extendedPlots.extend(_extendedPlots) + plotterExt.append(name, folders, TrackingPlotFolder(*extendedPlots, **commonForTPF), **common) summaryName = "" if name != "": @@ -1330,7 +1340,7 @@ def _appendTrackingPlots(lastDirName, name, algoPlots, onlyForPileup=False, only _appendTrackingPlots("TrackFromPVAllTP", "fromPVAllTP", _simBasedPlots+_recoBasedPlots, onlyForPileup=True) _appendTrackingPlots("TrackFromPVAllTP2", "fromPVAllTP2", _simBasedPlots+_recoBasedPlots, onlyForPileup=True) _appendTrackingPlots("TrackSeeding", "seeding", _seedingBuildingPlots, seeding=True) -_appendTrackingPlots("TrackBuilding", "building", _seedingBuildingPlots) +_appendTrackingPlots("TrackBuilding", "building", _seedingBuildingPlots, building=True) _appendTrackingPlots("TrackConversion", "conversion", _simBasedPlots+_recoBasedPlots, onlyForConversion=True, rawSummary=True, highPuritySummary=False) _appendTrackingPlots("TrackGsf", "gsf", _simBasedPlots+_recoBasedPlots, onlyForElectron=True, rawSummary=True, highPuritySummary=False) _appendTrackingPlots("TrackBHadron", "bhadron", _simBasedPlots+_recoBasedPlots, onlyForBHadron=True)