tools/root/jacobian_histogram.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 <TFile.h>
0012 #include <TGaxis.h>
0013 #include <TGraph.h>
0014 #include <TLatex.h>
0015 #include <TLegend.h>
0016 #include <TLegendEntry.h>
0017 #include <TMath.h>
0018 #include <TMultiGraph.h>
0019 #include <TROOT.h>
0020 #include <TStyle.h>
0021
0022 #include <ROOT/RCsvDS.hxx>
0023 #include <ROOT/RDataFrame.hxx>
0024
0025 // System include(s).
0026 #include <array>
0027 #include <iostream>
0028 #include <sstream>
0029 #include <string>
0030 #include <vector>
0031
0032 namespace {
0033
0034 const double x_pos1 = -13.4;
0035 const double x_pos2 = 0.06;
0036 const double y_pos1 = 2e5;
0037 const double y_pos2 = 6.75e4;
0038
0039 // const double title_x = x_pos;
0040 // const double title_y = 0.8197;
0041
0042 const int label_font = 132;
0043 const double label_font_size = 0.055;
0044 const double header_text_size = 0.055;
0045 const double geom_text_size = 0.0434028;
0046 const double titleX_font_size = 0.05;
0047 const double titleY_font_size = 0.055;
0048 const double x_title_offset = 1.75;
0049 const double y_title_offset = 1.34;
0050 const double x_label_offset = 0.015;
0051 const double y_label_offset = 0.015;
0052 const int fill_style = 1001;
0053 const double color_alpha = 0.5;
0054
0055 const std::array<float, 2> cdim{700, 600};
0056 const double maximum = 1e6;
0057
0058 const double pad_x0 = 0.00;
0059 const double pad_x1 = 1.;
0060 const double pad_y0 = 0.00;
0061 const double pad_y1 = 1.;
0062
0063 const double bin_width = 0.2;
0064 }  // namespace
0065
0066 void draw_text(double x1, double y1, double y2, double s1, double s2,
0067                std::string t1, std::string t2) {
0068     TLatex* ttext1 = new TLatex(x1, y1, t1.c_str());
0069     TLatex* ttext2 = new TLatex(x1, y2, t2.c_str());
0070     ttext1->SetTextFont(22);
0071     ttext1->SetTextSize(s1);
0072     ttext2->SetTextFont(132);
0073     ttext2->SetTextSize(s2);
0074
0075     ttext1->Draw();
0076     ttext2->Draw();
0077 }
0078
0079 void histo_setup(TH1D* histo) {
0080     histo->GetXaxis()->SetLabelFont(label_font);
0081     histo->GetYaxis()->SetLabelFont(label_font);
0082     histo->GetXaxis()->SetLabelSize(label_font_size);
0083     histo->GetYaxis()->SetLabelSize(label_font_size);
0084     histo->GetXaxis()->SetTitleSize(titleX_font_size);
0085     histo->GetYaxis()->SetTitleSize(titleY_font_size);
0086     histo->GetYaxis()->SetTitleOffset(y_title_offset);
0087     histo->GetXaxis()->SetTitleOffset(x_title_offset + 0.1);
0088     histo->GetXaxis()->SetLabelOffset(x_label_offset);
0089     histo->GetYaxis()->SetLabelOffset(y_label_offset);
0090     histo->GetYaxis()->SetNdivisions(504);
0091     histo->GetYaxis()->SetMaxDigits(1);
0092     histo->GetXaxis()->CenterTitle(true);
0093     histo->GetYaxis()->CenterTitle(true);
0094     histo->GetXaxis()->SetTitleFont(132);
0095     histo->GetYaxis()->SetTitleFont(132);
0096 }
0097
0098 void set_yaxis_title(TH1D* h, const double text_size) {
0099     double bin_width = h->GetBinWidth(0u);
0100     std::string str = std::to_string(bin_width);
0101     str.erase(str.find_last_not_of('0') + 1, std::string::npos);
0102     str.erase(str.find_last_not_of('.') + 1, std::string::npos);
0103     std::string y_axis_title = "Counts / (" + str + ")";
0104     h->GetYaxis()->SetTitle(y_axis_title.c_str());
0105     h->GetYaxis()->SetTitleSize(text_size);
0106 }
0107
0108 void draw_histogram(const std::string root_name, const int num) {
0109
0110     const std::string rect_title = "Bound-to-bound transport";
0111     const std::string geom_title = "RKN with the ODD magnetic field and CsI";
0112
0113     TFile* f = TFile::Open(root_name.c_str(), "read");
0114     TTree* t = (TTree*)f->Get("inhom_rect_material");
0115
0116     auto dthetadl0_canvas =
0117         new TCanvas("dthetadl0_canvas", "dthetadl0_canvas", cdim[0], cdim[1]);
0118     dthetadl0_canvas->SetLogy();
0119
0120     TPad* pad1 = new TPad("pad1", "pad1", pad_x0, pad_y0, pad_x1, pad_y1);
0121     pad1->Draw();
0122     pad1->cd();
0123     pad1->SetLeftMargin(110. / pad1->GetWw());
0124     pad1->SetBottomMargin(125. / pad1->GetWh());
0125     pad1->SetLogy();
0126
0127     t->Draw("log10(abs(dthetadl0_E)) >> htemp(100,-14,-4)");
0128     TH1D* dthetadl0_hist = (TH1D*)gPad->GetPrimitive("htemp");
0129     histo_setup(dthetadl0_hist);
0130     set_yaxis_title(dthetadl0_hist, titleY_font_size);
0131     dthetadl0_hist->GetXaxis()->SetNdivisions(505);
0132     dthetadl0_hist->GetXaxis()->SetTitle(
0133         "log_{10}( #left|#frac{#partial#theta_{F}}{#partiall_{0I}}#right| "
0134         "[mm^{-1}] )");
0135     dthetadl0_hist->SetMaximum(maximum);
0136     dthetadl0_hist->SetLineColor(kOrange + 3);
0137     dthetadl0_hist->SetFillStyle(fill_style);
0138     dthetadl0_hist->SetFillColorAlpha(kOrange + 2, color_alpha);
0139
0140     draw_text(x_pos1, y_pos1, y_pos2, header_text_size, geom_text_size,
0141               rect_title.c_str(), geom_title.c_str());
0142
0143     dthetadl0_canvas->Draw();
0144     dthetadl0_canvas->SaveAs("bound_to_bound_dthetadl0_E_histo.pdf");
0145
0146     auto dqopdqop_canvas =
0147         new TCanvas("dqopdqop_canvas", "dqopdqop_canvas", cdim[0], cdim[1]);
0148     dqopdqop_canvas->SetLogy();
0149
0150     TPad* pad2 = new TPad("pad2", "pad2", pad_x0, pad_y0, pad_x1, pad_y1);
0151     pad2->Draw();
0152     pad2->cd();
0153     pad2->SetLeftMargin(110. / pad2->GetWw());
0154     pad2->SetBottomMargin(125. / pad2->GetWh());
0155     pad2->SetLogy();
0156
0157     t->Draw("log10(dqopdqop_E) >> htemp2(100,0,1)");
0158     TH1D* dqopdqop_hist = (TH1D*)gPad->GetPrimitive("htemp2");
0159     histo_setup(dqopdqop_hist);
0160     set_yaxis_title(dqopdqop_hist, titleY_font_size);
0161
0162     histo_setup(dqopdqop_hist);
0163     dqopdqop_hist->GetXaxis()->SetNdivisions(505);
0164     dqopdqop_hist->GetXaxis()->SetTitle(
0165         "log_{10}( "
0166         "#left|#frac{#partial#lambda_{F}}{#partial#lambda_{I}}#right| )");
0167     dqopdqop_hist->SetMaximum(maximum);
0168     dqopdqop_hist->SetLineColor(kGreen + 3);
0169     dqopdqop_hist->SetFillStyle(fill_style);
0170     dqopdqop_hist->SetFillColorAlpha(kGreen + 2, color_alpha);
0171
0172     draw_text(x_pos2, y_pos1, y_pos2, header_text_size, geom_text_size,
0173               rect_title.c_str(), geom_title.c_str());
0174
0175     dqopdqop_canvas->Draw();
0176     dqopdqop_canvas->SaveAs("bound_to_bound_dqopdqop_E_histo.pdf");
0177 }
0178
0179 void jacobian_histogram(int num) {
0180     gStyle->SetOptTitle(0);
0181     gStyle->SetOptStat(0);
0182
0183     const std::string name = "inhom_rect_material_" + std::to_string(num);
0184     const std::string csv_name = name + ".csv";
0185     const std::string root_name = "residual_histogram.root";
0186
0187     std::clog << "Processing file: " << csv_name << std::endl;
0188
0189     auto rdf = ROOT::RDF::FromCSV(csv_name);
0190     // Create root file
0191     rdf.Snapshot("inhom_rect_material", root_name.c_str());
0192
0193     draw_histogram(root_name, num);
0194 }