17 #include <TDirectory.h>
25 #include <TProfile2D.h>
46 bool filtered =
true,
bool smoothed =
true,
47 bool fitFiltered =
true,
bool fitPredicted =
true,
48 bool fitSmoothed =
true,
51 gStyle->SetOptFit(0000);
52 gStyle->SetOptStat(0000);
53 gStyle->SetPadLeftMargin(0.20);
54 gStyle->SetPadRightMargin(0.1);
55 gStyle->SetPadTopMargin(0.1);
56 gStyle->SetPadBottomMargin(0.15);
68 }
else if (filtered) {
76 std::cout <<
"Opening file: " << inFile << std::endl;
77 TFile*
file = TFile::Open(inFile.c_str(),
"read");
80 if (file ==
nullptr) {
84 std::cout <<
"Reading tree: " << treeName << std::endl;
85 TTree*
tree = (TTree*)file->Get(treeName.c_str());
88 if (tree ==
nullptr) {
92 TFile*
out = TFile::Open(outFile.c_str(),
"recreate");
99 std::map<int, std::vector<int>> volLayIds;
100 TCanvas* geometryCanvas =
101 new TCanvas(
"volLayCanvas",
"Volume Layer Matrix", 10, 10, 450, 600);
102 geometryCanvas->Divide(1, 2);
104 geometryCanvas->cd(1);
105 tree->Draw(
"layer_id:volume_id>>volID_layID(100,0.5,100.5,100,0.5,100.5)",
"",
107 auto h2_volID_layID =
dynamic_cast<TH2F*
>(out->Get(
"volID_layID"));
109 h2_volID_layID->Write();
111 geometryCanvas->cd(2);
113 rz_draw_string += range_tag;
114 rz_draw_string +=
"*g_x_";
115 rz_draw_string += range_tag;
116 rz_draw_string +=
"+g_y_";
117 rz_draw_string += range_tag;
118 rz_draw_string +=
"*g_y_";
119 rz_draw_string += range_tag;
120 rz_draw_string +=
")):g_z_";
121 rz_draw_string += range_tag;
122 rz_draw_string +=
">>geo_dim";
123 tree->Draw(rz_draw_string.c_str());
124 auto h2_geo_dim =
dynamic_cast<TH2F*
>(out->Get(
"geo_dim"));
126 float detectorZ = 1.15 * h2_geo_dim->GetXaxis()->GetXmax();
127 float detectorR = 1.15 * h2_geo_dim->GetYaxis()->GetXmax();
131 if (not saveAs.empty()) {
132 geometryCanvas->SaveAs((
std::string(
"all_vol_lay_ids.") + saveAs).c_str());
136 int volBins = h2_volID_layID->GetXaxis()->GetNbins();
137 int layBins = h2_volID_layID->GetYaxis()->GetNbins();
139 volLayIds[-1] = {-1};
141 for (
int volID = 1; volID <= volBins; ++volID) {
142 for (
int layID = 1; layID <= layBins; ++layID) {
143 if (h2_volID_layID->GetBinContent(volID, layID) != 0.) {
144 if (volLayIds.find(volID) == volLayIds.end()) {
146 volLayIds[volID] = {-1, layID};
148 volLayIds[volID].push_back(layID);
161 TCanvas* rangeCanvas =
162 new TCanvas(
"rangeCanvas",
"Range Estimation", 10, 10, 900, 600);
163 rangeCanvas->Divide(3, 2);
165 std::vector<std::string> res_ranges = {
std::string(
"res_eLOC0_") + range_tag,
175 auto volLayIdCut = [](
int vol,
int lay) -> std::array<std::string, 2> {
176 if (vol < 0 and lay < 0) {
187 return {vlstr, vlcut};
194 [&](
const std::array<std::string, 2>& vlIdCut) -> std::array<float, 6> {
195 std::array<float, 6> ranges = {0., 0., 0., 0., 0., 0.};
196 for (
unsigned int ir = 0; ir < 6; ++ir) {
197 rangeCanvas->cd(ir + 1);
200 vlIdCut[0] +
std::string(
"range_") + res_ranges[ir];
201 tree->Draw((drawString + histString).c_str(), vlIdCut[1].c_str());
202 auto h1_range =
dynamic_cast<TH1F*
>(out->Get(histString.c_str()));
204 float range = pullRange * h1_range->GetRMS();
207 if (not saveAs.empty()) {
209 (vlIdCut[0] +
std::string(
"res_ranges.") + saveAs).c_str());
215 std::vector<std::string> paramNames = {
"loc0",
"loc1",
"#phi",
216 "#theta",
"q/p",
"t"};
221 std::array<TProfile2D*, 6> p2d_res_zr_prt{};
222 std::array<TProfile2D*, 6> p2d_res_zr_flt{};
223 std::array<TProfile2D*, 6> p2d_res_zr_smt{};
224 std::array<TProfile2D*, 6> p2d_pull_zr_prt{};
225 std::array<TProfile2D*, 6> p2d_pull_zr_flt{};
226 std::array<TProfile2D*, 6> p2d_pull_zr_smt{};
228 for (
unsigned int ipar = 0; ipar < paramNames.size(); ++ipar) {
229 const auto& par = paramNames[ipar];
232 p2d_res_zr_prt[ipar] =
233 new TProfile2D(Form(
"all_prof_res_prt_%s", par.c_str()),
234 Form(
"residual profile of %s", par.c_str()), 100,
235 -detectorZ, detectorZ, 50, 0., detectorR);
237 p2d_pull_zr_prt[ipar] =
238 new TProfile2D(Form(
"all_prof_pull_prt_%s", par.c_str()),
239 Form(
"pull profile of %s", par.c_str()), 100,
240 -detectorZ, detectorZ, 50, 0., detectorR);
242 p2d_res_zr_prt[ipar]->SetErrorOption(
"s");
243 p2d_res_zr_prt[ipar]->GetXaxis()->SetTitle(
"z [mm]");
244 p2d_res_zr_prt[ipar]->GetYaxis()->SetTitle(
"R [mm]");
245 p2d_res_zr_prt[ipar]->GetZaxis()->SetTitle(Form(
"r_{%s}", par.c_str()));
248 p2d_pull_zr_prt[ipar]->SetErrorOption(
"s");
249 p2d_pull_zr_prt[ipar]->GetXaxis()->SetTitle(
"z [mm]");
250 p2d_pull_zr_prt[ipar]->GetYaxis()->SetTitle(
"R [mm]");
251 p2d_pull_zr_prt[ipar]->GetZaxis()->SetTitle(
252 Form(
"pull_{%s}", par.c_str()));
257 p2d_res_zr_flt[ipar] =
258 new TProfile2D(Form(
"all_prof_res_flt_%s", par.c_str()),
259 Form(
"residual profile of %s", par.c_str()), 100,
260 -detectorZ, detectorZ, 50, 0., detectorR);
261 p2d_pull_zr_flt[ipar] =
262 new TProfile2D(Form(
"all_prof_pull_flt_%s", par.c_str()),
263 Form(
"pull profile of %s", par.c_str()), 100,
264 -detectorZ, detectorZ, 50, 0., detectorR);
266 p2d_res_zr_flt[ipar]->SetErrorOption(
"s");
267 p2d_res_zr_flt[ipar]->GetXaxis()->SetTitle(
"z [mm]");
268 p2d_res_zr_flt[ipar]->GetYaxis()->SetTitle(
"R [mm]");
269 p2d_res_zr_flt[ipar]->GetZaxis()->SetTitle(Form(
"r_{%s}", par.c_str()));
272 p2d_pull_zr_flt[ipar]->SetErrorOption(
"s");
273 p2d_pull_zr_flt[ipar]->GetXaxis()->SetTitle(
"z [mm]");
274 p2d_pull_zr_flt[ipar]->GetYaxis()->SetTitle(
"R [mm]");
275 p2d_pull_zr_flt[ipar]->GetZaxis()->SetTitle(
276 Form(
"pull_{%s}", par.c_str()));
281 p2d_res_zr_smt[ipar] =
282 new TProfile2D(Form(
"all_prof_res_smt_%s", par.c_str()),
283 Form(
"residual profile of %s", par.c_str()), 100,
284 -detectorZ, detectorZ, 50, 0., detectorR);
286 p2d_pull_zr_smt[ipar] =
287 new TProfile2D(Form(
"all_prof_pull_smt_%s", par.c_str()),
288 Form(
"pull profile of %s", par.c_str()), 100,
289 -detectorZ, detectorZ, 50, 0., detectorR);
291 p2d_res_zr_smt[ipar]->SetErrorOption(
"s");
292 p2d_pull_zr_smt[ipar]->SetErrorOption(
"s");
294 p2d_res_zr_smt[ipar]->GetXaxis()->SetTitle(
"z [mm]");
295 p2d_res_zr_smt[ipar]->GetYaxis()->SetTitle(
"R [mm]");
296 p2d_res_zr_smt[ipar]->GetZaxis()->SetTitle(Form(
"r_{%s}", par.c_str()));
299 p2d_pull_zr_smt[ipar]->GetXaxis()->SetTitle(
"z [mm]");
300 p2d_pull_zr_smt[ipar]->GetYaxis()->SetTitle(
"R [mm]");
301 p2d_pull_zr_smt[ipar]->GetZaxis()->SetTitle(
302 Form(
"pull_{%s}", par.c_str()));
308 std::map<std::string, TH1F*> res_prt;
309 std::map<std::string, TH1F*> res_flt;
310 std::map<std::string, TH1F*> res_smt;
311 std::map<std::string, TH1F*> pull_prt;
312 std::map<std::string, TH1F*> pull_flt;
313 std::map<std::string, TH1F*> pull_smt;
318 for (
auto [vol, layers] : volLayIds) {
319 for (
auto lay : layers) {
321 auto vlIdCut = volLayIdCut(vol, lay);
322 auto ranges = histRanges(vlIdCut);
325 std::map<std::pair<std::string, std::string>,
double> paramResidualRange =
326 {{{
"loc0",
"l_{0}"}, ranges[0]}, {{
"loc1",
"l_{1}"}, ranges[1]},
327 {{
"#phi",
"#phi"}, ranges[2]}, {{
"#theta",
"#theta"}, ranges[3]},
328 {{
"q/p",
"q/p"}, ranges[4]}, {{
"t",
"t"}, ranges[5]}};
331 for (
const auto& [partwin, resRange] : paramResidualRange) {
336 TString par_string(partwin.second.c_str());
338 par_string + TString(
"^{rec} - ") + par_string + TString(
"^{true}");
340 TString pull_string = TString(
"(") + res_string +
341 TString(
")/#sigma_{") + par_string + TString(
"}");
345 new TH1F(Form((vlIdCut[0] +
std::string(
"res_prt_%s")).c_str(),
347 Form(
"residual of %s", par.c_str()), 100, -1 * resRange,
349 res_prt[id_par]->GetXaxis()->SetTitle(res_string.Data());
350 res_prt[id_par]->GetYaxis()->SetTitle(
"Entries");
353 pull_prt[id_par] =
new TH1F(
354 Form((vlIdCut[0] +
std::string(
"pull_prt_%s")).c_str(),
356 Form(
"pull of %s", par.c_str()), 100, -1 * pullRange, pullRange);
357 pull_prt[id_par]->GetXaxis()->SetTitle(pull_string.Data());
358 pull_prt[id_par]->GetYaxis()->SetTitle(
"Arb. Units");
364 new TH1F(Form((vlIdCut[0] +
std::string(
"res_flt_%s")).c_str(),
366 Form(
"residual of %s", par.c_str()), 100, -1 * resRange,
368 res_flt[id_par]->GetXaxis()->SetTitle(res_string.Data());
369 res_flt[id_par]->GetYaxis()->SetTitle(
"Entries");
372 pull_flt[id_par] =
new TH1F(
373 Form((vlIdCut[0] +
std::string(
"pull_flt_%s")).c_str(),
375 Form(
"pull of %s", par.c_str()), 100, -1 * pullRange, pullRange);
376 pull_flt[id_par]->GetXaxis()->SetTitle(pull_string.Data());
377 pull_flt[id_par]->GetYaxis()->SetTitle(
"Arb. Units");
383 new TH1F(Form((vlIdCut[0] +
std::string(
"res_smt_%s")).c_str(),
385 Form(
"residual of %s", par.c_str()), 100, -1 * resRange,
388 res_smt[id_par]->GetXaxis()->SetTitle(res_string.Data());
389 res_smt[id_par]->GetYaxis()->SetTitle(
"Entries");
392 pull_smt[id_par] =
new TH1F(
393 Form((vlIdCut[0] +
std::string(
"pull_smt_%s")).c_str(),
395 Form(
"pull of %s", par.c_str()), 100, -1 * pullRange, pullRange);
397 pull_smt[id_par]->GetXaxis()->SetTitle(pull_string.Data());
398 pull_smt[id_par]->GetYaxis()->SetTitle(
"Arb. Units");
408 int entries = tree->GetEntries();
409 for (
int j = 0;
j < entries;
j++) {
414 if (predicted and tsReader.
predicted->at(
i)) {
415 auto x_prt = tsReader.
g_x_prt->at(
i);
416 auto y_prt = tsReader.
g_y_prt->at(
i);
417 auto r_prt = std::hypot(x_prt, y_prt);
418 auto z_prt = tsReader.
g_z_prt->at(
i);
419 p2d_res_zr_prt[0]->Fill(z_prt, r_prt, tsReader.
res_LOC0_prt->at(
i));
420 p2d_res_zr_prt[1]->Fill(z_prt, r_prt, tsReader.
res_LOC1_prt->at(
i));
421 p2d_res_zr_prt[2]->Fill(z_prt, r_prt, tsReader.
res_PHI_prt->at(
i));
422 p2d_res_zr_prt[3]->Fill(z_prt, r_prt, tsReader.
res_THETA_prt->at(
i));
423 p2d_res_zr_prt[4]->Fill(z_prt, r_prt, tsReader.
res_QOP_prt->at(
i));
424 p2d_res_zr_prt[5]->Fill(z_prt, r_prt, tsReader.
res_T_prt->at(
i));
425 p2d_pull_zr_prt[0]->Fill(z_prt, r_prt, tsReader.
pull_LOC0_prt->at(
i));
426 p2d_pull_zr_prt[1]->Fill(z_prt, r_prt, tsReader.
pull_LOC1_prt->at(
i));
427 p2d_pull_zr_prt[2]->Fill(z_prt, r_prt, tsReader.
pull_PHI_prt->at(
i));
428 p2d_pull_zr_prt[3]->Fill(z_prt, r_prt, tsReader.
pull_THETA_prt->at(
i));
429 p2d_pull_zr_prt[4]->Fill(z_prt, r_prt, tsReader.
pull_QOP_prt->at(
i));
430 p2d_pull_zr_prt[5]->Fill(z_prt, r_prt, tsReader.
pull_T_prt->at(
i));
432 if (filtered and tsReader.
filtered->at(
i)) {
433 auto x_flt = tsReader.
g_x_flt->at(
i);
434 auto y_flt = tsReader.
g_y_flt->at(
i);
435 auto r_flt = std::hypot(x_flt, y_flt);
436 auto z_flt = tsReader.
g_z_flt->at(
i);
437 p2d_res_zr_flt[0]->Fill(z_flt, r_flt, tsReader.
res_LOC0_flt->at(
i));
438 p2d_res_zr_flt[1]->Fill(z_flt, r_flt, tsReader.
res_LOC1_flt->at(
i));
439 p2d_res_zr_flt[2]->Fill(z_flt, r_flt, tsReader.
res_PHI_flt->at(
i));
440 p2d_res_zr_flt[3]->Fill(z_flt, r_flt, tsReader.
res_THETA_flt->at(
i));
441 p2d_res_zr_flt[4]->Fill(z_flt, r_flt, tsReader.
res_QOP_flt->at(
i));
442 p2d_res_zr_flt[5]->Fill(z_flt, r_flt, tsReader.
res_T_flt->at(
i));
443 p2d_pull_zr_flt[0]->Fill(z_flt, r_flt, tsReader.
pull_LOC0_flt->at(
i));
444 p2d_pull_zr_flt[1]->Fill(z_flt, r_flt, tsReader.
pull_LOC1_flt->at(
i));
445 p2d_pull_zr_flt[2]->Fill(z_flt, r_flt, tsReader.
pull_PHI_flt->at(
i));
446 p2d_pull_zr_flt[3]->Fill(z_flt, r_flt, tsReader.
pull_THETA_flt->at(
i));
447 p2d_pull_zr_flt[4]->Fill(z_flt, r_flt, tsReader.
pull_QOP_flt->at(
i));
448 p2d_pull_zr_flt[5]->Fill(z_flt, r_flt, tsReader.
pull_T_flt->at(
i));
450 if (smoothed and tsReader.
smoothed->at(
i)) {
451 auto x_smt = tsReader.
g_x_smt->at(
i);
452 auto y_smt = tsReader.
g_y_smt->at(
i);
453 auto r_smt = std::hypot(x_smt, y_smt);
454 auto z_smt = tsReader.
g_z_smt->at(
i);
455 p2d_res_zr_smt[0]->Fill(z_smt, r_smt, tsReader.
res_LOC0_smt->at(
i));
456 p2d_res_zr_smt[1]->Fill(z_smt, r_smt, tsReader.
res_LOC1_smt->at(
i));
457 p2d_res_zr_smt[2]->Fill(z_smt, r_smt, tsReader.
res_PHI_smt->at(
i));
458 p2d_res_zr_smt[3]->Fill(z_smt, r_smt, tsReader.
res_THETA_smt->at(
i));
459 p2d_res_zr_smt[4]->Fill(z_smt, r_smt, tsReader.
res_QOP_smt->at(
i));
460 p2d_res_zr_smt[5]->Fill(z_smt, r_smt, tsReader.
res_T_smt->at(
i));
461 p2d_pull_zr_smt[0]->Fill(z_smt, r_smt, tsReader.
pull_LOC0_smt->at(
i));
462 p2d_pull_zr_smt[1]->Fill(z_smt, r_smt, tsReader.
pull_LOC1_smt->at(
i));
463 p2d_pull_zr_smt[2]->Fill(z_smt, r_smt, tsReader.
pull_PHI_smt->at(
i));
464 p2d_pull_zr_smt[3]->Fill(z_smt, r_smt, tsReader.
pull_THETA_smt->at(
i));
465 p2d_pull_zr_smt[4]->Fill(z_smt, r_smt, tsReader.
pull_QOP_smt->at(
i));
466 p2d_pull_zr_smt[5]->Fill(z_smt, r_smt, tsReader.
pull_T_smt->at(
i));
473 std::vector<std::array<int, 2>> fillIds = {
474 {-1, -1}, {vol, -1}, {vol, lay}};
476 for (
const auto&
fid : fillIds) {
477 auto vlID = volLayIdCut(
fid[0],
fid[1])[0];
479 if (predicted and tsReader.
predicted->at(
i)) {
480 res_prt[vlID + paramNames[0]]->Fill(tsReader.
res_LOC0_prt->at(
i), 1);
481 res_prt[vlID + paramNames[1]]->Fill(tsReader.
res_LOC1_prt->at(
i), 1);
482 res_prt[vlID + paramNames[2]]->Fill(tsReader.
res_PHI_prt->at(
i), 1);
483 res_prt[vlID + paramNames[3]]->Fill(tsReader.
res_THETA_prt->at(
i), 1);
484 res_prt[vlID + paramNames[4]]->Fill(tsReader.
res_QOP_prt->at(
i), 1);
485 res_prt[vlID + paramNames[5]]->Fill(tsReader.
res_T_prt->at(
i), 1);
486 pull_prt[vlID + paramNames[0]]->Fill(tsReader.
pull_LOC0_prt->at(
i),
488 pull_prt[vlID + paramNames[1]]->Fill(tsReader.
pull_LOC1_prt->at(
i),
490 pull_prt[vlID + paramNames[2]]->Fill(tsReader.
pull_PHI_prt->at(
i), 1);
493 pull_prt[vlID + paramNames[4]]->Fill(tsReader.
pull_QOP_prt->at(
i), 1);
494 pull_prt[vlID + paramNames[5]]->Fill(tsReader.
pull_T_prt->at(
i), 1);
497 if (filtered and tsReader.
filtered->at(
i)) {
498 res_flt[vlID + paramNames[0]]->Fill(tsReader.
res_LOC0_flt->at(
i), 1);
499 res_flt[vlID + paramNames[1]]->Fill(tsReader.
res_LOC1_flt->at(
i), 1);
500 res_flt[vlID + paramNames[2]]->Fill(tsReader.
res_PHI_flt->at(
i), 1);
501 res_flt[vlID + paramNames[3]]->Fill(tsReader.
res_THETA_flt->at(
i), 1);
502 res_flt[vlID + paramNames[4]]->Fill(tsReader.
res_QOP_flt->at(
i), 1);
503 res_flt[vlID + paramNames[5]]->Fill(tsReader.
res_T_flt->at(
i), 1);
504 pull_flt[vlID + paramNames[0]]->Fill(tsReader.
pull_LOC0_flt->at(
i),
506 pull_flt[vlID + paramNames[1]]->Fill(tsReader.
pull_LOC1_flt->at(
i),
508 pull_flt[vlID + paramNames[2]]->Fill(tsReader.
pull_PHI_flt->at(
i), 1);
511 pull_flt[vlID + paramNames[4]]->Fill(tsReader.
pull_QOP_flt->at(
i), 1);
512 pull_flt[vlID + paramNames[5]]->Fill(tsReader.
pull_T_flt->at(
i), 1);
515 if (smoothed and tsReader.
smoothed->at(
i)) {
516 res_smt[vlID + paramNames[0]]->Fill(tsReader.
res_LOC0_smt->at(
i), 1);
517 res_smt[vlID + paramNames[1]]->Fill(tsReader.
res_LOC1_smt->at(
i), 1);
518 res_smt[vlID + paramNames[2]]->Fill(tsReader.
res_PHI_smt->at(
i), 1);
519 res_smt[vlID + paramNames[3]]->Fill(tsReader.
res_THETA_smt->at(
i), 1);
520 res_smt[vlID + paramNames[4]]->Fill(tsReader.
res_QOP_smt->at(
i), 1);
521 res_smt[vlID + paramNames[5]]->Fill(tsReader.
res_T_smt->at(
i), 1);
522 pull_smt[vlID + paramNames[0]]->Fill(tsReader.
pull_LOC0_smt->at(
i),
524 pull_smt[vlID + paramNames[1]]->Fill(tsReader.
pull_LOC1_smt->at(
i),
526 pull_smt[vlID + paramNames[2]]->Fill(tsReader.
pull_PHI_smt->at(
i), 1);
529 pull_smt[vlID + paramNames[4]]->Fill(tsReader.
pull_QOP_smt->at(
i), 1);
530 pull_smt[vlID + paramNames[5]]->Fill(tsReader.
pull_T_smt->at(
i), 1);
539 auto respull_mean_prf =
540 new TCanvas(
"respull_mean_prf",
541 "Residual/Pull Distributions - mean profiles", 1800, 800);
542 respull_mean_prf->Divide(3, 2);
544 auto respull_var_prf =
545 new TCanvas(
"respull_var_prf",
546 "Residual/Pull Distributions - variance profiles", 1800, 800);
547 respull_var_prf->Divide(3, 2);
549 auto plotProfiles = [&](std::array<TProfile2D*, 6>& profiles,
553 for (
size_t ipar = 0; ipar < paramNames.size(); ++ipar) {
554 respull_mean_prf->cd(ipar + 1);
556 if (res_pull ==
"pull") {
557 profiles[ipar]->GetZaxis()->SetRangeUser(-1. * pullRange, pullRange);
561 profiles[ipar]->Draw(
"colz");
562 profiles[ipar]->Write();
565 if (not saveAs.empty()) {
566 respull_mean_prf->SaveAs((
std::string(
"all_") + res_pull +
573 for (
size_t ipar = 0; ipar < paramNames.size(); ++ipar) {
574 respull_var_prf->cd(ipar + 1);
575 auto zAxis = profiles[ipar]->GetXaxis();
576 auto rAxis = profiles[ipar]->GetYaxis();
577 auto binsZ = zAxis->GetNbins();
578 auto binsR = rAxis->GetNbins();
580 hist_name += res_pull;
582 hist_name += paramNames[ipar];
586 new TH2F(hist_name.c_str(),
587 (res_pull +
std::string(
" ") + paramNames[ipar]).c_str(),
588 binsZ, zAxis->GetXmin(), zAxis->GetXmax(), binsR,
589 rAxis->GetXmin(), rAxis->GetXmax());
590 for (
int iz = 1; iz <= binsZ; ++iz) {
591 for (
int ir = 1; ir <= binsR; ++ir) {
592 if (profiles[ipar]->GetBinContent(iz, ir) != 0.) {
593 var->SetBinContent(iz, ir, profiles[ipar]->GetBinError(iz, ir));
597 if (res_pull ==
"pull") {
598 var->GetZaxis()->SetRangeUser(0., pullRange);
605 if (not saveAs.empty()) {
606 respull_var_prf->SaveAs((
std::string(
"all_") + res_pull +
615 plotProfiles(p2d_res_zr_prt,
"res",
"prt");
616 plotProfiles(p2d_pull_zr_prt,
"pull",
"prt");
619 plotProfiles(p2d_res_zr_flt,
"res",
"flt");
620 plotProfiles(p2d_pull_zr_flt,
"pull",
"flt");
623 plotProfiles(p2d_res_zr_smt,
"res",
"smt");
624 plotProfiles(p2d_pull_zr_smt,
"pull",
"smt");
629 new TCanvas(
"residuals",
"Residual Distributions", 1200, 800);
630 residuals->Divide(3, 2);
632 auto pulls =
new TCanvas(
"pulls",
"Pull distributions", 1200, 800);
635 for (
auto [vol, layers] : volLayIds) {
636 for (
auto lay : layers) {
637 auto vlID = volLayIdCut(vol, lay)[0];
640 for (
size_t ipar = 0; ipar < paramNames.size(); ipar++) {
641 residuals->cd(ipar + 1);
643 auto legend =
new TLegend(0.7, 0.7, 0.9, 0.9);
649 res_smt[
name]->DrawNormalized(
"");
650 res_smt[
name]->Write();
651 legend->AddEntry(res_smt[name],
"smoothed",
"l");
655 res_flt[
name]->DrawNormalized(drawOptions.c_str());
656 res_flt[
name]->Write();
657 legend->AddEntry(res_flt[name],
"filtered",
"l");
661 res_prt[
name]->DrawNormalized(drawOptions.c_str());
662 res_prt[
name]->Write();
663 legend->AddEntry(res_prt[name],
"predicted",
"l");
666 legend->SetBorderSize(0);
667 legend->SetFillStyle(0);
668 legend->SetTextFont(42);
671 if (not saveAs.empty()) {
672 residuals->SaveAs((vlID +
std::string(
"residuals.") + saveAs).c_str());
676 for (
size_t ipar = 0; ipar < paramNames.size(); ipar++) {
681 auto legend =
new TLegend(0.7, 0.5, 0.9, 0.9);
687 auto scale = 1. / pull_smt[
name]->Integral(
"width");
688 pull_smt[
name]->Scale(scale);
690 pull_smt[
name]->Write();
692 legend->AddEntry(pull_smt[name],
"smoothed",
"pe");
695 pull_smt[
name]->Fit(
"gaus",
"q");
696 TF1*
gauss = pull_smt[
name]->GetFunction(
"gaus");
697 gauss->SetLineColorAlpha(kBlack, 0.5);
699 auto mu = gauss->GetParameter(1);
700 auto sigma = gauss->GetParameter(2);
701 auto mu_fit_info = TString::Format(
"#mu = %.3f",
mu);
702 auto su_fit_info = TString::Format(
"#sigma = %.3f",
sigma);
704 legend->AddEntry(gauss, mu_fit_info.Data(),
"l");
705 legend->AddEntry(pull_prt[name], su_fit_info.Data(),
"");
712 auto scale = 1. / pull_flt[
name]->Integral(
"width");
713 pull_flt[
name]->Scale(scale);
715 pull_flt[
name]->Write();
717 legend->AddEntry(pull_flt[name],
"filtered",
"pe");
720 pull_flt[
name]->Fit(
"gaus",
"q",
"same");
721 TF1*
gauss = pull_flt[
name]->GetFunction(
"gaus");
722 gauss->SetLineColorAlpha(kBlue, 0.5);
724 auto mu = gauss->GetParameter(1);
725 auto sigma = gauss->GetParameter(2);
726 auto mu_fit_info = TString::Format(
"#mu = %.3f",
mu);
727 auto su_fit_info = TString::Format(
"#sigma = %.3f",
sigma);
729 legend->AddEntry(gauss, mu_fit_info.Data(),
"l");
730 legend->AddEntry(pull_prt[name], su_fit_info.Data(),
"");
735 auto drawOptions = (smoothed or filtered) ?
"pe same" :
"pe";
737 auto scale = 1. / pull_prt[
name]->Integral(
"width");
738 pull_prt[
name]->Scale(scale);
740 pull_prt[
name]->Write();
742 legend->AddEntry(pull_prt[name],
"predicted",
"pe");
745 pull_prt[
name]->Fit(
"gaus",
"q",
"same");
746 TF1*
gauss = pull_prt[
name]->GetFunction(
"gaus");
747 gauss->SetLineColorAlpha(kRed, 0.5);
749 auto mu = gauss->GetParameter(1);
750 auto sigma = gauss->GetParameter(2);
751 auto mu_fit_info = TString::Format(
"#mu = %.3f",
mu);
752 auto su_fit_info = TString::Format(
"#sigma = %.3f",
sigma);
754 legend->AddEntry(gauss, mu_fit_info.Data(),
"l");
755 legend->AddEntry(pull_prt[name], su_fit_info.Data(),
"");
760 auto ref =
new TF1(
"ref",
"ROOT::Math::normal_pdf(x,1,0)", -5, 5);
761 ref->SetLineColor(kGreen);
763 legend->AddEntry(
ref,
"#mu = 0 #sigma = 1",
"l");
765 legend->SetBorderSize(0);
766 legend->SetFillStyle(0);
767 legend->SetTextFont(42);
772 if (not saveAs.empty()) {
773 pulls->SaveAs((vlID +
std::string(
"pulls.") + saveAs).c_str());
778 if (out !=
nullptr) {