tools/root/rk_tolerance_comparison.C

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 <TCanvas.h>
0011 #include <TGaxis.h>
0012 #include <TGraph.h>
0013 #include <TLatex.h>
0014 #include <TLegend.h>
0015 #include <TLegendEntry.h>
0016 #include <TMath.h>
0017 #include <TMultiGraph.h>
0018 #include <TROOT.h>
0019 #include <TStyle.h>
0020
0021 #include <ROOT/RCsvDS.hxx>
0022 #include <ROOT/RDataFrame.hxx>
0023
0024 // System include(s).
0025 #include <array>
0026 #include <iostream>
0027 #include <sstream>
0028 #include <string>
0029 #include <vector>
0030
0031 namespace {
0032 const double step_title_x_offset = 1.4;
0033 const double step_title_x = 0.175;
0034 const double step_title_y = 0.835;
0035 const double step_ygap = -0.0505;
0036 const double x_label_offset = 0.015;
0037 const double x_margin = 1.;
0038 const double y_label_offset = 0.015;
0039 const double rk_title_x_offset = 0.9;
0040 const double rk_title_offset_fraction = 0.0218;
0041 const double rk_title_y = 14.69;
0042 const double rk_ygap = -0.83;
0043 const double rk_header_text_size = 0.046;
0044 const double rk_geom_text_size = 0.0362903;
0045 const double step_header_text_size = 0.055;
0046 const double step_geom_text_size = 0.0434028;
0047 const double label_font_size_step = 0.055;
0048 const double title_font_size_step = 0.055;
0049 const double label_font_size_rk_tol = 0.046;
0050 const double title_font_size_rk_tol = 0.046;
0051 const int title_font = 132;
0052 const int label_font = 132;
0053 const int legend_font = 132;
0054 const double pad_x0 = 0.005;
0055 const double pad_x1 = 1.;
0056 const double pad_y0 = 0.005;
0057 const double pad_y1 = 1.;
0058 const double ymin = -2.;
0059 const double ymax = 4.;
0060 std::array<double, 4u> ldim{0.20985, 0.533, 0.87997, 0.939};
0061 }  // namespace
0062
0063 std::vector<std::string> create_labels() {
0064
0065     std::vector<std::string> varI = {
0066         "{#partiall_{0I}}", "{#partiall_{1I}}", "{#partial#phi_{I}}",
0067         "{#partial#theta_{I}}", "{#partial#lambda_{I}}"};
0068
0069     std::vector<std::string> varF = {
0070         "{#partiall_{0F}}", "{#partiall_{1F}}", "{#partial#phi_{F}}",
0071         "{#partial#theta_{F}}", "{#partial#lambda_{F}}"};
0072
0073     std::vector<std::string> labels;
0074
0075     for (int i = 0; i < 5; i++) {
0076         for (int j = 0; j < 5; j++) {
0077             std::string head = "#frac";
0078             labels.push_back(head + varF[i] + varI[j]);
0079         }
0080     }
0081
0082     return labels;
0083 }
0084
0085 void draw_text(double x1, double y1, double y_delta, double s1, double s2,
0086                std::string t1, std::string t2) {
0087     TLatex* ttext1 = new TLatex(x1, y1, t1.c_str());
0088     TLatex* ttext2 = new TLatex(x1, y1 + y_delta, t2.c_str());
0089     ttext1->SetTextFont(22);
0090     ttext1->SetTextSize(s1);
0091     ttext2->SetTextFont(132);
0092     ttext2->SetTextSize(s2);
0093
0094     ttext1->Draw();
0095     ttext2->Draw();
0096 }
0097
0098 std::vector<std::string> create_columns(const std::string& tag) {
0099     std::vector<std::string> varI = {"dl0", "dl1", "dphi", "dtheta", "dqop"};
0100
0101     std::vector<std::string> varF = {"dl0", "dl1", "dphi", "dtheta", "dqop"};
0102
0103     std::vector<std::string> columns;
0104
0105     for (int i = 0; i < 5; i++) {
0106         for (int j = 0; j < 5; j++) {
0107             columns.push_back(varF[i] + varI[j] + "_" + tag);
0108         }
0109     }
0110
0111     return columns;
0112 }
0113
0114 std::array<double, 25u> get_means(ROOT::RDataFrame& rdf) {
0115
0116     std::array<double, 25u> ret;
0117
0118     auto col_names_R = create_columns("R");
0119
0120     for (std::size_t i = 0u; i < 25u; i++) {
0121         const double residual_mean = *rdf.Mean<double>(col_names_R[i]);
0122         ret[i] = residual_mean;
0123     }
0124
0125     return ret;
0126 }
0127
0128 std::map<std::string, std::vector<double>> get_means(
0129     const std::vector<std::string> labels, const std::string tag, const int min,
0130     const int max, std::vector<double>& mean_step_sizes) {
0131
0132     std::map<std::string, std::vector<double>> ret;
0133
0134     int num = min;
0135     while (num <= max + 1e-3) {
0136
0137         const std::string name = tag + "_" + std::to_string(num);
0138         const std::string csv_name = name + ".csv";
0139
0140         const std::string sign = num >= 0 ? "p" : "m";
0141         const std::string root_name = name + ".root";
0142
0143         std::clog << "Processing file: " << csv_name << std::endl;
0144
0145         auto rdf = ROOT::RDF::FromCSV(csv_name);
0146
0147         const std::array<double, 25u> means = get_means(rdf);
0148
0149         for (unsigned int i = 0; i < 25u; i++) {
0150             ret[labels[i]].push_back(TMath::Log10(means[i]));
0151         }
0152
0153         mean_step_sizes.push_back(
0154             TMath::Log10(*rdf.Mean<double>("average_step_size")));
0155
0156         // Create root file
0157         rdf.Snapshot(tag.c_str(), root_name.c_str());
0158
0159         num = num + 2;
0160     }
0161
0162     return ret;
0163 }
0164
0165 std::vector<double> get_x_vector(const int min, const int max) {
0166     std::vector<double> ret;
0167
0168     int num = min;
0169     while (num <= max + 1e-3) {
0170         ret.push_back(num);
0171         num = num + 2;
0172     }
0173
0174     return ret;
0175 }
0176
0177 void draw_graphs(const std::string header_title, const std::string geom_title,
0178                  const std::vector<std::string> labels,
0179                  const std::vector<double> x_vec,
0180                  std::map<std::string, std::vector<double>> means) {
0181
0182     TPad* gr_pad = new TPad("gr_pad", "gr_pad", 0, 0, 1, 1);
0183     gr_pad->Draw();
0184     gr_pad->cd();
0185     gr_pad->SetLeftMargin(122. / gr_pad->GetWw());
0186     gr_pad->SetTopMargin(50. / gr_pad->GetWh());
0187     gr_pad->SetBottomMargin(120. / gr_pad->GetWh());
0188
0189     TGraph* gr[25];
0190     TMultiGraph* mg = new TMultiGraph();
0191
0192     const std::array<int, 5u> marker_styles = {7, 2, 5, 27, 25};
0193     const std::array<double, 5u> marker_sizes = {2.135, 2.135, 2.135, 2.135,
0194                                                  1.49};
0195     const std::array<int, 5u> line_styles = {1, 3, 2, 7, 4};
0196     const std::array<int, 5u> hues = {kOrange + 2, kPink + 5, kBlue + 2,
0197                                       kCyan + 2, kGreen + 2};
0198
0199     auto legend = new TLegend(ldim[0u], ldim[1u], ldim[2u], ldim[3u]);
0200
0201     // legend->SetHeader(header_title.c_str());
0202     legend->SetHeader("");
0203     legend->SetNColumns(5);
0204     legend->SetColumnSeparation(-0.3);
0205     legend->SetFillStyle(0);
0206     legend->SetMargin(0.55);
0207     legend->SetTextFont(legend_font);
0208
0209     for (int i = 0; i < 25; i++) {
0210         gr[i] = new TGraph(x_vec.size(), &x_vec[0], &means[labels[i]][0]);
0211
0212         const int n = i / 5;
0213         const int m = i % 5;
0214
0215         gr[i]->SetMarkerStyle(marker_styles[n]);
0216         gr[i]->SetMarkerSize(marker_sizes[n]);
0217         gr[i]->SetLineStyle(line_styles[m]);
0218         gr[i]->SetMarkerColor(hues[m]);
0219         gr[i]->SetLineColor(hues[m]);
0220
0221         mg->Add(gr[i]);
0222         legend->AddEntry(gr[i], labels[i].c_str(), "lp");
0223     }
0224
0225     mg->GetXaxis()->SetTitle("log_{10}(#font[12]{#tau} [mm])");
0226     mg->GetXaxis()->SetLabelOffset(-0.005);
0227     mg->GetXaxis()->SetTitleOffset(rk_title_x_offset);
0228     mg->GetXaxis()->SetTitleSize(title_font_size_rk_tol);
0229     mg->GetXaxis()->SetLabelSize(label_font_size_rk_tol);
0230     mg->GetXaxis()->CenterTitle(true);
0231     mg->GetYaxis()->CenterTitle(true);
0232
0233     double yaxis_min = -10;
0234     double yaxis_max = 15;
0235     double yaxis_margin = 1.;
0236
0237     std::clog << "Vec size: " << x_vec.size() << std::endl;
0238     mg->GetYaxis()->SetTitle("log_{10}(Mean of #font[12]{#Omega_{R}})");
0239     mg->GetYaxis()->SetTitleOffset(1.48);
0240     mg->GetYaxis()->SetTitleSize(title_font_size_rk_tol);
0241     mg->GetYaxis()->SetLabelSize(label_font_size_rk_tol);
0242     mg->GetYaxis()->SetRangeUser(yaxis_min - yaxis_margin,
0243                                  yaxis_max + yaxis_margin);
0244     mg->GetYaxis()->SetLabelOffset(0.01);
0245     mg->GetYaxis()->SetTitleFont(title_font);
0246     mg->GetXaxis()->SetTitleFont(title_font);
0247
0248     double x_min = x_vec.front();
0249     double x_max = x_vec.back();
0250
0251     if (x_vec.size() > 10) {
0252         if (x_vec.size() == 13) {
0253             x_margin = 2;
0254         }
0255
0256         x_min = x_min - x_margin;
0257         x_max = x_max + x_margin;
0258         mg->GetXaxis()->SetLimits(x_min, x_max);
0259         mg->GetXaxis()->SetLabelSize(0);
0260         mg->GetXaxis()->SetTickLength(0);
0261
0262         mg->Draw("APL");
0263
0264         if (x_vec.size() == 11) {
0265             auto ga = new TGaxis(x_vec.front(), yaxis_min - yaxis_margin,
0266                                  x_vec.back(), yaxis_min - yaxis_margin,
0267                                  x_vec.front(), x_vec.back(), 405, "N");
0268             ga->SetLabelFont(label_font);
0269             ga->SetLabelSize(label_font_size_rk_tol);
0270             ga->SetLabelOffset(-0.0065);
0271             ga->Draw();
0272         } else if (x_vec.size() == 13) {
0273             auto ga = new TGaxis(x_vec.front(), yaxis_min - yaxis_margin,
0274                                  x_vec.back(), yaxis_min - yaxis_margin,
0275                                  x_vec.front(), x_vec.back(), 304, "N");
0276             ga->SetLabelFont(label_font);
0277             ga->SetLabelSize(label_font_size_rk_tol);
0278             ga->SetLabelOffset(-0.0065);
0279             ga->Draw();
0280         }
0281
0282         mg->GetYaxis()->SetLabelSize(0);
0283         mg->GetYaxis()->SetTickLength(0);
0284         auto ga_y = new TGaxis(x_min, yaxis_min, x_min, yaxis_max, yaxis_min,
0285                                yaxis_max, 505, "N");
0286         ga_y->SetLabelFont(label_font);
0287         ga_y->SetLabelOffset(0.02);
0288         ga_y->SetLabelSize(label_font_size_rk_tol);
0289         ga_y->Draw();
0290
0291     } else {
0292         mg->Draw("APL");
0293     }
0294
0295     TLegendEntry* header =
0296         (TLegendEntry*)legend->GetListOfPrimitives()->First();
0297     header->SetTextFont(22);
0298     header->SetTextSize(.033);
0299
0300     double rk_title_deltaX =
0301         (x_vec.back() - x_vec.front()) * rk_title_offset_fraction;
0302
0303     if (x_vec.size() == 13) {
0304         rk_title_deltaX = -0.2;
0305     }
0306
0307     legend->Draw();
0308     draw_text(rk_title_deltaX + x_vec.front(), rk_title_y, rk_ygap,
0309               rk_header_text_size, rk_geom_text_size, header_title, geom_title);
0310 }
0311
0312 void draw_mean_step_size(const std::string header_title,
0313                          const std::string geom_title,
0314                          const std::vector<double>& x_vec,
0315                          const std::vector<double>& means) {
0316     TPad* step_pad =
0317         new TPad("step_pad", "step_pad", pad_x0, pad_y0, pad_x1, pad_y1);
0318     step_pad->Draw();
0319     step_pad->cd();
0320
0321     step_pad->SetLeftMargin(93. / step_pad->GetWw());
0322     step_pad->SetBottomMargin(90. / step_pad->GetWh());
0323     step_pad->SetTopMargin(50. / step_pad->GetWh());
0324
0325     TGraph* gr = new TGraph(x_vec.size(), &x_vec[0], &means[0]);
0326
0327     gr->SetMarkerStyle(4);
0328     gr->SetMarkerSize(1.2);
0329     gr->GetXaxis()->SetTitle("log_{10}(#font[12]{#tau} [mm])");
0330     gr->GetYaxis()->SetTitle("log_{10}(Mean of avg. step size [mm])");
0331     gr->GetXaxis()->SetLimits(x_vec.front() - 0.5, x_vec.back() + 0.5);
0332
0333     if (x_vec.size() == 13) {
0334         ymin = -4;
0335     }
0336
0337     gr->GetYaxis()->SetRangeUser(ymin, ymax);
0338     gr->GetYaxis()->SetNdivisions(505);
0339     gr->GetXaxis()->SetLabelSize(label_font_size_step);
0340     gr->GetYaxis()->SetLabelSize(label_font_size_step);
0341     gr->GetXaxis()->SetTitleSize(title_font_size_step);
0342     gr->GetYaxis()->SetTitleSize(title_font_size_step);
0343     gr->GetXaxis()->SetTitleOffset(step_title_x_offset);
0344     gr->GetYaxis()->SetTitleOffset(1.15);
0345     gr->GetYaxis()->SetLabelOffset(y_label_offset);
0346     gr->GetYaxis()->SetTitleFont(title_font);
0347     gr->GetXaxis()->SetTitleFont(title_font);
0348     gr->GetXaxis()->SetLabelOffset(x_label_offset);
0349     gr->GetXaxis()->CenterTitle(true);
0350     gr->GetYaxis()->CenterTitle(true);
0351     gr->GetXaxis()->SetLabelFont(label_font);
0352     gr->GetYaxis()->SetLabelFont(label_font);
0353
0354     if (x_vec.size() > 10) {
0355         if (x_vec.size() == 13) {
0356             x_margin = 2;
0357         }
0358
0359         gr->GetXaxis()->SetLimits(x_vec.front() - x_margin,
0360                                   x_vec.back() + x_margin);
0361         gr->GetXaxis()->SetLabelSize(0);
0362         gr->GetXaxis()->SetTickLength(0);
0363
0364         gr->Draw("APL");
0365
0366         if (x_vec.size() == 11) {
0367             auto ga = new TGaxis(x_vec.front(), ymin, x_vec.back(), ymin,
0368                                  x_vec.front(), x_vec.back(), 405, "N");
0369             ga->SetLabelFont(label_font);
0370             ga->SetLabelSize(label_font_size_step);
0371             ga->SetLabelOffset(x_label_offset);
0372             ga->Draw();
0373         } else if (x_vec.size() == 13) {
0374             auto ga = new TGaxis(x_vec.front(), ymin, x_vec.back(), ymin,
0375                                  x_vec.front(), x_vec.back(), 304, "N");
0376             ga->SetLabelFont(label_font);
0377             ga->SetLabelSize(label_font_size_step);
0378             ga->SetLabelOffset(x_label_offset);
0379             ga->Draw();
0380         }
0381
0382     } else {
0383         gr->Draw();
0384     }
0385
0386     TPad* text_pad =
0387         new TPad("text_pad", "text_pad", pad_x0, pad_y0, pad_x1, pad_y1);
0388     text_pad->SetFillStyle(4000);
0389     text_pad->Draw();
0390     text_pad->cd();
0391
0392     draw_text(step_title_x, step_title_y, step_ygap, step_header_text_size,
0393               step_geom_text_size, header_title, geom_title);
0394 }
0395
0396 // ROOT Script for jacboain file reading
0397 void rk_tolerance_comparison(int min, int max) {
0398     gStyle->SetOptTitle(0);
0399     gStyle->SetLegendBorderSize(0);
0400     gStyle->SetLegendTextSize(0.0333);
0401
0402     const std::array<float, 2> cdim1{800, 1300};
0403     const std::array<float, 2> cdim2{700, 600};
0404
0405     auto labels = create_labels();
0406
0407     const auto x_vec = get_x_vector(min, max);
0408
0409     const std::string rect_header = "Bound-to-bound transport";
0410     const std::string wire_header = "Perigee-to-perigee transport";
0411     const std::string geom_header = "RKN with the ODD magnetic field and CsI";
0412
0413     /************************
0414      *  Rectangular
0415      * **********************/
0416
0417     const std::string rect_pdf = "bound_to_bound_rk_tolerance.pdf";
0418
0419     auto rect_canvas =
0420         new TCanvas("rect_canvas", "rect_canvas", cdim1[0], cdim1[1]);
0421
0422     std::vector<double> rect_mean_step_sizes;
0423     const auto rect_y_means = get_means(labels, "inhom_rect_material", min, max,
0424                                         rect_mean_step_sizes);
0425     draw_graphs(rect_header, geom_header, labels, x_vec, rect_y_means);
0426
0427     rect_canvas->SaveAs(rect_pdf.c_str());
0428
0429     auto rect_canvas2 =
0430         new TCanvas("rect_canvas2", "rect_canvas2", cdim2[0], cdim2[1]);
0431     const std::string rect_mean_step_pdf = "bound_to_bound_mean_step_size.pdf";
0432
0433     draw_mean_step_size(rect_header, geom_header, x_vec, rect_mean_step_sizes);
0434     rect_canvas2->SaveAs(rect_mean_step_pdf.c_str());
0435
0436     /************************
0437      *  Wire
0438      * **********************/
0439
0440     const std::string wire_pdf = "perigee_to_perigee_rk_tolerance.pdf";
0441
0442     auto wire_canvas =
0443         new TCanvas("wire_canvas", "wire_canvas", cdim1[0], cdim1[1]);
0444     std::vector<double> wire_mean_step_sizes;
0445     const auto wire_y_means = get_means(labels, "inhom_wire_material", min, max,
0446                                         wire_mean_step_sizes);
0447     draw_graphs(wire_header, geom_header, labels, x_vec, wire_y_means);
0448
0449     wire_canvas->SaveAs(wire_pdf.c_str());
0450
0451     auto wire_canvas2 =
0452         new TCanvas("wire_canvas2", "wire_canvas2", cdim2[0], cdim2[1]);
0453     const std::string wire_mean_step_pdf =
0454         "perigee_to_perigee_mean_step_size.pdf";
0455
0456     draw_mean_step_size(wire_header, geom_header, x_vec, wire_mean_step_sizes);
0457     wire_canvas2->SaveAs(wire_mean_step_pdf.c_str());
0458 }