Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-05-27 07:24:18 

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008 
0009 // ROOT include(s).
0010 #include <Math/ProbFuncMathCore.h>
0011 #include <TCanvas.h>
0012 #include <TF1.h>
0013 #include <TFile.h>
0014 #include <TFitResult.h>
0015 #include <TH1D.h>
0016 #include <TLatex.h>
0017 #include <TLine.h>
0018 #include <TMath.h>
0019 #include <TROOT.h>
0020 #include <TStyle.h>
0021 #include <TTree.h>
0022 
0023 #include <ROOT/RCsvDS.hxx>
0024 #include <ROOT/RDataFrame.hxx>
0025 
0026 // System include(s).
0027 #include <fstream>
0028 #include <iomanip>
0029 #include <iostream>
0030 #include <sstream>
0031 #include <vector>
0032 
0033 namespace {
0034 const double x_pos = 0.205;
0035 const double title_x = x_pos;
0036 const double title_y = 0.8197;
0037 const double y_gap = -0.0505;
0038 const double header_text_size = 0.055;
0039 const double geom_text_size = 0.0434028;
0040 
0041 const double pull_fit_title_x = x_pos;
0042 const double pull_fit_title_y = 0.700;
0043 const double pval_fit_title_x = x_pos;
0044 const double pval_fit_title_y = 0.700;
0045 const double gaus_fit_par_x = x_pos;
0046 const double number_offset = 0.125;
0047 const double gaus_fit_par_y = pull_fit_title_y - 0.065;
0048 const double const_fit_par_x = x_pos;
0049 const double const_fit_par_y = pval_fit_title_y - 0.0459;
0050 const double tolerance_x = 0.7;
0051 const double tolerance_y = 0.67;
0052 const double pull_text_size = 0.0434028;
0053 const double pval_text_size = 0.0434028;
0054 const double pad_x0 = 0.00;
0055 const double pad_x1 = 1.;
0056 const double pad_y0 = 0.00;
0057 const double pad_y1 = 1.;
0058 const int label_font = 132;
0059 const double label_font_size = 0.055;
0060 const double titleX_font_size = 0.055;
0061 const double titleY_font_size = 0.055;
0062 const double x_title_offset = 1.25;
0063 const double y_title_offset = 1.34;
0064 const double y_title_offset_pval = 0.9;
0065 const double x_label_offset = 0.015;
0066 const double y_label_offset = 0.015;
0067 const double pull_min = -6.5;
0068 const double pull_max = 6.5;
0069 
0070 }  // namespace
0071 
0072 auto get_tree(std::string name) {
0073 
0074     const std::string csv_name = name + ".csv";
0075     const std::string root_name = name + ".root";
0076 
0077     auto rdf = ROOT::RDF::FromCSV(csv_name);
0078 
0079     // Create root file
0080     rdf.Snapshot(name, root_name);
0081 
0082     auto f = TFile::Open(root_name.c_str(), "update");
0083     auto t = (TTree*)f->Get(name.c_str());
0084 
0085     return t;
0086 }
0087 
0088 void draw_text(double x1, double y1, double y_delta, double s1, double s2,
0089                std::string t1, std::string t2) {
0090     TLatex* ttext1 = new TLatex(x1, y1, t1.c_str());
0091     TLatex* ttext2 = new TLatex(x1, y1 + y_delta, t2.c_str());
0092     ttext1->SetTextFont(22);
0093     ttext1->SetTextSize(s1);
0094     ttext2->SetTextFont(132);
0095     ttext2->SetTextSize(s2);
0096 
0097     ttext1->Draw();
0098     ttext2->Draw();
0099 }
0100 
0101 std::pair<std::array<double, 3u>, std::array<double, 3u>> fit_pull(
0102     TH1D* h_pull, std::array<double, 14u>& arr_pull) {
0103 
0104     // Function used for the fit.
0105     TF1 gaus{"gaus", "gaus", pull_min, pull_max};
0106     double fit_par[3];
0107     double fit_par_error[3];
0108 
0109     // Set the mean seed to 0
0110     gaus.SetParameters(1, 0);
0111     gaus.SetParLimits(1, -1., 1.);
0112     // Set the standard deviation seed to 1
0113     gaus.SetParameters(2, 1.0);
0114     gaus.SetParLimits(2, 0.5, 2.);
0115 
0116     auto res = h_pull->Fit("gaus", "Q0S");
0117     gaus.GetParameters(&fit_par[0]);
0118 
0119     std::array<double, 3u> par{fit_par[0], fit_par[1], fit_par[2]};
0120     std::array<double, 3u> error;
0121     error[0] = gaus.GetParError(0);
0122     error[1] = gaus.GetParError(1);
0123     error[2] = gaus.GetParError(2);
0124 
0125     arr_pull[0u] = fit_par[0u];
0126     arr_pull[1u] = fit_par[1u];
0127     arr_pull[2u] = fit_par[2u];
0128     arr_pull[3u] = error[0u];
0129     arr_pull[4u] = error[1u];
0130     arr_pull[5u] = error[2u];
0131     arr_pull[6u] = res->Ndf();
0132     arr_pull[7u] = res->Chi2();
0133     arr_pull[8u] = arr_pull[6u] / arr_pull[7u];
0134     arr_pull[9u] = ROOT::Math::chisquared_cdf_c(arr_pull[7u], arr_pull[6u]);
0135 
0136     return {par, error};
0137 }
0138 
0139 std::pair<double, double> fit_pval(TH1D* h_pval) {
0140 
0141     // Function used for the fit.
0142     TF1 unif{"uniform", "[0]", 0.f, 1.f};
0143     double fit_par[1];
0144 
0145     auto res = h_pval->Fit("uniform", "Q0S");
0146     unif.GetParameters(&fit_par[0]);
0147     double error = unif.GetParError(0);
0148 
0149     return {fit_par[0], error};
0150 }
0151 
0152 void set_yaxis_title(TH1D* h, const double text_size) {
0153     double bin_width = h->GetBinWidth(0u);
0154     std::string str = std::to_string(bin_width);
0155     str.erase(str.find_last_not_of('0') + 1, std::string::npos);
0156     str.erase(str.find_last_not_of('.') + 1, std::string::npos);
0157     std::string y_axis_title = "Counts / (" + str + ")";
0158     h->GetYaxis()->SetTitle(y_axis_title.c_str());
0159     h->GetYaxis()->SetTitleSize(text_size);
0160 }
0161 
0162 void set_xaxis_title(TH1D* h, const double text_size) {
0163 
0164     std::string x_axis_title;
0165 
0166     const TString h_name = h->GetName();
0167 
0168     if (h_name.Contains("l0")) {
0169         x_axis_title = "PL(l_{0F})";
0170     } else if (h_name.Contains("l1")) {
0171         x_axis_title = "PL(l_{1F})";
0172     } else if (h_name.Contains("phi")) {
0173         x_axis_title = "PL(#phi_{F})";
0174     } else if (h_name.Contains("theta")) {
0175         x_axis_title = "PL(#theta_{F})";
0176     } else if (h_name.Contains("qop")) {
0177         x_axis_title = "PL(#lambda_{F})";
0178     } else if (h_name.Contains("pval")) {
0179         x_axis_title = "p-value";
0180     }
0181 
0182     h->GetXaxis()->SetTitle(x_axis_title.c_str());
0183     h->GetXaxis()->SetTitleSize(text_size);
0184 }
0185 
0186 void draw_fit_title(const std::string title, const double x, const double y,
0187                     const double text_size) {
0188 
0189     TLatex* ttext = new TLatex(x, y, title.c_str());
0190     ttext->SetTextFont(22);
0191     ttext->SetTextSize(text_size);
0192     ttext->Draw();
0193     gPad->cd();
0194 }
0195 
0196 void draw_gaus_fit_par(const std::array<double, 3u>& fit_par,
0197                        const std::array<double, 3u>& fit_par_error,
0198                        const double x, const double y, const double text_size) {
0199     TLatex* ttext = new TLatex(x, y, "#splitline{Mean}{Std Dev}");
0200     ttext->SetTextFont(132);
0201     ttext->SetTextSize(text_size);
0202     ttext->Draw();
0203 
0204     std::stringstream mean_stream;
0205     mean_stream << std::fixed << std::setprecision(3) << fit_par[1] << " #pm "
0206                 << fit_par_error[1];
0207     std::stringstream sigma_stream;
0208     sigma_stream << std::fixed << std::setprecision(3) << fit_par[2] << " #pm "
0209                  << fit_par_error[2];
0210 
0211     TLatex* ttext2 = new TLatex(x + number_offset, y,
0212                                 "#splitline{" + TString(mean_stream.str()) +
0213                                     "}{" + TString(sigma_stream.str()) + "}");
0214     ttext2->SetTextFont(132);
0215     ttext2->SetTextSize(text_size);
0216     ttext2->Draw();
0217 }
0218 
0219 void draw_const_fit_par(const double fit_par, const double fit_par_error,
0220                         const double x, const double y,
0221                         const double text_size) {
0222     std::stringstream val_stream;
0223     val_stream << std::fixed << std::setprecision(3) << fit_par << " #pm "
0224                << fit_par_error;
0225 
0226     TLatex* ttext = new TLatex(x, y, "Value  " + TString(val_stream.str()));
0227     ttext->SetTextFont(132);
0228     ttext->SetTextSize(text_size);
0229     ttext->Draw();
0230 }
0231 
0232 void draw_tolerance(const double log10_rk_tolerance, const double x,
0233                     const double y) {
0234     std::stringstream val_stream;
0235     val_stream << "#tau = 10^{" << int(log10_rk_tolerance) << "} mm";
0236 
0237     TLatex* ttext = new TLatex(x, y, TString(val_stream.str()));
0238     ttext->SetTextFont(132);
0239     ttext->SetTextSize(0.052);
0240     ttext->Draw();
0241 }
0242 
0243 void draw_pull(TH1D* h_pull, const std::string& header_text,
0244                const std::string& geom_text, const double log10_rk_tol,
0245                std::array<double, 14u>& arr_pull) {
0246 
0247     TPad* pull_pad =
0248         new TPad("pull_pad", "pull_pad", pad_x0, pad_y0, pad_x1, pad_y1);
0249     pull_pad->SetLogy();
0250     pull_pad->Draw();
0251     pull_pad->cd();
0252 
0253     pull_pad->SetLeftMargin(110. / pull_pad->GetWw());
0254     pull_pad->SetBottomMargin(95. / pull_pad->GetWh());
0255 
0256     auto fit_res = fit_pull(h_pull, arr_pull);
0257     auto fit_pars = fit_res.first;
0258     auto fit_errors = fit_res.second;
0259 
0260     set_xaxis_title(h_pull, pull_text_size);
0261     set_yaxis_title(h_pull, pull_text_size);
0262     const double y_axis_max = h_pull->GetEntries() * 50.;
0263     h_pull->GetYaxis()->SetRangeUser(0.5f, y_axis_max);
0264     h_pull->GetXaxis()->SetLabelFont(label_font);
0265     h_pull->GetYaxis()->SetLabelFont(label_font);
0266     h_pull->GetXaxis()->SetLabelSize(label_font_size);
0267     h_pull->GetYaxis()->SetLabelSize(label_font_size);
0268     h_pull->GetXaxis()->SetTitleSize(titleX_font_size);
0269     h_pull->GetYaxis()->SetTitleSize(titleY_font_size);
0270     h_pull->GetYaxis()->SetTitleOffset(y_title_offset);
0271     h_pull->GetXaxis()->SetTitleOffset(x_title_offset + 0.1);
0272     h_pull->GetXaxis()->SetLabelOffset(x_label_offset);
0273     h_pull->GetYaxis()->SetLabelOffset(y_label_offset);
0274     h_pull->GetYaxis()->SetNdivisions(504);
0275     h_pull->GetYaxis()->SetMaxDigits(1);
0276     h_pull->GetXaxis()->CenterTitle(true);
0277     h_pull->GetYaxis()->CenterTitle(true);
0278     h_pull->GetXaxis()->SetTitleFont(132);
0279     h_pull->GetYaxis()->SetTitleFont(132);
0280     h_pull->Draw();
0281     TF1* gaus = new TF1(h_pull->GetName(), "gaus", pull_min, pull_max);
0282     gaus->SetParameters(fit_pars[0], fit_pars[1], fit_pars[2]);
0283     gaus->Draw("same");
0284 
0285     TPad* text_pad = new TPad("gaus_text_pad", "gaus_text_pad", pad_x0, pad_y0,
0286                               pad_x1, pad_y1);
0287     text_pad->SetFillStyle(4000);
0288     text_pad->Draw();
0289     text_pad->cd();
0290 
0291     draw_text(title_x, title_y, y_gap, header_text_size, geom_text_size,
0292               header_text.c_str(), geom_text.c_str());
0293     draw_fit_title("Gaussian fit", pull_fit_title_x, pull_fit_title_y,
0294                    pull_text_size);
0295     draw_gaus_fit_par(fit_pars, fit_errors, gaus_fit_par_x, gaus_fit_par_y,
0296                       pull_text_size);
0297     // draw_tolerance(log10_rk_tol, tolerance_x, tolerance_y);
0298 }
0299 
0300 void draw_pval(TH1D* h_pval, const std::string& header_text,
0301                const std::string& geom_text, const double log10_rk_tol) {
0302 
0303     TPad* pval_pad =
0304         new TPad("pval_pad", "pval_pad", pad_x0, pad_y0, pad_x1, pad_y1);
0305     pval_pad->Draw();
0306     pval_pad->cd();
0307 
0308     pval_pad->SetLeftMargin(110. / pval_pad->GetWw());
0309     pval_pad->SetBottomMargin(95. / pval_pad->GetWh());
0310 
0311     auto fit_res = fit_pval(h_pval);
0312     auto fit_par = fit_res.first;
0313     auto fit_error = fit_res.second;
0314     set_xaxis_title(h_pval, pval_text_size);
0315     set_yaxis_title(h_pval, pval_text_size);
0316     const double y_axis_max = 2. * h_pval->GetEntries() / h_pval->GetNbinsX();
0317     h_pval->GetYaxis()->SetRangeUser(0.f, y_axis_max);
0318     h_pval->GetXaxis()->SetLabelFont(label_font);
0319     h_pval->GetYaxis()->SetLabelFont(label_font);
0320     h_pval->GetXaxis()->SetLabelSize(label_font_size);
0321     h_pval->GetYaxis()->SetLabelSize(label_font_size);
0322     h_pval->GetXaxis()->SetTitleSize(titleX_font_size);
0323     h_pval->GetYaxis()->SetTitleSize(titleY_font_size);
0324     h_pval->GetYaxis()->SetTitleOffset(y_title_offset_pval);
0325     h_pval->GetXaxis()->SetTitleOffset(x_title_offset);
0326     h_pval->GetYaxis()->SetLabelOffset(y_label_offset);
0327     h_pval->GetXaxis()->SetLabelOffset(x_label_offset);
0328     h_pval->GetXaxis()->SetNdivisions(505);
0329     h_pval->GetYaxis()->SetMaxDigits(2);
0330     h_pval->GetYaxis()->SetNdivisions(505);
0331     h_pval->GetYaxis()->SetDecimals();
0332     h_pval->GetXaxis()->CenterTitle(true);
0333     h_pval->GetYaxis()->CenterTitle(true);
0334     h_pval->GetXaxis()->SetTitleFont(132);
0335     h_pval->GetYaxis()->SetTitleFont(132);
0336 
0337     h_pval->Draw();
0338     TLine* line = new TLine(0.f, fit_par, 1.f, fit_par);
0339     line->SetLineColor(kRed);
0340     line->SetLineWidth(2);
0341     line->Draw();
0342 
0343     TPad* text_pad = new TPad("const_text_pad", "const_text_pad", 0, 0, 1, 1);
0344     text_pad->SetFillStyle(4000);
0345     text_pad->Draw();
0346     text_pad->cd();
0347 
0348     draw_text(title_x, title_y, y_gap, header_text_size, geom_text_size,
0349               header_text.c_str(), geom_text.c_str());
0350     draw_fit_title("Constant fit", pval_fit_title_x, pval_fit_title_y,
0351                    pval_text_size);
0352     draw_const_fit_par(fit_par, fit_error, const_fit_par_x, const_fit_par_y,
0353                        pval_text_size);
0354     // draw_tolerance(log10_rk_tol, tolerance_x, tolerance_y);
0355 }
0356 
0357 std::string to_pdf(const std::string& name) {
0358     return name + ".pdf";
0359 }
0360 
0361 void read_tree(TTree* t, const std::string& tag,
0362                const std::string& header_title, const std::string& geom_title) {
0363     const std::array<float, 2> cdim1{700, 600};
0364     const std::array<float, 2> cdim2{700, 600};
0365 
0366     int n_bins = (pull_max - pull_min) / 0.1;
0367 
0368     double pull_l0;
0369     double pull_l1;
0370     double pull_phi;
0371     double pull_theta;
0372     double pull_qop;
0373     double chi2;
0374     double log10_rk_tolerance;
0375     double p_value;
0376 
0377     t->SetBranchAddress("pull_l0", &pull_l0);
0378     t->SetBranchAddress("pull_l1", &pull_l1);
0379     t->SetBranchAddress("pull_phi", &pull_phi);
0380     t->SetBranchAddress("pull_theta", &pull_theta);
0381     t->SetBranchAddress("pull_qop", &pull_qop);
0382     t->SetBranchAddress("chi2", &chi2);
0383     t->SetBranchAddress("log10_rk_tolerance", &log10_rk_tolerance);
0384     auto b_pval = t->Branch("p_value", &p_value, "p_value/D");
0385 
0386     std::string l0_name = tag + "_pull_l0";
0387     std::string l1_name = tag + "_pull_l1";
0388     std::string phi_name = tag + "_pull_phi";
0389     std::string theta_name = tag + "_pull_theta";
0390     std::string qop_name = tag + "_pull_qop";
0391     std::string chi2_name = tag + "_chi2";
0392     std::string pval_name = tag + "_pval";
0393 
0394     TH1D* h_l0 =
0395         new TH1D(l0_name.c_str(), l0_name.c_str(), n_bins, pull_min, pull_max);
0396     TH1D* h_l1 =
0397         new TH1D(l1_name.c_str(), l1_name.c_str(), n_bins, pull_min, pull_max);
0398     TH1D* h_phi = new TH1D(phi_name.c_str(), phi_name.c_str(), n_bins, pull_min,
0399                            pull_max);
0400     TH1D* h_theta = new TH1D(theta_name.c_str(), theta_name.c_str(), n_bins,
0401                              pull_min, pull_max);
0402     TH1D* h_qop = new TH1D(qop_name.c_str(), qop_name.c_str(), n_bins, pull_min,
0403                            pull_max);
0404     TH1D* h_chi2 =
0405         new TH1D(chi2_name.c_str(), chi2_name.c_str(), 50, 0.f, 50.f);
0406     TH1D* h_pval =
0407         new TH1D(pval_name.c_str(), pval_name.c_str(), 50, 0.f, 1.0f);
0408 
0409     // N, mean, stddev, N_error, mean_error, stddev_error, ndf, chi2, ndf/chi2,
0410     // pval, N_4sigma, N_4sigma_fraction, N_expected_4sigma,
0411     // N_4sigma/N_expected_4sigma
0412     std::array<double, 14u> finfo_l0;
0413     std::array<double, 14u> finfo_l1;
0414     std::array<double, 14u> finfo_phi;
0415     std::array<double, 14u> finfo_theta;
0416     std::array<double, 14u> finfo_qop;
0417     finfo_l0[10u] = 0;
0418     finfo_l1[10u] = 0;
0419     finfo_phi[10u] = 0;
0420     finfo_theta[10u] = 0;
0421     finfo_qop[10u] = 0;
0422 
0423     unsigned int n_outliers{0u};
0424 
0425     // Fill the histograms
0426     for (int i = 0; i < t->GetEntries(); i++) {
0427         t->GetEntry(i);
0428         h_l0->Fill(pull_l0);
0429         h_l1->Fill(pull_l1);
0430         h_phi->Fill(pull_phi);
0431         h_theta->Fill(pull_theta);
0432         h_qop->Fill(pull_qop);
0433         h_chi2->Fill(chi2);
0434         p_value = ROOT::Math::chisquared_cdf_c(chi2, 5.f);
0435         h_pval->Fill(p_value);
0436         b_pval->Fill();
0437 
0438         if ((pull_l0 < pull_min) || (pull_l0 > pull_max) ||
0439             (pull_l1 < pull_min) || (pull_l1 > pull_max) ||
0440             (pull_phi < pull_min) || (pull_phi > pull_max) ||
0441             (pull_theta < pull_min) || (pull_theta > pull_max) ||
0442             (pull_qop < pull_min) || (pull_qop > pull_max)) {
0443             n_outliers++;
0444         }
0445 
0446         if (abs(pull_l0) > 4) {
0447             finfo_l0[10u]++;
0448         }
0449         if (abs(pull_l1) > 4) {
0450             finfo_l1[10u]++;
0451         }
0452         if (abs(pull_phi) > 4) {
0453             finfo_phi[10u]++;
0454         }
0455         if (abs(pull_theta) > 4) {
0456             finfo_theta[10u]++;
0457         }
0458         if (abs(pull_qop) > 4) {
0459             finfo_qop[10u]++;
0460         }
0461     }
0462 
0463     t->Write("", TObject::kOverwrite);
0464 
0465     std::clog << tag << " n outliers: " << n_outliers << std::endl;
0466 
0467     finfo_l0[11u] = double(finfo_l0[10u]) / double(t->GetEntries());
0468     finfo_l1[11u] = double(finfo_l1[10u]) / double(t->GetEntries());
0469     finfo_phi[11u] = double(finfo_phi[10u]) / double(t->GetEntries());
0470     finfo_theta[11u] = double(finfo_theta[10u]) / double(t->GetEntries());
0471     finfo_qop[11u] = double(finfo_qop[10u]) / double(t->GetEntries());
0472 
0473     finfo_l0[12u] = t->GetEntries() * 0.000063342484;
0474     finfo_l1[12u] = t->GetEntries() * 0.000063342484;
0475     finfo_phi[12u] = t->GetEntries() * 0.000063342484;
0476     finfo_theta[12u] = t->GetEntries() * 0.000063342484;
0477     finfo_qop[12u] = t->GetEntries() * 0.000063342484;
0478 
0479     finfo_l0[13u] = finfo_l0[10u] / finfo_l0[12u];
0480     finfo_l1[13u] = finfo_l1[10u] / finfo_l1[12u];
0481     finfo_phi[13u] = finfo_phi[10u] / finfo_phi[12u];
0482     finfo_theta[13u] = finfo_theta[10u] / finfo_theta[12u];
0483     finfo_qop[13u] = finfo_qop[10u] / finfo_qop[12u];
0484 
0485     auto c_l0 =
0486         new TCanvas(h_l0->GetName(), h_l0->GetName(), cdim1[0], cdim1[1]);
0487     c_l0->SetLogy();
0488 
0489     draw_pull(h_l0, header_title, geom_title, log10_rk_tolerance, finfo_l0);
0490 
0491     c_l0->SaveAs(to_pdf(l0_name).c_str());
0492 
0493     auto c_l1 =
0494         new TCanvas(h_l1->GetName(), h_l1->GetName(), cdim1[0], cdim1[1]);
0495     c_l1->SetLogy();
0496     draw_pull(h_l1, header_title, geom_title, log10_rk_tolerance, finfo_l1);
0497     c_l1->SaveAs(to_pdf(l1_name).c_str());
0498 
0499     auto c_phi =
0500         new TCanvas(h_phi->GetName(), h_phi->GetName(), cdim1[0], cdim1[1]);
0501     c_phi->SetLogy();
0502     draw_pull(h_phi, header_title, geom_title, log10_rk_tolerance, finfo_phi);
0503     c_phi->SaveAs(to_pdf(phi_name).c_str());
0504 
0505     auto c_theta =
0506         new TCanvas(h_theta->GetName(), h_theta->GetName(), cdim1[0], cdim1[1]);
0507     c_theta->SetLogy();
0508 
0509     draw_pull(h_theta, header_title, geom_title, log10_rk_tolerance,
0510               finfo_theta);
0511     c_theta->SaveAs(to_pdf(theta_name).c_str());
0512 
0513     auto c_qop =
0514         new TCanvas(h_qop->GetName(), h_qop->GetName(), cdim1[0], cdim1[1]);
0515     c_qop->SetLogy();
0516     draw_pull(h_qop, header_title, geom_title, log10_rk_tolerance, finfo_qop);
0517     c_qop->SaveAs(to_pdf(qop_name).c_str());
0518 
0519     auto c_chi2 =
0520         new TCanvas(h_chi2->GetName(), h_chi2->GetName(), cdim1[0], cdim1[1]);
0521     h_chi2->Draw();
0522     c_chi2->SaveAs(to_pdf(chi2_name).c_str());
0523 
0524     auto c_pval =
0525         new TCanvas(h_pval->GetName(), h_pval->GetName(), cdim2[0], cdim2[1]);
0526     draw_pval(h_pval, header_title, geom_title, log10_rk_tolerance);
0527     c_pval->SaveAs(to_pdf(pval_name).c_str());
0528 
0529     std::ofstream f;
0530 
0531     std::string fname = tag + ".txt";
0532     f.open(fname.c_str());
0533 
0534     f << "N, mean, stddev, N_error, mean_error, stddev_error, ndf, chi2, "
0535          "ndf/chi2, pval, N_4sigma, N_4sigma_fraction, N_4sigma_expected, "
0536          "N_4sigma/N_4sigma_expected";
0537     f << std::endl;
0538 
0539     for (const auto& n : finfo_l0) {
0540         f << n << ",";
0541     }
0542     f << std::endl;
0543 
0544     for (const auto& n : finfo_l1) {
0545         f << n << ",";
0546     }
0547     f << std::endl;
0548 
0549     for (const auto& n : finfo_phi) {
0550         f << n << ",";
0551     }
0552     f << std::endl;
0553 
0554     for (const auto& n : finfo_theta) {
0555         f << n << ",";
0556     }
0557     f << std::endl;
0558 
0559     for (const auto& n : finfo_qop) {
0560         f << n << ",";
0561     }
0562     f << std::endl;
0563 
0564     f.close();
0565 }
0566 
0567 // ROOT Script for covariance file reading
0568 void covariance_validation() {
0569 
0570     gStyle->SetOptTitle(0);
0571     gStyle->SetOptStat(0);
0572 
0573     const std::string rect_title = "Bound-to-bound transport";
0574     const std::string wire_title = "Perigee-to-perigee transport";
0575     const std::string geom_title = "RKN with the ODD magnetic field and CsI";
0576 
0577     /************************
0578      *  Rectangular
0579      * **********************/
0580 
0581     std::string rect_name = "rect_cov_transport";
0582     auto rect_tree = get_tree(rect_name);
0583     read_tree(rect_tree, "bound_to_bound", rect_title, geom_title);
0584 
0585     /************************
0586      *  Wire
0587      * **********************/
0588 
0589     std::string wire_name = "wire_cov_transport";
0590     auto wire_tree = get_tree(wire_name);
0591     read_tree(wire_tree, "perigee_to_perigee", wire_title, geom_title);
0592 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!