File indexing completed on 2024-04-06 12:28:40
0001 #include <iostream>
0002
0003 #include <Eigen/Core>
0004 #include <Eigen/Eigenvalues>
0005
0006 #ifdef USE_BL
0007 #include "RecoTracker/PixelTrackFitting/interface/BrokenLine.h"
0008 #else
0009 #include "RecoTracker/PixelTrackFitting/interface/RiemannFit.h"
0010 #endif
0011
0012 #include "test_common.h"
0013
0014 using namespace Eigen;
0015
0016 namespace riemannFit {
0017 constexpr uint32_t maxNumberOfTracks() { return 5 * 1024; }
0018 constexpr uint32_t stride() { return maxNumberOfTracks(); }
0019
0020 template <int N>
0021 using Matrix3xNd = Eigen::Matrix<double, 3, N>;
0022 template <int N>
0023 using Map3xNd = Eigen::Map<Matrix3xNd<N>, 0, Eigen::Stride<3 * stride(), stride()> >;
0024
0025 template <int N>
0026 using Matrix6xNf = Eigen::Matrix<float, 6, N>;
0027 template <int N>
0028 using Map6xNf = Eigen::Map<Matrix6xNf<N>, 0, Eigen::Stride<6 * stride(), stride()> >;
0029
0030 using Map4d = Eigen::Map<Vector4d, 0, Eigen::InnerStride<stride()> >;
0031
0032 }
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047 template <typename M3xN, typename M6xN>
0048 void fillHitsAndHitsCov(M3xN& hits, M6xN& hits_ge) {
0049 constexpr uint32_t N = M3xN::ColsAtCompileTime;
0050
0051 if (N == 5) {
0052 hits << 2.934787, 6.314229, 8.936963, 10.360559, 12.856387, 0.773211, 1.816356, 2.765734, 3.330824, 4.422212,
0053 -10.980247, -23.162731, -32.759060, -38.061260, -47.518867;
0054 hits_ge.col(0) << 1.424715e-07, -4.996975e-07, 1.752614e-06, 3.660689e-11, 1.644638e-09, 7.346080e-05;
0055 hits_ge.col(1) << 6.899177e-08, -1.873414e-07, 5.087101e-07, -2.078806e-10, -2.210498e-11, 4.346079e-06;
0056 hits_ge.col(2) << 1.406273e-06, 4.042467e-07, 6.391180e-07, -3.141497e-07, 6.513821e-08, 1.163863e-07;
0057 hits_ge.col(3) << 1.176358e-06, 2.154100e-07, 5.072816e-07, -8.161219e-08, 1.437878e-07, 5.951832e-08;
0058 hits_ge.col(4) << 2.852843e-05, 7.956492e-06, 3.117701e-06, -1.060541e-06, 8.777413e-09, 1.426417e-07;
0059 return;
0060 }
0061
0062 if (N > 3)
0063 hits << 1.98645, 4.72598, 7.65632, 11.3151, 2.18002, 4.88864, 7.75845, 11.3134, 2.46338, 6.99838, 11.808, 17.793;
0064 else
0065 hits << 1.98645, 4.72598, 7.65632, 2.18002, 4.88864, 7.75845, 2.46338, 6.99838, 11.808;
0066
0067 hits_ge.col(0)[0] = 7.14652e-06;
0068 hits_ge.col(1)[0] = 2.15789e-06;
0069 hits_ge.col(2)[0] = 1.63328e-06;
0070 if (N > 3)
0071 hits_ge.col(3)[0] = 6.27919e-06;
0072 hits_ge.col(0)[2] = 6.10348e-06;
0073 hits_ge.col(1)[2] = 2.08211e-06;
0074 hits_ge.col(2)[2] = 1.61672e-06;
0075 if (N > 3)
0076 hits_ge.col(3)[2] = 6.28081e-06;
0077 hits_ge.col(0)[5] = 5.184e-05;
0078 hits_ge.col(1)[5] = 1.444e-05;
0079 hits_ge.col(2)[5] = 6.25e-06;
0080 if (N > 3)
0081 hits_ge.col(3)[5] = 3.136e-05;
0082 hits_ge.col(0)[1] = -5.60077e-06;
0083 hits_ge.col(1)[1] = -1.11936e-06;
0084 hits_ge.col(2)[1] = -6.24945e-07;
0085 if (N > 3)
0086 hits_ge.col(3)[1] = -5.28e-06;
0087 }
0088
0089 template <int N>
0090 void testFit() {
0091 constexpr double B = 0.0113921;
0092 riemannFit::Matrix3xNd<N> hits;
0093 riemannFit::Matrix6xNf<N> hits_ge = MatrixXf::Zero(6, N);
0094
0095 fillHitsAndHitsCov(hits, hits_ge);
0096
0097 std::cout << "sizes " << N << ' ' << sizeof(hits) << ' ' << sizeof(hits_ge) << ' ' << sizeof(Vector4d) << std::endl;
0098
0099 std::cout << "Generated hits:\n" << hits << std::endl;
0100 std::cout << "Generated cov:\n" << hits_ge << std::endl;
0101
0102
0103 #ifdef USE_BL
0104 Vector4d fast_fit_results;
0105 brokenline::fastFit(hits, fast_fit_results);
0106 #else
0107 Vector4d fast_fit_results;
0108 riemannFit::fastFit(hits, fast_fit_results);
0109 #endif
0110 std::cout << "Fitted values (FastFit, [X0, Y0, R, tan(theta)]):\n" << fast_fit_results << std::endl;
0111
0112
0113
0114 #ifdef USE_BL
0115 brokenline::PreparedBrokenLineData<N> data;
0116 brokenline::karimaki_circle_fit circle_fit_results;
0117 riemannFit::Matrix3d Jacob;
0118
0119 brokenline::prepareBrokenLineData(hits, fast_fit_results, B, data);
0120 riemannFit::LineFit line_fit_results;
0121 brokenline::lineFit(hits_ge, fast_fit_results, B, data, line_fit_results);
0122 brokenline::circleFit(hits, hits_ge, fast_fit_results, B, data, circle_fit_results);
0123 Jacob << 1., 0, 0, 0, 1., 0, 0, 0,
0124 -B / std::copysign(riemannFit::sqr(circle_fit_results.par(2)), circle_fit_results.par(2));
0125 circle_fit_results.par(2) = B / std::abs(circle_fit_results.par(2));
0126 circle_fit_results.cov = Jacob * circle_fit_results.cov * Jacob.transpose();
0127 #else
0128 riemannFit::VectorNd<N> rad = (hits.block(0, 0, 2, N).colwise().norm());
0129 riemannFit::Matrix2Nd<N> hits_cov = riemannFit::Matrix2Nd<N>::Zero();
0130 riemannFit::loadCovariance2D(hits_ge, hits_cov);
0131 riemannFit::CircleFit circle_fit_results =
0132 riemannFit::circleFit(hits.block(0, 0, 2, N), hits_cov, fast_fit_results, rad, B, true);
0133
0134 riemannFit::LineFit line_fit_results =
0135 riemannFit::lineFit(hits, hits_ge, circle_fit_results, fast_fit_results, B, true);
0136 riemannFit::par_uvrtopak(circle_fit_results, B, true);
0137
0138 #endif
0139
0140 std::cout << "Fitted values (CircleFit):\n"
0141 << circle_fit_results.par << "\nchi2 " << circle_fit_results.chi2 << std::endl;
0142 std::cout << "Fitted values (LineFit):\n" << line_fit_results.par << "\nchi2 " << line_fit_results.chi2 << std::endl;
0143
0144 std::cout << "Fitted cov (CircleFit) CPU:\n" << circle_fit_results.cov << std::endl;
0145 std::cout << "Fitted cov (LineFit): CPU\n" << line_fit_results.cov << std::endl;
0146 }
0147
0148 int main(int argc, char* argv[]) {
0149 testFit<4>();
0150 testFit<3>();
0151 testFit<5>();
0152
0153 return 0;
0154 }