File indexing completed on 2023-03-31 23:02:24
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 };
0030
0031
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
0107 constexpr double rad[8] = {2.95, 6.8, 10.9, 16., 3.1, 7., 11., 16.2};
0108
0109
0110
0111
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
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
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
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;
0362 gen_par(1) = 0.1;
0363 gen_par(2) = -1.;
0364 gen_par(3) = 45.;
0365 gen_par(4) = 10.;
0366 gen_par(5) = 1.;
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
0386
0387
0388
0389
0390
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 }