Skip to content

Commit ed97d20

Browse files
authored
Merge pull request #765 from beomki-yeo/improve-script
Tune the jacobian validation parameters
2 parents b414920 + 42584df commit ed97d20

File tree

7 files changed

+233
-29
lines changed

7 files changed

+233
-29
lines changed

.gitignore

+2
Original file line numberDiff line numberDiff line change
@@ -18,3 +18,5 @@
1818
**/*C_ACLiC*
1919
**/thread_*
2020
**/merged/
21+
**/*.err
22+
**/*.log

core/include/detray/propagator/rk_stepper.ipp

-2
Original file line numberDiff line numberDiff line change
@@ -6,8 +6,6 @@
66
*/
77

88
// Project include(s).
9-
#include <iostream>
10-
119
#include "detray/geometry/tracking_volume.hpp"
1210

1311
template <typename magnetic_field_t, typename algebra_t, typename constraint_t,

tests/integration_tests/cpu/propagator/jacobian_validation.cpp

+7-6
Original file line numberDiff line numberDiff line change
@@ -56,13 +56,13 @@ namespace {
5656
const vector3 B_z{0.f, 0.f, 1.996f * unit<scalar>::T};
5757

5858
// Initial delta for numerical differentiaion
59-
const std::array<scalar, 5u> h_sizes_rect{3e0f, 3e0f, 2e-2f, 1e-3f, 1e-3f};
60-
const std::array<scalar, 5u> h_sizes_wire{3e0f, 3e0f, 2e-2f, 1e-3f, 1e-3f};
59+
const std::array<scalar, 5u> h_sizes_rect{1e0f, 1e0f, 2e-2f, 1e-3f, 1e-3f};
60+
const std::array<scalar, 5u> h_sizes_wire{1e0f, 1e0f, 2e-2f, 1e-3f, 1e-3f};
6161

6262
// Ridders' algorithm setup
6363
constexpr const unsigned int Nt = 100u;
6464
const std::array<scalar, 5u> safe{5.0f, 5.0f, 5.0f, 5.0f, 5.0f};
65-
const std::array<scalar, 5u> con{1.2f, 1.2f, 1.2f, 1.2f, 1.2f};
65+
const std::array<scalar, 5u> con{1.1f, 1.1f, 1.1f, 1.1f, 1.1f};
6666
constexpr const scalar big = std::numeric_limits<scalar>::max();
6767

6868
std::random_device rd;
@@ -180,8 +180,9 @@ struct ridders_derivative {
180180
Arr[j][q][p];
181181
/*
182182
// Please leave this for debug
183-
if (j == e_bound_theta && i == e_bound_loc1) {
184-
std::cout << getter::element(
183+
if (j == e_bound_theta && i == e_bound_loc0) {
184+
std::cout << q << " " << p << " "
185+
<< getter::element(
185186
differentiated_jacobian, j, i)
186187
<< " " << math::abs(Arr[j][q][p])
187188
<< std::endl;
@@ -195,7 +196,7 @@ struct ridders_derivative {
195196
for (unsigned int j = 0; j < 5u; j++) {
196197
/*
197198
// Please leave this for debug
198-
if (j == e_bound_theta && i == e_bound_loc1) {
199+
if (j == e_bound_loc0 && i == e_bound_theta) {
199200
std::cout << getter::element(differentiated_jacobian, j, i)
200201
<< " " << Arr[j][p][p] << " "
201202
<< Arr[j][p - 1][p - 1] << " "

tests/scripts/run_jacobian_validation.sh

+19-11
Original file line numberDiff line numberDiff line change
@@ -66,57 +66,58 @@ while getopts "hd:n:t:m:c:p:q:i:s:f:r:v:" arg; do
6666
;;
6767
d)
6868
dir=$OPTARG
69-
echo "Directory of detray_integration_test_jacobian_validation: ${dir}"
7069
;;
7170
n)
7271
n_threads=$OPTARG
73-
echo "Number of threads: ${n_threads}"
7472
;;
7573
t)
7674
n_tracks_per_thread=$OPTARG
77-
echo "Number of tracks per thread: ${n_tracks_per_thread}"
7875
;;
7976
m)
8077
log10_rk_tol_dis=$OPTARG
8178
echo "log10(rk_error_tolerance_in_mm_for_displaced_tracks): ${log10_rk_tol_dis}"
8279
;;
8380
c)
8481
log10_rk_tol_cov=$OPTARG
85-
echo "log10(rk_error_tolerance_in_mm_for_covariance_transport): ${log10_rk_tol_cov}"
8682
;;
8783
p)
8884
log10_min_rk_tol=$OPTARG
8985
log10_rk_tol_jac=${log10_min_rk_tol}
90-
echo "log10(min_rk_error_tolerance_in_mm): ${log10_min_rk_tol}"
9186
;;
9287
q)
9388
log10_max_rk_tol=$OPTARG
94-
echo "log10(max_rk_error_tolerance_in_mm): ${log10_max_rk_tol}"
9589
;;
9690
i)
9791
log10_helix_tol=$OPTARG
9892
log10_on_surface_tol=$OPTARG
99-
echo "log10(intersection_tolerance_in_mm): ${log10_helix_tol}"
10093
;;
10194
s)
10295
mc_seed=$OPTARG
103-
echo "Monte-Carlo seed: ${mc_seed}"
10496
;;
10597
f)
10698
skip_first_phase=$OPTARG
107-
echo "Skip the first phase: ${skip_first_phase}"
10899
;;
109100
r)
110101
skip_second_phase=$OPTARG
111-
echo "Skip the second phase: ${skip_second_phase}"
112102
;;
113103
v)
114104
verbose_level=$OPTARG
115-
echo "Set the verbose level: ${verbose_level}"
116105
;;
117106
esac
118107
done
119108

109+
echo "Directory of detray_integration_test_jacobian_validation: ${dir}"
110+
echo "Number of threads: ${n_threads}"
111+
echo "Number of tracks per thread: ${n_tracks_per_thread}"
112+
echo "log10(rk_error_tolerance_in_mm_for_covariance_transport): ${log10_rk_tol_cov}"
113+
echo "log10(min_rk_error_tolerance_in_mm): ${log10_min_rk_tol}"
114+
echo "log10(max_rk_error_tolerance_in_mm): ${log10_max_rk_tol}"
115+
echo "log10(intersection_tolerance_in_mm): ${log10_helix_tol}"
116+
echo "Monte-Carlo seed: ${mc_seed}"
117+
echo "Skip the first phase: ${skip_first_phase}"
118+
echo "Skip the second phase: ${skip_second_phase}"
119+
echo "Set the verbose level: ${verbose_level}"
120+
120121
echo ""
121122

122123
if [ -z "${dir}" ]; then
@@ -241,6 +242,9 @@ if [ "$skip_first_phase" = false ] && [ "$skip_second_phase" = false ]; then
241242
# Run rk_tolerance_comparision.C
242243
root -q '../../../tests/validation/root/rk_tolerance_comparison.C+O('${log10_min_rk_tol}','${log10_max_rk_tol}')'
243244

245+
# Run jacobian_histogram.C
246+
root -q -l '../../../tests/validation/root/jacobian_histogram.C+O('${log10_min_rk_tol}')'
247+
244248
# Run jacobian_comparison.C
245249
root -q -l ../../../tests/validation/root/jacobian_comparison.C+O
246250

@@ -256,4 +260,8 @@ if [ "$skip_first_phase" = false ] && [ "$skip_second_phase" = false ]; then
256260

257261
# Run rk_tolerance_comparision.C
258262
root '../../../tests/validation/root/rk_tolerance_comparison.C+O('${log10_min_rk_tol}','${log10_max_rk_tol}')'
263+
264+
# Run jacobian_histogram.C
265+
root '../../../tests/validation/root/jacobian_histogram.C+O('${log10_min_rk_tol}')'
266+
259267
fi

tests/validation/root/covariance_validation.C

+5-5
Original file line numberDiff line numberDiff line change
@@ -165,15 +165,15 @@ void set_xaxis_title(TH1D* h, const double text_size) {
165165
const TString h_name = h->GetName();
166166

167167
if (h_name.Contains("l0")) {
168-
x_axis_title = "#font[12]{PL}(l_{0F})";
168+
x_axis_title = "PL(l_{0F})";
169169
} else if (h_name.Contains("l1")) {
170-
x_axis_title = "#font[12]{PL}(l_{1F})";
170+
x_axis_title = "PL(l_{1F})";
171171
} else if (h_name.Contains("phi")) {
172-
x_axis_title = "#font[12]{PL}(#phi_{F})";
172+
x_axis_title = "PL(#phi_{F})";
173173
} else if (h_name.Contains("theta")) {
174-
x_axis_title = "#font[12]{PL}(#theta_{F})";
174+
x_axis_title = "PL(#theta_{F})";
175175
} else if (h_name.Contains("qop")) {
176-
x_axis_title = "#font[12]{PL}(#lambda_{F})";
176+
x_axis_title = "PL(#lambda_{F})";
177177
} else if (h_name.Contains("pval")) {
178178
x_axis_title = "p-value";
179179
}

tests/validation/root/jacobian_comparison.C

+9-5
Original file line numberDiff line numberDiff line change
@@ -41,7 +41,7 @@ int legend_font = 132;
4141
double legend_font_size = 0.045;
4242
double y_min = -15;
4343
double y_max = 10;
44-
double y_margin = 0;
44+
double y_margin = 1;
4545
double header_size = 0.05;
4646
std::array<float, 4> ldim{0.59015, 0.62395, 0.942404, 0.880252};
4747
double pad_x0 = 0.00;
@@ -132,8 +132,12 @@ void fill_histo(TH1D* hist, const std::array<double, 25u>& means,
132132
hist->SetFillColor(38);
133133
hist->SetCanExtend(TH1::kAllAxes);
134134

135-
for (int i = 0; i < n_labels; i++) {
136-
hist->Fill(labels[i].c_str(), TMath::Log10(means[i]));
135+
for (int i = 0; i < 25; i++) {
136+
if (i < n_labels || i % 6 == 0) {
137+
hist->Fill(labels[i].c_str(), TMath::Log10(means[i]));
138+
} else {
139+
hist->Fill(labels[i].c_str(), 1e20);
140+
}
137141
}
138142

139143
hist->LabelsDeflate();
@@ -233,7 +237,7 @@ void draw_pad(const std::string& pad_name) {
233237
void draw_text(const std::string& text) {
234238

235239
const float x1 = 1.23427;
236-
const float y1 = 6.44408;
240+
const float y1 = 7.37948;
237241

238242
TLatex* ttext = new TLatex(0.f, 0.f, text.c_str());
239243
ttext->SetTextFont(132);
@@ -247,7 +251,7 @@ void draw_text(const std::string& text) {
247251
new TPaveLabel(x1, y1, x1 + float(w) / gPad->GetWw() * 0.62,
248252
y1 + float(h) / gPad->GetWh() * 1.15, text.c_str());
249253

250-
plabel->SetTextFont(132);
254+
plabel->SetTextFont(22);
251255
plabel->SetFillColor(kWhite);
252256
plabel->Draw();
253257
}
+191
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,191 @@
1+
/** Detray library, part of the ACTS project (R&D line)
2+
*
3+
* (c) 2024 CERN for the benefit of the ACTS project
4+
*
5+
* Mozilla Public License Version 2.0
6+
*/
7+
8+
// ROOT include(s).
9+
#include <TCanvas.h>
10+
#include <TFile.h>
11+
#include <TGaxis.h>
12+
#include <TGraph.h>
13+
#include <TLatex.h>
14+
#include <TLegend.h>
15+
#include <TLegendEntry.h>
16+
#include <TMath.h>
17+
#include <TMultiGraph.h>
18+
#include <TROOT.h>
19+
#include <TStyle.h>
20+
21+
#include <ROOT/RCsvDS.hxx>
22+
#include <ROOT/RDataFrame.hxx>
23+
24+
// System include(s).
25+
#include <array>
26+
#include <iostream>
27+
#include <sstream>
28+
#include <string>
29+
#include <vector>
30+
31+
namespace {
32+
33+
double x_pos1 = -13.4f;
34+
double x_pos2 = 0.06f;
35+
double y_pos1 = 2e5;
36+
double y_pos2 = 6.75e4;
37+
38+
// double title_x = x_pos;
39+
// double title_y = 0.8197f;
40+
41+
int label_font = 132;
42+
double label_font_size = 0.055;
43+
double header_text_size = 0.055;
44+
double geom_text_size = 0.0434028;
45+
double titleX_font_size = 0.05;
46+
double titleY_font_size = 0.055;
47+
double x_title_offset = 1.75;
48+
double y_title_offset = 1.34;
49+
double x_label_offset = 0.015;
50+
double y_label_offset = 0.015;
51+
52+
const std::array<float, 2> cdim{700, 600};
53+
double maximum = 1e6;
54+
55+
double pad_x0 = 0.00f;
56+
double pad_x1 = 1.f;
57+
double pad_y0 = 0.00f;
58+
double pad_y1 = 1.f;
59+
60+
double bin_width = 0.2f;
61+
} // namespace
62+
63+
void draw_text(double x1, double y1, double y2, double s1, double s2,
64+
std::string t1, std::string t2) {
65+
TLatex* ttext1 = new TLatex(x1, y1, t1.c_str());
66+
TLatex* ttext2 = new TLatex(x1, y2, t2.c_str());
67+
ttext1->SetTextFont(22);
68+
ttext1->SetTextSize(s1);
69+
ttext2->SetTextFont(132);
70+
ttext2->SetTextSize(s2);
71+
72+
ttext1->Draw();
73+
ttext2->Draw();
74+
}
75+
76+
void histo_setup(TH1D* histo) {
77+
histo->GetXaxis()->SetLabelFont(label_font);
78+
histo->GetYaxis()->SetLabelFont(label_font);
79+
histo->GetXaxis()->SetLabelSize(label_font_size);
80+
histo->GetYaxis()->SetLabelSize(label_font_size);
81+
histo->GetXaxis()->SetTitleSize(titleX_font_size);
82+
histo->GetYaxis()->SetTitleSize(titleY_font_size);
83+
histo->GetYaxis()->SetTitleOffset(y_title_offset);
84+
histo->GetXaxis()->SetTitleOffset(x_title_offset + 0.1);
85+
histo->GetXaxis()->SetLabelOffset(x_label_offset);
86+
histo->GetYaxis()->SetLabelOffset(y_label_offset);
87+
histo->GetYaxis()->SetNdivisions(504);
88+
histo->GetYaxis()->SetMaxDigits(1);
89+
histo->GetXaxis()->CenterTitle(true);
90+
histo->GetYaxis()->CenterTitle(true);
91+
histo->GetXaxis()->SetTitleFont(132);
92+
histo->GetYaxis()->SetTitleFont(132);
93+
}
94+
95+
void set_yaxis_title(TH1D* h, const double text_size) {
96+
double bin_width = h->GetBinWidth(0u);
97+
std::string str = std::to_string(bin_width);
98+
str.erase(str.find_last_not_of('0') + 1, std::string::npos);
99+
str.erase(str.find_last_not_of('.') + 1, std::string::npos);
100+
std::string y_axis_title = "Counts / (" + str + ")";
101+
h->GetYaxis()->SetTitle(y_axis_title.c_str());
102+
h->GetYaxis()->SetTitleSize(text_size);
103+
}
104+
105+
void draw_histogram(const std::string root_name, const int num) {
106+
107+
const std::string rect_title = "Bound-to-bound transport";
108+
const std::string geom_title = "RKN with the ODD magnetic field and CsI";
109+
110+
TFile* f = TFile::Open(root_name.c_str(), "read");
111+
TTree* t = (TTree*)f->Get("inhom_rect_material");
112+
113+
auto dthetadl0_canvas =
114+
new TCanvas("dthetadl0_canvas", "dthetadl0_canvas", cdim[0], cdim[1]);
115+
dthetadl0_canvas->SetLogy();
116+
117+
TPad* pad1 = new TPad("pad1", "pad1", pad_x0, pad_y0, pad_x1, pad_y1);
118+
pad1->Draw();
119+
pad1->cd();
120+
pad1->SetLeftMargin(110. / pad1->GetWw());
121+
pad1->SetBottomMargin(125. / pad1->GetWh());
122+
pad1->SetLogy();
123+
124+
t->Draw("log10(abs(dthetadl0_E)) >> htemp(100,-14,-4)");
125+
TH1D* dthetadl0_hist = (TH1D*)gPad->GetPrimitive("htemp");
126+
histo_setup(dthetadl0_hist);
127+
set_yaxis_title(dthetadl0_hist, titleY_font_size);
128+
dthetadl0_hist->GetXaxis()->SetNdivisions(505);
129+
dthetadl0_hist->GetXaxis()->SetTitle(
130+
"log_{10}( #left|#frac{#partial#theta_{F}}{#partiall_{0I}}#right| "
131+
"[mm^{-1}] )");
132+
dthetadl0_hist->SetMaximum(maximum);
133+
dthetadl0_hist->SetLineColor(kOrange + 3);
134+
dthetadl0_hist->SetFillColor(kOrange + 2);
135+
dthetadl0_hist->SetFillStyle(3001);
136+
137+
draw_text(x_pos1, y_pos1, y_pos2, header_text_size, geom_text_size,
138+
rect_title.c_str(), geom_title.c_str());
139+
140+
dthetadl0_canvas->Draw();
141+
dthetadl0_canvas->SaveAs("bound_to_bound_dthetadl0_E_histo.pdf");
142+
143+
auto dqopdqop_canvas =
144+
new TCanvas("dqopdqop_canvas", "dqopdqop_canvas", cdim[0], cdim[1]);
145+
dqopdqop_canvas->SetLogy();
146+
147+
TPad* pad2 = new TPad("pad2", "pad2", pad_x0, pad_y0, pad_x1, pad_y1);
148+
pad2->Draw();
149+
pad2->cd();
150+
pad2->SetLeftMargin(110. / pad2->GetWw());
151+
pad2->SetBottomMargin(125. / pad2->GetWh());
152+
pad2->SetLogy();
153+
154+
t->Draw("log10(dqopdqop_E) >> htemp2(100,0,1)");
155+
TH1D* dqopdqop_hist = (TH1D*)gPad->GetPrimitive("htemp2");
156+
histo_setup(dqopdqop_hist);
157+
set_yaxis_title(dqopdqop_hist, titleY_font_size);
158+
159+
histo_setup(dqopdqop_hist);
160+
dqopdqop_hist->GetXaxis()->SetNdivisions(505);
161+
dqopdqop_hist->GetXaxis()->SetTitle(
162+
"log_{10}( "
163+
"#left|#frac{#partial#lambda_{F}}{#partial#lambda_{I}}#right| )");
164+
dqopdqop_hist->SetMaximum(maximum);
165+
dqopdqop_hist->SetLineColor(kGreen + 3);
166+
dqopdqop_hist->SetFillColor(kGreen + 2);
167+
dqopdqop_hist->SetFillStyle(3001);
168+
169+
draw_text(x_pos2, y_pos1, y_pos2, header_text_size, geom_text_size,
170+
rect_title.c_str(), geom_title.c_str());
171+
172+
dqopdqop_canvas->Draw();
173+
dqopdqop_canvas->SaveAs("bound_to_bound_dqopdqop_E_histo.pdf");
174+
}
175+
176+
void jacobian_histogram(int num) {
177+
gStyle->SetOptTitle(0);
178+
gStyle->SetOptStat(0);
179+
180+
const std::string name = "inhom_rect_material_" + std::to_string(num);
181+
const std::string csv_name = name + ".csv";
182+
const std::string root_name = "residual_histogram.root";
183+
184+
std::cout << "Processing file: " << csv_name << std::endl;
185+
186+
auto rdf = ROOT::RDF::FromCSV(csv_name);
187+
// Create root file
188+
rdf.Snapshot("inhom_rect_material", root_name.c_str());
189+
190+
draw_histogram(root_name, num);
191+
}

0 commit comments

Comments
 (0)