PixelTrackFitting/test/PixelTrackFits.cc

0001 #define _USE_MATH_DEFINES
0002
0003 #include <chrono>
0004 #include <cmath>
0005 #include <iomanip>
0006 #include <iostream>
0007 #include <memory>
0008 #include <random>
0009
0010 #include <TFile.h>
0011 #include <TH1F.h>
0012
0013 #ifdef USE_BL
0014 #include "RecoTracker/PixelTrackFitting/interface/BrokenLine.h"
0015 #else
0016 #include "RecoTracker/PixelTrackFitting/interface/RiemannFit.h"
0017 #endif
0018
0019 using namespace std;
0020 using namespace Eigen;
0021 using namespace riemannFit;
0022 using std::unique_ptr;
0023
0024 namespace riemannFit {
0025   using Vector3i = Eigen::Matrix<int, 3, 1>;
0026   using Vector4i = Eigen::Matrix<int, 4, 1>;
0027   using Vector6d = Eigen::Matrix<double, 6, 1>;
0028   using Vector8d = Eigen::Matrix<double, 8, 1>;
0029 };  // namespace riemannFit
0030
0031 // quadruplets...
0032 struct hits_gen {
0033   Matrix3xNd<4> hits;
0034   Eigen::Matrix<float, 6, 4> hits_ge;
0035   Vector5d true_par;
0036 };
0037
0038 struct geometry {
0039   Vector8d barrel;
0040   Vector4i barrel_2;
0041   Vector8d R_err;
0042   Vector8d Rp_err;
0043   Vector8d z_err;
0044   Vector6d hand;
0045   Vector3i hand_2;
0046   Vector6d xy_err;
0047   Vector6d zh_err;
0048   double z_max;
0049   double r_max;
0050 };
0051
0052 void test_helix_fit();
0053
0054 constexpr int c_speed = 299792458;
0055 constexpr double pi = M_PI;
0056 default_random_engine generator(1);
0057
0058 void smearing(const Vector5d& err, const bool& isbarrel, double& x, double& y, double& z) {
0059   normal_distribution<double> dist_R(0., err[0]);
0060   normal_distribution<double> dist_Rp(0., err[1]);
0061   normal_distribution<double> dist_z(0., err[2]);
0062   normal_distribution<double> dist_xyh(0., err[3]);
0063   normal_distribution<double> dist_zh(0., err[4]);
0064   if (isbarrel) {
0065     double dev_Rp = dist_Rp(generator);
0066     double dev_R = dist_R(generator);
0067     double R = sqrt(riemannFit::sqr(x) + riemannFit::sqr(y));
0068     x += dev_Rp * +y / R + dev_R * -x / R;
0069     y += dev_Rp * -x / R + dev_R * -y / R;
0070     z += dist_z(generator);
0071   } else {
0072     x += dist_xyh(generator);
0073     y += dist_xyh(generator);
0074     z += dist_zh(generator);
0075   }
0076 }
0077
0078 template <int N>
0079 void Hits_cov(Eigen::Matrix<float, 6, 4>& V,
0080               const unsigned int& i,
0081               const unsigned int& n,
0082               const Matrix3xNd<N>& hits,
0083               const Vector5d& err,
0084               bool isbarrel) {
0085   if (isbarrel) {
0086     double R2 = riemannFit::sqr(hits(0, i)) + riemannFit::sqr(hits(1, i));
0087     V.col(i)[0] = (riemannFit::sqr(err[1]) * riemannFit::sqr(hits(1, i)) +
0088                    riemannFit::sqr(err[0]) * riemannFit::sqr(hits(0, i))) /
0089                   R2;
0090     V.col(i)[2] = (riemannFit::sqr(err[1]) * riemannFit::sqr(hits(0, i)) +
0091                    riemannFit::sqr(err[0]) * riemannFit::sqr(hits(1, i))) /
0092                   R2;
0093     V.col(i)[1] = (riemannFit::sqr(err[0]) - riemannFit::sqr(err[1])) * hits(1, i) * hits(0, i) / R2;
0094     V.col(i)[5] = riemannFit::sqr(err[2]);
0095   } else {
0096     V.col(i)[0] = riemannFit::sqr(err[3]);
0097     V.col(i)[2] = riemannFit::sqr(err[3]);
0098     V.col(i)[5] = riemannFit::sqr(err[4]);
0099   }
0100 }
0101
0102 hits_gen Hits_gen(const unsigned int& n, const Matrix<double, 6, 1>& gen_par) {
0103   hits_gen gen;
0104   gen.hits = MatrixXd::Zero(3, n);
0105   gen.hits_ge = Eigen::Matrix<float, 6, 4>::Zero();
0106   // err /= 10000.;
0107   constexpr double rad[8] = {2.95, 6.8, 10.9, 16., 3.1, 7., 11., 16.2};
0108   // constexpr double R_err[8] = {5./10000, 5./10000, 5./10000, 5./10000, 5./10000,
0109   // 5./10000, 5./10000, 5./10000};  constexpr double Rp_err[8] = {35./10000, 18./10000,
0110   // 15./10000, 34./10000, 35./10000, 18./10000, 15./10000, 34./10000};  constexpr double z_err[8] =
0111   // {72./10000, 38./10000, 25./10000, 56./10000, 72./10000, 38./10000, 25./10000, 56./10000};
0112   constexpr double R_err[8] = {
0113       10. / 10000, 10. / 10000, 10. / 10000, 10. / 10000, 10. / 10000, 10. / 10000, 10. / 10000, 10. / 10000};
0114   constexpr double Rp_err[8] = {
0115       35. / 10000, 18. / 10000, 15. / 10000, 34. / 10000, 35. / 10000, 18. / 10000, 15. / 10000, 34. / 10000};
0116   constexpr double z_err[8] = {
0117       72. / 10000, 38. / 10000, 25. / 10000, 56. / 10000, 72. / 10000, 38. / 10000, 25. / 10000, 56. / 10000};
0118   const double x2 = gen_par(0) + gen_par(4) * cos(gen_par(3) * pi / 180);
0119   const double y2 = gen_par(1) + gen_par(4) * sin(gen_par(3) * pi / 180);
0120   const double alpha = atan2(y2, x2);
0121
0122   for (unsigned int i = 0; i < n; ++i) {
0123     const double a = gen_par(4);
0124     const double b = rad[i];
0125     const double c = sqrt(riemannFit::sqr(x2) + riemannFit::sqr(y2));
0126     const double beta = acos((riemannFit::sqr(a) - riemannFit::sqr(b) - riemannFit::sqr(c)) / (-2. * b * c));
0127     const double gamma = alpha + beta;
0128     gen.hits(0, i) = rad[i] * cos(gamma);
0129     gen.hits(1, i) = rad[i] * sin(gamma);
0130     gen.hits(2, i) =
0131         gen_par(2) +
0132         1 / tan(gen_par(5) * pi / 180) * 2. *
0133             asin(sqrt(riemannFit::sqr((gen_par(0) - gen.hits(0, i))) + riemannFit::sqr((gen_par(1) - gen.hits(1, i)))) /
0134                  (2. * gen_par(4))) *
0135             gen_par(4);
0136     // isbarrel(i) = ??
0137     Vector5d err;
0138     err << R_err[i], Rp_err[i], z_err[i], 0, 0;
0139     smearing(err, true, gen.hits(0, i), gen.hits(1, i), gen.hits(2, i));
0140     Hits_cov(gen.hits_ge, i, n, gen.hits, err, true);
0141   }
0142
0143   return gen;
0144 }
0145
0146 Vector5d True_par(const Matrix<double, 6, 1>& gen_par, const int& charge, const double& B_field) {
0147   Vector5d true_par;
0148   const double x0 = gen_par(0) + gen_par(4) * cos(gen_par(3) * pi / 180);
0149   const double y0 = gen_par(1) + gen_par(4) * sin(gen_par(3) * pi / 180);
0150   CircleFit circle;
0151   circle.par << x0, y0, gen_par(4);
0152   circle.qCharge = 1;
0153   riemannFit::par_uvrtopak(circle, B_field, false);
0154   true_par.block(0, 0, 3, 1) = circle.par;
0155   true_par(3) = 1 / tan(gen_par(5) * pi / 180);
0156   const int dir = ((gen_par(0) - cos(true_par(0) - pi / 2) * true_par(1)) * (gen_par(1) - y0) -
0157                        (gen_par(1) - sin(true_par(0) - pi / 2) * true_par(1)) * (gen_par(0) - x0) >
0158                    0)
0159                       ? -1
0160                       : 1;
0161   true_par(4) = gen_par(2) + 1 / tan(gen_par(5) * pi / 180) * dir * 2.f *
0162                                  asin(sqrt(riemannFit::sqr((gen_par(0) - cos(true_par(0) - pi / 2) * true_par(1))) +
0163                                            riemannFit::sqr((gen_par(1) - sin(true_par(0) - pi / 2) * true_par(1)))) /
0164                                       (2.f * gen_par(4))) *
0165                                  gen_par(4);
0166   return true_par;
0167 }
0168
0169 Matrix<double, 6, 1> New_par(const Matrix<double, 6, 1>& gen_par, const int& charge, const double& B_field) {
0170   Matrix<double, 6, 1> new_par;
0171   new_par.block(0, 0, 3, 1) = gen_par.block(0, 0, 3, 1);
0172   new_par(3) = gen_par(3) - charge * 90;
0173   new_par(4) = gen_par(4) / B_field;
0174   //  new_par(5) = atan(sinh(gen_par(5))) * 180 / pi;
0175   new_par(5) = 2. * atan(exp(-gen_par(5))) * 180 / pi;
0176   return new_par;
0177 }
0178
0179 template <typename Fit, size_t N>
0180 void computePull(std::array<Fit, N>& fit, const char* label, int n_, int iteration, const Vector5d& true_par) {
0181   Eigen::Matrix<double, 41, Eigen::Dynamic, 1> score(41, iteration);
0182
0183   std::string histo_name("Phi Pull");
0184   histo_name += label;
0185   TH1F phi_pull(histo_name.data(), histo_name.data(), 100, -10., 10.);
0186   histo_name = "dxy Pull ";
0187   histo_name += label;
0188   TH1F dxy_pull(histo_name.data(), histo_name.data(), 100, -10., 10.);
0189   histo_name = "dz Pull ";
0190   histo_name += label;
0191   TH1F dz_pull(histo_name.data(), histo_name.data(), 100, -10., 10.);
0192   histo_name = "Theta Pull ";
0193   histo_name += label;
0194   TH1F theta_pull(histo_name.data(), histo_name.data(), 100, -10., 10.);
0195   histo_name = "Pt Pull ";
0196   histo_name += label;
0197   TH1F pt_pull(histo_name.data(), histo_name.data(), 100, -10., 10.);
0198   histo_name = "Phi Error ";
0199   histo_name += label;
0200   TH1F phi_error(histo_name.data(), histo_name.data(), 100, 0., 0.1);
0201   histo_name = "dxy error ";
0202   histo_name += label;
0203   TH1F dxy_error(histo_name.data(), histo_name.data(), 100, 0., 0.1);
0204   histo_name = "dz error ";
0205   histo_name += label;
0206   TH1F dz_error(histo_name.data(), histo_name.data(), 100, 0., 0.1);
0207   histo_name = "Theta error ";
0208   histo_name += label;
0209   TH1F theta_error(histo_name.data(), histo_name.data(), 100, 0., 0.1);
0210   histo_name = "Pt error ";
0211   histo_name += label;
0212   TH1F pt_error(histo_name.data(), histo_name.data(), 100, 0., 0.1);
0213   for (int x = 0; x < iteration; x++) {
0214     // Compute PULLS information
0215     score(0, x) = (fit[x].par(0) - true_par(0)) / sqrt(fit[x].cov(0, 0));
0216     score(1, x) = (fit[x].par(1) - true_par(1)) / sqrt(fit[x].cov(1, 1));
0217     score(2, x) = (fit[x].par(2) - true_par(2)) / sqrt(fit[x].cov(2, 2));
0218     score(3, x) = (fit[x].par(3) - true_par(3)) / sqrt(fit[x].cov(3, 3));
0219     score(4, x) = (fit[x].par(4) - true_par(4)) / sqrt(fit[x].cov(4, 4));
0220     phi_pull.Fill(score(0, x));
0221     dxy_pull.Fill(score(1, x));
0222     pt_pull.Fill(score(2, x));
0223     theta_pull.Fill(score(3, x));
0224     dz_pull.Fill(score(4, x));
0225     phi_error.Fill(sqrt(fit[x].cov(0, 0)));
0226     dxy_error.Fill(sqrt(fit[x].cov(1, 1)));
0227     pt_error.Fill(sqrt(fit[x].cov(2, 2)));
0228     theta_error.Fill(sqrt(fit[x].cov(3, 3)));
0229     dz_error.Fill(sqrt(fit[x].cov(4, 4)));
0230     score(5, x) = (fit[x].par(0) - true_par(0)) * (fit[x].par(1) - true_par(1)) / (fit[x].cov(0, 1));
0231     score(6, x) = (fit[x].par(0) - true_par(0)) * (fit[x].par(2) - true_par(2)) / (fit[x].cov(0, 2));
0232     score(7, x) = (fit[x].par(1) - true_par(1)) * (fit[x].par(2) - true_par(2)) / (fit[x].cov(1, 2));
0233     score(8, x) = (fit[x].par(3) - true_par(3)) * (fit[x].par(4) - true_par(4)) / (fit[x].cov(3, 4));
0234     score(9, x) = fit[x].chi2_circle;
0235     score(25, x) = fit[x].chi2_line;
0236     score(10, x) = sqrt(fit[x].cov(0, 0)) / fit[x].par(0) * 100;
0237     score(13, x) = sqrt(fit[x].cov(3, 3)) / fit[x].par(3) * 100;
0238     score(14, x) = sqrt(fit[x].cov(4, 4)) / fit[x].par(4) * 100;
0239     score(15, x) =
0240         (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));
0241     score(16, x) =
0242         (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));
0243     score(17, x) =
0244         (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));
0245     score(18, x) =
0246         (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));
0247     score(19, x) =
0248         (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));
0249     score(20, x) =
0250         (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));
0251     score(21, x) =
0252         (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));
0253     score(22, x) =
0254         (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));
0255     score(23, x) =
0256         (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));
0257     score(24, x) =
0258         (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));
0259     score(30, x) = fit[x].par(0);
0260     score(31, x) = fit[x].par(1);
0261     score(32, x) = fit[x].par(2);
0262     score(33, x) = fit[x].par(3);
0263     score(34, x) = fit[x].par(4);
0264     score(35, x) = sqrt(fit[x].cov(0, 0));
0265     score(36, x) = sqrt(fit[x].cov(1, 1));
0266     score(37, x) = sqrt(fit[x].cov(2, 2));
0267     score(38, x) = sqrt(fit[x].cov(3, 3));
0268     score(39, x) = sqrt(fit[x].cov(4, 4));
0269   }
0270
0271   double phi_ = score.row(0).mean();
0272   double a_ = score.row(1).mean();
0273   double pt_ = score.row(2).mean();
0274   double coT_ = score.row(3).mean();
0275   double Zip_ = score.row(4).mean();
0276   std::cout << std::setprecision(5) << std::scientific << label << " AVERAGE FITTED VALUES: \n"
0277             << "phi: " << score.row(30).mean() << " +/- " << score.row(35).mean() << " [+/-] "
0278             << sqrt(score.row(35).array().abs2().mean() - score.row(35).mean() * score.row(35).mean()) << std::endl
0279             << "d0:  " << score.row(31).mean() << " +/- " << score.row(36).mean() << " [+/-] "
0280             << sqrt(score.row(36).array().abs2().mean() - score.row(36).mean() * score.row(36).mean()) << std::endl
0281             << "pt:  " << score.row(32).mean() << " +/- " << score.row(37).mean() << " [+/-] "
0282             << sqrt(score.row(37).array().abs2().mean() - score.row(37).mean() * score.row(37).mean()) << std::endl
0283             << "coT: " << score.row(33).mean() << " +/- " << score.row(38).mean() << " [+/-] "
0284             << sqrt(score.row(38).array().abs2().mean() - score.row(38).mean() * score.row(38).mean()) << std::endl
0285             << "Zip: " << score.row(34).mean() << " +/- " << score.row(39).mean() << " [+/-] "
0286             << sqrt(score.row(39).array().abs2().mean() - score.row(39).mean() * score.row(39).mean()) << std::endl;
0287
0288   Matrix5d correlation;
0289   correlation << 1., score.row(21).mean(), score.row(22).mean(), score.row(15).mean(), score.row(20).mean(),
0290       score.row(21).mean(), 1., score.row(23).mean(), score.row(16).mean(), score.row(19).mean(), score.row(22).mean(),
0291       score.row(23).mean(), 1., score.row(17).mean(), score.row(20).mean(), score.row(15).mean(), score.row(16).mean(),
0292       score.row(17).mean(), 1., score.row(24).mean(), score.row(18).mean(), score.row(19).mean(), score.row(20).mean(),
0293       score.row(24).mean(), 1.;
0294
0295   cout << "\n"
0296        << label << " PULLS (mean, sigma, relative_error):\n"
0297        << "phi:  " << phi_ << "     " << sqrt((score.row(0).array() - phi_).square().sum() / (iteration - 1)) << "   "
0298        << abs(score.row(10).mean()) << "%\n"
0299        << "a0 :  " << a_ << "     " << sqrt((score.row(1).array() - a_).square().sum() / (iteration - 1)) << "   "
0300        << abs(score.row(11).mean()) << "%\n"
0301        << "pt :  " << pt_ << "     " << sqrt((score.row(2).array() - pt_).square().sum() / (iteration - 1)) << "   "
0302        << abs(score.row(12).mean()) << "%\n"
0303        << "coT:  " << coT_ << "     " << sqrt((score.row(3).array() - coT_).square().sum() / (iteration - 1)) << "   "
0304        << abs(score.row(13).mean()) << "%\n"
0305        << "Zip:  " << Zip_ << "     " << sqrt((score.row(4).array() - Zip_).square().sum() / (iteration - 1)) << "   "
0306        << abs(score.row(14).mean()) << "%\n\n"
0307        << "cov(phi,a0)_:  " << score.row(5).mean() << "\n"
0308        << "cov(phi,pt)_:  " << score.row(6).mean() << "\n"
0309        << "cov(a0,pt)_:   " << score.row(7).mean() << "\n"
0310        << "cov(coT,Zip)_: " << score.row(8).mean() << "\n\n"
0311        << "chi2_circle:  " << score.row(9).mean() << " vs " << n_ - 3 << "\n"
0312        << "chi2_line:    " << score.row(25).mean() << " vs " << n_ - 2 << "\n\n"
0313        << "correlation matrix:\n"
0314        << correlation << "\n\n"
0315        << endl;
0316
0317   phi_pull.Fit("gaus", "Q");
0318   dxy_pull.Fit("gaus", "Q");
0319   dz_pull.Fit("gaus", "Q");
0320   theta_pull.Fit("gaus", "Q");
0321   pt_pull.Fit("gaus", "Q");
0322   phi_pull.Write();
0323   dxy_pull.Write();
0324   dz_pull.Write();
0325   theta_pull.Write();
0326   pt_pull.Write();
0327   phi_error.Write();
0328   dxy_error.Write();
0329   dz_error.Write();
0330   theta_error.Write();
0331   pt_error.Write();
0332 }
0333
0334 void test_helix_fit(bool getcin) {
0335   int n_;
0336   const double B_field = 3.8 * c_speed / pow(10, 9) / 100;
0337   Matrix<double, 6, 1> gen_par;
0338   Vector5d true_par;
0339   Vector5d err;
0340   generator.seed(1);
0341   std::cout << std::setprecision(6);
0342   cout << "_________________________________________________________________________\n";
0343   cout << "n x(cm) y(cm) z(cm) phi(grad) R(Gev/c) eta iteration debug" << endl;
0344   if (getcin) {
0345     cout << "hits: ";
0346     cin >> n_;
0347     cout << "x: ";
0348     cin >> gen_par(0);
0349     cout << "y: ";
0350     cin >> gen_par(1);
0351     cout << "z: ";
0352     cin >> gen_par(2);
0353     cout << "phi: ";
0354     cin >> gen_par(3);
0355     cout << "p_t: ";
0356     cin >> gen_par(4);
0357     cout << "eta: ";
0358     cin >> gen_par(5);
0359   } else {
0360     n_ = 4;
0361     gen_par(0) = -0.1;  // x
0362     gen_par(1) = 0.1;   // y
0363     gen_par(2) = -1.;   // z
0364     gen_par(3) = 45.;   // phi
0365     gen_par(4) = 10.;   // R (p_t)
0366     gen_par(5) = 1.;    // eta
0367   }
0368
0369   const int iteration = 5000;
0370   gen_par = New_par(gen_par, 1, B_field);
0371   true_par = True_par(gen_par, 1, B_field);
0372   std::array<HelixFit, iteration> helixRiemann_fit;
0373
0374   std::cout << "\nTrue parameters: "
0375             << "phi: " << true_par(0) << " "
0376             << "dxy: " << true_par(1) << " "
0377             << "pt: " << true_par(2) << " "
0378             << "CotT: " << true_par(3) << " "
0379             << "Zip: " << true_par(4) << " " << std::endl;
0380   auto start = std::chrono::high_resolution_clock::now();
0381   auto delta = start - start;
0382   for (int i = 0; i < 100 * iteration; i++) {
0383     hits_gen gen;
0384     gen = Hits_gen(n_, gen_par);
0385     //      gen.hits = MatrixXd::Zero(3, 4);
0386     //      gen.hits_cov = MatrixXd::Zero(3 * 4, 3 * 4);
0387     //      gen.hits.col(0) << 1.82917642593, 2.0411875248, 7.18495464325;
0388     //      gen.hits.col(1) << 4.47041416168, 4.82704305649, 18.6394691467;
0389     //      gen.hits.col(2) << 7.25991010666, 7.74653434753, 30.6931324005;
0390     //      gen.hits.col(3) << 8.99161434174, 9.54262828827, 38.1338043213;
0391     delta -= std::chrono::high_resolution_clock::now() - start;
0392     helixRiemann_fit[i % iteration] =
0393 #ifdef USE_BL
0394         brokenline::helixFit(gen.hits, gen.hits_ge, B_field);
0395 #else
0396         riemannFit::helixFit(gen.hits, gen.hits_ge, B_field, true);
0397 #endif
0398     delta += std::chrono::high_resolution_clock::now() - start;
0399
0400     if (helixRiemann_fit[i % iteration].par(0) > 10.)
0401       std::cout << "error" << std::endl;
0402     if (0 == i)
0403       cout << std::setprecision(6) << "phi:  " << helixRiemann_fit[i].par(0) << " +/- "
0404            << sqrt(helixRiemann_fit[i].cov(0, 0)) << " vs " << true_par(0) << endl
0405            << "Tip:  " << helixRiemann_fit[i].par(1) << " +/- " << sqrt(helixRiemann_fit[i].cov(1, 1)) << " vs "
0406            << true_par(1) << endl
0407            << "p_t:  " << helixRiemann_fit[i].par(2) << " +/- " << sqrt(helixRiemann_fit[i].cov(2, 2)) << " vs "
0408            << true_par(2) << endl
0409            << "theta:" << helixRiemann_fit[i].par(3) << " +/- " << sqrt(helixRiemann_fit[i].cov(3, 3)) << " vs "
0410            << true_par(3) << endl
0411            << "Zip:  " << helixRiemann_fit[i].par(4) << " +/- " << sqrt(helixRiemann_fit[i].cov(4, 4)) << " vs "
0412            << true_par(4) << endl
0413            << "charge:" << helixRiemann_fit[i].qCharge << " vs 1" << endl
0414            << "covariance matrix:" << endl
0415            << helixRiemann_fit[i].cov << endl
0416            << "Initial hits:\n"
0417            << gen.hits << endl
0418            << "Initial Covariance:\n"
0419            << gen.hits_ge << endl;
0420   }
0421   std::cout << "elapsted time " << double(std::chrono::duration_cast<std::chrono::nanoseconds>(delta).count()) / 1.e6
0422             << std::endl;
0423   computePull(helixRiemann_fit, "Riemann", n_, iteration, true_par);
0424 }
0425
0426 int main(int nargs, char**) {
0427   TFile f("TestFitResults.root", "RECREATE");
0428   test_helix_fit(nargs > 1);
0429   f.Close();
0430   return 0;
0431 }