PixelTrackFitting/test/testFits.cpp

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   // hits
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   // errors
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   // fast fit
0030   using Map4d = Eigen::Map<Vector4d, 0, Eigen::InnerStride<stride()> >;
0031
0032 }  // namespace riemannFit
0033
0034 /*
0035 Hit global: 641,0 2: 2.934787,0.773211,-10.980247
0036 Error: 641,0 2: 1.424715e-07,-4.996975e-07,1.752614e-06,3.660689e-11,1.644638e-09,7.346080e-05
0037 Hit global: 641,1 104: 6.314229,1.816356,-23.162731
0038 Error: 641,1 104: 6.899177e-08,-1.873414e-07,5.087101e-07,-2.078806e-10,-2.210498e-11,4.346079e-06
0039 Hit global: 641,2 1521: 8.936963,2.765734,-32.759060
0040 Error: 641,2 1521: 1.406273e-06,4.042467e-07,6.391180e-07,-3.141497e-07,6.513821e-08,1.163863e-07
0041 Hit global: 641,3 1712: 10.360559,3.330824,-38.061260
0042 Error: 641,3 1712: 1.176358e-06,2.154100e-07,5.072816e-07,-8.161219e-08,1.437878e-07,5.951832e-08
0043 Hit global: 641,4 1824: 12.856387,4.422212,-47.518867
0044 Error: 641,4 1824: 2.852843e-05,7.956492e-06,3.117701e-06,-1.060541e-06,8.777413e-09,1.426417e-07
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   // FAST_FIT_CPU
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   // CIRCLE_FIT CPU
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   // LINE_FIT CPU
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 }