tools/root/jacobian_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 <TLatex.h>
0013 #include <TLegend.h>
0014 #include <TLegendEntry.h>
0015 #include <TLine.h>
0016 #include <TMath.h>
0017 #include <TPaveLabel.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 <iostream>
0026 #include <limits>
0027 #include <sstream>
0028 #include <vector>
0029
0030 namespace {
0031 const double labelx_font_size = 0.055;
0032 const double labelx_offset = 0.0065;
0033 const double labely_font_size = 0.055;
0034 const double labely_offset = 0.01;
0035 const double title_font_size = 0.055;
0036 const double title_offset = 0.71;
0037 const double marker_size = 1.3875;
0038 const double legend_margin = 0.12;
0039 const int title_font = 132;
0040 const int label_font = 132;
0041 const int legend_font = 132;
0042 const double legend_font_size = 0.045;
0043 const double y_min = -15;
0044 const double y_max = 10;
0045 const double y_margin = 1;
0046 const double header_size = 0.05;
0047 const std::array<float, 4> ldim{0.59015, 0.62395, 0.942404, 0.880252};
0048 const double pad_x0 = 0.00;
0049 const double pad_x1 = 1;
0050 const double pad_y0 = 0.00;
0051 const double pad_y1 = 1;
0052
0053 }  // namespace
0054
0055 std::vector<std::string> create_labels() {
0056
0057     std::vector<std::string> varI = {
0058         "{#partiall_{0I}}", "{#partiall_{1I}}", "{#partial#phi_{I}}",
0059         "{#partial#theta_{I}}", "{#partial#lambda_{I}}"};
0060
0061     std::vector<std::string> varF = {
0062         "{#partiall_{0F}}", "{#partiall_{1F}}", "{#partial#phi_{F}}",
0063         "{#partial#theta_{F}}", "{#partial#lambda_{F}}"};
0064
0065     std::vector<std::string> labels;
0066
0067     for (int i = 0; i < 5; i++) {
0068         for (int j = 0; j < 5; j++) {
0069             std::string head = "#frac";
0070             labels.push_back(head + varF[i] + varI[j]);
0071         }
0072     }
0073
0074     return labels;
0075 }
0076
0077 std::vector<std::string> create_columns(const std::string& tag) {
0078     std::vector<std::string> varI = {"dl0", "dl1", "dphi", "dtheta", "dqop"};
0079
0080     std::vector<std::string> varF = {"dl0", "dl1", "dphi", "dtheta", "dqop"};
0081
0082     std::vector<std::string> columns;
0083
0084     for (int i = 0; i < 5; i++) {
0085         for (int j = 0; j < 5; j++) {
0086             columns.push_back(varF[i] + varI[j] + "_" + tag);
0087         }
0088     }
0089
0090     return columns;
0091 }
0092
0093 std::array<double, 25u> get_means(ROOT::RDataFrame& rdf) {
0094
0095     // Residual columns
0096     auto col_names_R = create_columns("R");
0097
0098     // Evaluate columns
0099     auto col_names_E = create_columns("E");
0100
0101     std::array<double, 25u> ret;
0102
0103     for (std::size_t i = 0u; i < 25u; i++) {
0104
0105         const double residual_mean = *rdf.Mean<double>(col_names_R[i]);
0106         const double evaluate_mean = *rdf.Mean<double>(col_names_E[i]);
0107
0108         // dtheta/d(l0,l1,phi,qop) is analytically zero for homogeneous field
0109         // dqop(dl0,dl1,dphi,dtheta) is analytically zero for free space without
0110         // material so those should not appear on the plots
0111         if (std::abs(evaluate_mean) < 1e-20f) {
0112             ret[i] = std::numeric_limits<double>::min();
0113         } else {
0114             ret[i] = residual_mean;
0115         }
0116     }
0117
0118     return ret;
0119 }
0120
0121 TString get_legend_header(const double log10_rk_tol) {
0122     std::stringstream val_stream;
0123     val_stream << "#tau = 10^{" << int(log10_rk_tol)
0124                << "} mm as an RKN error tolerance";
0125
0126     return TString(val_stream.str());
0127 }
0128
0129 void fill_histo(TH1D* hist, const std::array<double, 25u>& means,
0130                 const std::vector<string>& labels, const int n_labels,
0131                 const int marker_style) {
0132     hist->SetStats(0);
0133     hist->SetFillColor(38);
0134     hist->SetCanExtend(TH1::kAllAxes);
0135
0136     for (int i = 0; i < 25; i++) {
0137         if (i < n_labels || i % 6 == 0) {
0138             hist->Fill(labels[i].c_str(), TMath::Log10(means[i]));
0139         } else {
0140             hist->Fill(labels[i].c_str(), 1e20);
0141         }
0142     }
0143
0144     hist->LabelsDeflate();
0145     hist->SetMarkerSize(marker_size);
0146     hist->SetMarkerStyle(marker_style);
0147 }
0148
0149 void draw_lines() {
0150     for (int i = 1; i < 25; i++) {
0151         if (i % 5 == 0) {
0152             TLine* line = new TLine(i, y_min - y_margin, i, y_max + y_margin);
0153             line->SetLineColor(kBlack);
0154             line->SetLineWidth(1);
0155             line->Draw();
0156         } else {
0157             TLine* line = new TLine(i, y_min - y_margin, i, y_max + y_margin);
0158             line->SetLineColor(kBlack);
0159             line->SetLineWidth(1);
0160             line->SetLineStyle(3);
0161             line->Draw();
0162         }
0163     }
0164 }
0165
0166 TH1D* get_histogram(std::string name, const int n_labels,
0167                     const int marker_style, double& log10_rk_tol) {
0168
0169     std::clog << "Generating histogram for " << name << std::endl;
0170
0171     auto labels = create_labels();
0172
0173     const std::string csv_name = name + ".csv";
0174     const std::string root_name = name + ".root";
0175     const std::string histo_name = name + "_histo";
0176
0177     auto rdf = ROOT::RDF::FromCSV(csv_name);
0178     auto rdf_means = get_means(rdf);
0179
0180     // Count the number of non-convergence event
0181     auto conv_success = rdf.Filter("total_convergence == 1").Count();
0182     auto conv_failure = rdf.Filter("total_convergence == 0").Count();
0183
0184     std::clog << "Convergence events: " << *conv_success << std::endl;
0185     std::clog << "Non-convergence events: " << *conv_failure << std::endl;
0186
0187     TH1D* histo = new TH1D(histo_name.c_str(), histo_name.c_str(), 3, 0, 3);
0188     histo->GetYaxis()->SetRangeUser(y_min - y_margin, y_max + y_margin);
0189     histo->GetYaxis()->SetLabelSize(0);
0190     histo->GetYaxis()->SetTickLength(0);
0191
0192     histo->GetYaxis()->SetTitle("log_{10}(Mean of #font[12]{#Omega_{R}})");
0193     histo->GetYaxis()->SetTitleOffset(0.68);
0194     histo->GetYaxis()->SetTitleSize(title_font_size);
0195     histo->GetYaxis()->SetTitleFont(title_font);
0196     histo->GetYaxis()->SetTitleOffset(title_offset);
0197     histo->GetXaxis()->SetLabelSize(labelx_font_size);
0198     histo->GetXaxis()->SetLabelOffset(labelx_offset);
0199     histo->GetXaxis()->SetLabelFont(label_font);
0200     histo->GetYaxis()->CenterTitle(true);
0201
0202     if (TString(name).Contains("helix")) {
0203         auto ga_y = new TGaxis(0, y_min, 0, y_max, y_min, y_max, 505, "N");
0204         ga_y->SetLabelFont(42);
0205         ga_y->SetLabelOffset(labely_offset);
0206         ga_y->SetLabelSize(labely_font_size);
0207         ga_y->SetLabelFont(label_font);
0208         ga_y->Draw();
0209     }
0210
0211     fill_histo(histo, rdf_means, labels, n_labels, marker_style);
0212
0213     // No color
0214     histo->SetFillColor(kWhite);
0215
0216     // Create root file
0217     rdf.Snapshot(name, root_name);
0218
0219     // Get log10(rk_tolerance)
0220     if (!TString(name).Contains("helix")) {
0221         log10_rk_tol = *rdf.Mean<double>("log10_rk_tolerance");
0222     }
0223
0224     return histo;
0225 }
0226
0227 void draw_pad(const std::string& pad_name) {
0228     TPad* apad = new TPad(pad_name.c_str(), pad_name.c_str(), pad_x0, pad_y0,
0229                           pad_x1, pad_y1);
0230     apad->Draw();
0231     apad->cd();
0232
0233     apad->SetLeftMargin(105. / apad->GetWw());
0234     apad->SetRightMargin(60. / apad->GetWw());
0235     apad->SetBottomMargin(60. / apad->GetWh());
0236 }
0237
0238 void draw_text(const std::string& text) {
0239
0240     const float x1 = 1.23427;
0241     const float y1 = 7.37948;
0242
0243     TLatex* ttext = new TLatex(0.f, 0.f, text.c_str());
0244     ttext->SetTextFont(132);
0245     ttext->SetTextSize(3.);
0246
0247     UInt_t w;
0248     UInt_t h;
0249     ttext->GetBoundingBox(w, h);
0250
0251     TPaveLabel* plabel =
0252         new TPaveLabel(x1, y1, x1 + float(w) / gPad->GetWw() * 0.62,
0253                        y1 + float(h) / gPad->GetWh() * 1.15, text.c_str());
0254
0255     plabel->SetTextFont(22);
0256     plabel->SetFillColor(kWhite);
0257     plabel->Draw();
0258 }
0259
0260 // ROOT Script for jacboain file reading
0261 void jacobian_comparison() {
0262
0263     gStyle->SetOptTitle(0);
0264     const std::array<float, 2> cdim{1200, 500};
0265     const std::array<int, 4> markers{kOpenCross, kOpenSquare, kOpenTriangleUp,
0266                                      kOpenCircle};
0267     const std::array<int, 4> hues{kOrange + 8, kMagenta + 1, kAzure,
0268                                   kGreen + 2};
0269     const std::string rect_pdf = "bound_to_bound_jacobian_comparison.pdf";
0270     const std::string wire_pdf = "perigee_to_perigee_jacobian_comparison.pdf";
0271
0272     double log10_rk_tolerance_rect;
0273     double log10_rk_tolerance_wire;
0274     double dummy;
0275
0276     /************************
0277      *  Rectangular
0278      * **********************/
0279
0280     auto rect_canvas =
0281         new TCanvas("rect_canvas", "rect_canvas", cdim[0], cdim[1]);
0282     draw_pad("rect_pad");
0283
0284     auto rect_legend = new TLegend(ldim[0], ldim[1], ldim[2], ldim[3]);
0285     rect_legend->SetMargin(legend_margin);
0286     rect_legend->SetTextFont(legend_font);
0287     rect_legend->SetTextSize(legend_font_size);
0288     rect_legend->SetBorderSize(4);
0289
0290     const std::string RKN_ODD_CsI = "RKN with the ODD magnetic field and CsI";
0291     const std::string RKN_ODD = "RKN with the ODD magnetic field";
0292     const std::string RKN_homogeneous =
0293         "RKN with the homogeneous magnetic field";
0294     const std::string Helix_homogeneous =
0295         "Helix with the homogeneous magnetic field";
0296
0297     std::string rect_text = "Bound-to-bound transport";
0298     auto inhom_rect_material_histo = get_histogram(
0299         "inhom_rect_material", 25, markers[0u], log10_rk_tolerance_rect);
0300     // rect_legend->SetHeader("Bound-to-bound transport");
0301     inhom_rect_material_histo->SetMarkerColor(hues[0u]);
0302     inhom_rect_material_histo->Draw("hist P ");
0303     rect_legend->AddEntry(inhom_rect_material_histo, RKN_ODD_CsI.c_str(), "p");
0304
0305     auto inhom_rect_histo = get_histogram("inhom_rect", 20, markers[1u], dummy);
0306     inhom_rect_histo->SetMarkerColor(hues[1u]);
0307     inhom_rect_histo->Draw("hist P same");
0308     rect_legend->AddEntry(inhom_rect_histo, RKN_ODD.c_str(), "p");
0309
0310     auto const_rect_histo = get_histogram("const_rect", 15, markers[2u], dummy);
0311     const_rect_histo->SetMarkerColor(hues[2u]);
0312     const_rect_histo->Draw("hist P same");
0313     rect_legend->AddEntry(const_rect_histo, RKN_homogeneous.c_str(), "p");
0314
0315     auto helix_rect_histo = get_histogram("helix_rect", 15, markers[3u], dummy);
0316     helix_rect_histo->SetMarkerColor(hues[3u]);
0317     helix_rect_histo->Draw("hist P same");
0318     rect_legend->AddEntry(helix_rect_histo, Helix_homogeneous.c_str(), "p");
0319
0320     draw_lines();
0321     rect_legend->Draw();
0322     draw_text(rect_text);
0323     rect_canvas->Draw();
0324
0325     rect_canvas->SaveAs(rect_pdf.c_str());
0326
0327     /************************
0328      *  Wire
0329      * **********************/
0330
0331     auto wire_canvas =
0332         new TCanvas("wire_canvas", "wire_canvas", cdim[0], cdim[1]);
0333     draw_pad("wire_pad");
0334
0335     auto wire_legend = new TLegend(ldim[0], ldim[1], ldim[2], ldim[3]);
0336     wire_legend->SetMargin(legend_margin);
0337     wire_legend->SetTextFont(legend_font);
0338     wire_legend->SetTextSize(legend_font_size);
0339     wire_legend->SetBorderSize(4);
0340
0341     std::string wire_text = "Perigee-to-perigee transport";
0342     auto inhom_wire_material_histo = get_histogram(
0343         "inhom_wire_material", 25, markers[0u], log10_rk_tolerance_wire);
0344     // wire_legend->SetHeader("Perigee-to-perigee transport");
0345     inhom_wire_material_histo->SetMarkerColor(hues[0u]);
0346     inhom_wire_material_histo->Draw("hist P ");
0347     wire_legend->AddEntry(inhom_wire_material_histo, RKN_ODD_CsI.c_str(), "p");
0348
0349     auto inhom_wire_histo = get_histogram("inhom_wire", 20, markers[1u], dummy);
0350     inhom_wire_histo->SetMarkerColor(hues[1u]);
0351     inhom_wire_histo->Draw("hist P same");
0352     wire_legend->AddEntry(inhom_wire_histo, RKN_ODD.c_str(), "p");
0353
0354     auto const_wire_histo = get_histogram("const_wire", 15, markers[2u], dummy);
0355     const_wire_histo->SetMarkerColor(hues[2u]);
0356     const_wire_histo->Draw("hist P same");
0357     wire_legend->AddEntry(const_wire_histo, RKN_homogeneous.c_str(), "p");
0358
0359     auto helix_wire_histo = get_histogram("helix_wire", 15, markers[3u], dummy);
0360     helix_wire_histo->SetMarkerColor(hues[3u]);
0361     helix_wire_histo->Draw("hist P same");
0362     wire_legend->AddEntry(helix_wire_histo, Helix_homogeneous.c_str(), "p");
0363
0364     draw_lines();
0365     wire_legend->Draw();
0366     draw_text(wire_text);
0367     wire_canvas->Draw();
0368
0369     wire_canvas->SaveAs(wire_pdf.c_str());
0370 }