Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
boundParamResolution.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file boundParamResolution.C
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2019-2022 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
9 #include <array>
10 #include <iostream>
11 #include <map>
12 #include <string>
13 #include <vector>
14 
15 #include <TCanvas.h>
16 #include <TColor.h>
17 #include <TDirectory.h>
18 #include <TError.h>
19 #include <TF1.h>
20 #include <TFile.h>
21 #include <TH1F.h>
22 #include <TH2F.h>
23 #include <TLegend.h>
24 #include <TMath.h>
25 #include <TProfile2D.h>
26 #include <TStyle.h>
27 #include <TTree.h>
28 
29 #include "CommonUtils.h"
30 #include "TreeReader.h"
31 
32 using namespace ROOT;
33 
44 int boundParamResolution(const std::string& inFile, const std::string& treeName,
45  const std::string& outFile, bool predicted = true,
46  bool filtered = true, bool smoothed = true,
47  bool fitFiltered = true, bool fitPredicted = true,
48  bool fitSmoothed = true,
49  const std::string& saveAs = "") {
50  // Some style options
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);
57 
58  // Pull range, residual ranges are automatically determined
59  double pullRange = 5;
60 
61  //
62  // Residual ranges should be largest in predicted and smallest in smoothed,
63  // hence reverse the order.
64  //
65  std::string range_tag = "smt";
66  if (predicted) {
67  range_tag = "prt";
68  } else if (filtered) {
69  range_tag = "flt";
70  }
71 
72  // Section 0: file handling ---------------------------------------------
73  //
74  // Open root file written by RootTrajectoryWriter
75  // Create output root file
76  std::cout << "Opening file: " << inFile << std::endl;
77  TFile* file = TFile::Open(inFile.c_str(), "read");
78 
79  // Bail out if no tree was found
80  if (file == nullptr) {
81  return -1;
82  }
83 
84  std::cout << "Reading tree: " << treeName << std::endl;
85  TTree* tree = (TTree*)file->Get(treeName.c_str());
86 
87  // Bail out if no tree was found
88  if (tree == nullptr) {
89  return -2;
90  }
91 
92  TFile* out = TFile::Open(outFile.c_str(), "recreate");
93  out->cd();
94 
95  // Section 1: geometry parsing ------------------------------------------
96  //
97  // Gathers accessible volume_id and layer_id values
98  // Draw the volume_id : layer_id correlation matrix
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);
103  // Volume/Layer Id
104  geometryCanvas->cd(1);
105  tree->Draw("layer_id:volume_id>>volID_layID(100,0.5,100.5,100,0.5,100.5)", "",
106  "colz");
107  auto h2_volID_layID = dynamic_cast<TH2F*>(out->Get("volID_layID"));
108  setHistStyle(h2_volID_layID, 2);
109  h2_volID_layID->Write();
110  // Geometry size
111  geometryCanvas->cd(2);
112  std::string rz_draw_string = "(sqrt(g_x_";
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"));
125  setHistStyle(h2_geo_dim, 2);
126  float detectorZ = 1.15 * h2_geo_dim->GetXaxis()->GetXmax();
127  float detectorR = 1.15 * h2_geo_dim->GetYaxis()->GetXmax();
128  h2_geo_dim->Write();
129 
130  // Save the plots on screen
131  if (not saveAs.empty()) {
132  geometryCanvas->SaveAs((std::string("all_vol_lay_ids.") + saveAs).c_str());
133  }
134 
135  // Now determine the valid volume / layer bins
136  int volBins = h2_volID_layID->GetXaxis()->GetNbins();
137  int layBins = h2_volID_layID->GetYaxis()->GetNbins();
138  // Add the overall plots
139  volLayIds[-1] = {-1};
140  // Go through volumes and add plots per volume
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()) {
145  // First occurrence of this layer, add -1 for volume plots
146  volLayIds[volID] = {-1, layID};
147  } else {
148  volLayIds[volID].push_back(layID);
149  }
150  }
151  }
152  }
153 
154  // Section 2: Branch assignment -----------------------------------------
155  //
156  // Helper for assigning branches to the input tree
157  TrackStatesReader tsReader(tree, false);
158 
159  // Section 3: Histogram booking -----------------------------------------
160 
161  TCanvas* rangeCanvas =
162  new TCanvas("rangeCanvas", "Range Estimation", 10, 10, 900, 600);
163  rangeCanvas->Divide(3, 2);
164 
165  std::vector<std::string> res_ranges = {std::string("res_eLOC0_") + range_tag,
166  std::string("res_eLOC1_") + range_tag,
167  std::string("res_ePHI_") + range_tag,
168  std::string("res_eTHETA_") + range_tag,
169  std::string("res_eQOP_") + range_tag,
170  std::string("res_eT_") + range_tag};
171 
175  auto volLayIdCut = [](int vol, int lay) -> std::array<std::string, 2> {
176  if (vol < 0 and lay < 0) {
177  return {std::string("all_"), ""};
178  }
179 
180  std::string vlstr = "vol" + std::to_string(vol);
181  std::string vlcut = "volume_id == " + std::to_string(vol);
182  if (lay > 0) {
183  vlstr += (std::string("_lay") + std::to_string(lay));
184  vlcut += (std::string(" && layer_id == ") + std::to_string(lay));
185  }
186  vlstr += std::string("_");
187  return {vlstr, vlcut};
188  };
189 
193  auto histRanges =
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);
198  std::string drawString = res_ranges[ir] + ">>";
199  std::string histString =
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()));
203  h1_range->Write();
204  float range = pullRange * h1_range->GetRMS();
205  ranges[ir] = range;
206  }
207  if (not saveAs.empty()) {
208  rangeCanvas->SaveAs(
209  (vlIdCut[0] + std::string("res_ranges.") + saveAs).c_str());
210  }
211  return ranges;
212  };
213 
214  // Track parameter name
215  std::vector<std::string> paramNames = {"loc0", "loc1", "#phi",
216  "#theta", "q/p", "t"};
217 
218  // Book histograms (with adapted range):
219  //
220  // Global profile histograms : residuals/pulls
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{};
227 
228  for (unsigned int ipar = 0; ipar < paramNames.size(); ++ipar) {
229  const auto& par = paramNames[ipar];
230 
231  if (predicted) {
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);
236 
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);
241 
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()));
246  setHistStyle(p2d_res_zr_prt[ipar], 1);
247 
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()));
253  setHistStyle(p2d_pull_zr_prt[ipar], 1);
254  }
255 
256  if (filtered) {
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);
265 
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()));
270  setHistStyle(p2d_res_zr_flt[ipar], 2);
271 
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()));
277  setHistStyle(p2d_pull_zr_flt[ipar], 2);
278  }
279 
280  if (smoothed) {
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);
285 
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);
290 
291  p2d_res_zr_smt[ipar]->SetErrorOption("s");
292  p2d_pull_zr_smt[ipar]->SetErrorOption("s");
293 
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()));
297  setHistStyle(p2d_res_zr_smt[ipar], 4);
298 
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()));
303  setHistStyle(p2d_pull_zr_smt[ipar], 4);
304  }
305  }
306 
307  // Residual / Pull histograms
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;
314 
315  // - per layer (identified by vol, lay)
316  // - per volume (identified by vol, -1)
317  // - overall (identified by lay)
318  for (auto [vol, layers] : volLayIds) {
319  for (auto lay : layers) {
320  // Estimate the ranges from smoothed
321  auto vlIdCut = volLayIdCut(vol, lay);
322  auto ranges = histRanges(vlIdCut);
323 
324  // Residual range
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]}};
329 
330  // Create the hists and set up for them
331  for (const auto& [partwin, resRange] : paramResidualRange) {
332  // histogram creation
333  std::string par = partwin.first;
334  std::string id_par = vlIdCut[0] + par;
335 
336  TString par_string(partwin.second.c_str());
337  TString res_string =
338  par_string + TString("^{rec} - ") + par_string + TString("^{true}");
339 
340  TString pull_string = TString("(") + res_string +
341  TString(")/#sigma_{") + par_string + TString("}");
342 
343  if (predicted) {
344  res_prt[id_par] =
345  new TH1F(Form((vlIdCut[0] + std::string("res_prt_%s")).c_str(),
346  par.c_str()),
347  Form("residual of %s", par.c_str()), 100, -1 * resRange,
348  resRange);
349  res_prt[id_par]->GetXaxis()->SetTitle(res_string.Data());
350  res_prt[id_par]->GetYaxis()->SetTitle("Entries");
351  setHistStyle(res_prt[id_par], kRed);
352 
353  pull_prt[id_par] = new TH1F(
354  Form((vlIdCut[0] + std::string("pull_prt_%s")).c_str(),
355  par.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");
359  setHistStyle(pull_prt[id_par], kRed);
360  }
361 
362  if (filtered) {
363  res_flt[id_par] =
364  new TH1F(Form((vlIdCut[0] + std::string("res_flt_%s")).c_str(),
365  par.c_str()),
366  Form("residual of %s", par.c_str()), 100, -1 * resRange,
367  resRange);
368  res_flt[id_par]->GetXaxis()->SetTitle(res_string.Data());
369  res_flt[id_par]->GetYaxis()->SetTitle("Entries");
370  setHistStyle(res_flt[id_par], kBlue);
371 
372  pull_flt[id_par] = new TH1F(
373  Form((vlIdCut[0] + std::string("pull_flt_%s")).c_str(),
374  par.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");
378  setHistStyle(pull_flt[id_par], kBlue);
379  }
380 
381  if (smoothed) {
382  res_smt[id_par] =
383  new TH1F(Form((vlIdCut[0] + std::string("res_smt_%s")).c_str(),
384  par.c_str()),
385  Form("residual of %s", par.c_str()), 100, -1 * resRange,
386  resRange);
387 
388  res_smt[id_par]->GetXaxis()->SetTitle(res_string.Data());
389  res_smt[id_par]->GetYaxis()->SetTitle("Entries");
390  setHistStyle(res_smt[id_par], kBlack);
391 
392  pull_smt[id_par] = new TH1F(
393  Form((vlIdCut[0] + std::string("pull_smt_%s")).c_str(),
394  par.c_str()),
395  Form("pull of %s", par.c_str()), 100, -1 * pullRange, pullRange);
396 
397  pull_smt[id_par]->GetXaxis()->SetTitle(pull_string.Data());
398  pull_smt[id_par]->GetYaxis()->SetTitle("Arb. Units");
399  setHistStyle(pull_smt[id_par], kBlack);
400  }
401  }
402  }
403  }
404 
405  // Section 4: Histogram filling -----------------------------------------
406  //
407  // - Running through the entries and filling the histograms
408  int entries = tree->GetEntries();
409  for (int j = 0; j < entries; j++) {
410  tsReader.getEntry(j);
411 
412  for (unsigned int i = 0; i < tsReader.nMeasurements; i++) {
413  // global profile filling
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));
431  }
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));
449  }
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));
467  }
468 
469  int vol = tsReader.volume_id->at(i);
470  int lay = tsReader.layer_id->at(i);
471 
473  std::vector<std::array<int, 2>> fillIds = {
474  {-1, -1}, {vol, -1}, {vol, lay}};
475 
476  for (const auto& fid : fillIds) {
477  auto vlID = volLayIdCut(fid[0], fid[1])[0];
478  // Fill predicated parameters
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),
487  1);
488  pull_prt[vlID + paramNames[1]]->Fill(tsReader.pull_LOC1_prt->at(i),
489  1);
490  pull_prt[vlID + paramNames[2]]->Fill(tsReader.pull_PHI_prt->at(i), 1);
491  pull_prt[vlID + paramNames[3]]->Fill(tsReader.pull_THETA_prt->at(i),
492  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);
495  }
496  // Fill filtered parameters
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),
505  1);
506  pull_flt[vlID + paramNames[1]]->Fill(tsReader.pull_LOC1_flt->at(i),
507  1);
508  pull_flt[vlID + paramNames[2]]->Fill(tsReader.pull_PHI_flt->at(i), 1);
509  pull_flt[vlID + paramNames[3]]->Fill(tsReader.pull_THETA_flt->at(i),
510  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);
513  }
514  // Fill smoothed parameters
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),
523  1);
524  pull_smt[vlID + paramNames[1]]->Fill(tsReader.pull_LOC1_smt->at(i),
525  1);
526  pull_smt[vlID + paramNames[2]]->Fill(tsReader.pull_PHI_smt->at(i), 1);
527  pull_smt[vlID + paramNames[3]]->Fill(tsReader.pull_THETA_smt->at(i),
528  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);
531  }
532  }
533  }
534  }
535 
536  // Section 5: Histogram plotting
537 
538  // Plotting global profiles
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);
543 
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);
548 
549  auto plotProfiles = [&](std::array<TProfile2D*, 6>& profiles,
550  const std::string& res_pull,
551  const std::string& type) -> void {
552  // Mean
553  for (size_t ipar = 0; ipar < paramNames.size(); ++ipar) {
554  respull_mean_prf->cd(ipar + 1);
555 
556  if (res_pull == "pull") {
557  profiles[ipar]->GetZaxis()->SetRangeUser(-1. * pullRange, pullRange);
558  }
559  adaptColorPalette(profiles[ipar], -1. * pullRange, pullRange, 0., 0.25,
560  104);
561  profiles[ipar]->Draw("colz");
562  profiles[ipar]->Write();
563  }
564  // Save the canvas: mean
565  if (not saveAs.empty()) {
566  respull_mean_prf->SaveAs((std::string("all_") + res_pull +
567  std::string("_mean_prf_") + type +
568  std::string(".") + saveAs)
569  .c_str());
570  }
571 
572  // Variance
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();
579  std::string hist_name = "all_";
580  hist_name += res_pull;
581  hist_name += "_";
582  hist_name += paramNames[ipar];
583  hist_name += "_";
584  hist_name += type;
585  TH2F* var =
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));
594  }
595  }
596  }
597  if (res_pull == "pull") {
598  var->GetZaxis()->SetRangeUser(0., pullRange);
599  }
600  adaptColorPalette(var, 0., pullRange, 1., 0.5, 104);
601  var->Draw("colz");
602  var->Write();
603  }
604  // Save the canvas: pulls
605  if (not saveAs.empty()) {
606  respull_var_prf->SaveAs((std::string("all_") + res_pull +
607  std::string("_var_prf_") + type +
608  std::string(".") + saveAs)
609  .c_str());
610  }
611  };
612 
613  // Plotting profiles: res/pulls
614  if (predicted) {
615  plotProfiles(p2d_res_zr_prt, "res", "prt");
616  plotProfiles(p2d_pull_zr_prt, "pull", "prt");
617  }
618  if (filtered) {
619  plotProfiles(p2d_res_zr_flt, "res", "flt");
620  plotProfiles(p2d_pull_zr_flt, "pull", "flt");
621  }
622  if (smoothed) {
623  plotProfiles(p2d_res_zr_smt, "res", "smt");
624  plotProfiles(p2d_pull_zr_smt, "pull", "smt");
625  }
626 
627  // Plotting residual
628  auto residuals =
629  new TCanvas("residuals", "Residual Distributions", 1200, 800);
630  residuals->Divide(3, 2);
631 
632  auto pulls = new TCanvas("pulls", "Pull distributions", 1200, 800);
633  pulls->Divide(3, 2);
634 
635  for (auto [vol, layers] : volLayIds) {
636  for (auto lay : layers) {
637  auto vlID = volLayIdCut(vol, lay)[0];
638 
639  // Residual plotting
640  for (size_t ipar = 0; ipar < paramNames.size(); ipar++) {
641  residuals->cd(ipar + 1);
642 
643  auto legend = new TLegend(0.7, 0.7, 0.9, 0.9);
644 
645  const std::string name = vlID + paramNames.at(ipar);
646 
647  // Draw them
648  if (smoothed) {
649  res_smt[name]->DrawNormalized("");
650  res_smt[name]->Write();
651  legend->AddEntry(res_smt[name], "smoothed", "l");
652  }
653  if (filtered) {
654  std::string drawOptions = smoothed ? "same" : "";
655  res_flt[name]->DrawNormalized(drawOptions.c_str());
656  res_flt[name]->Write();
657  legend->AddEntry(res_flt[name], "filtered", "l");
658  }
659  if (predicted) {
660  std::string drawOptions = (smoothed or filtered) ? "same" : "";
661  res_prt[name]->DrawNormalized(drawOptions.c_str());
662  res_prt[name]->Write();
663  legend->AddEntry(res_prt[name], "predicted", "l");
664  }
665 
666  legend->SetBorderSize(0);
667  legend->SetFillStyle(0);
668  legend->SetTextFont(42);
669  legend->Draw();
670  }
671  if (not saveAs.empty()) {
672  residuals->SaveAs((vlID + std::string("residuals.") + saveAs).c_str());
673  }
674 
675  // Pull plotting & writing
676  for (size_t ipar = 0; ipar < paramNames.size(); ipar++) {
677  const std::string name = vlID + paramNames.at(ipar);
678 
679  pulls->cd(ipar + 1);
680 
681  auto legend = new TLegend(0.7, 0.5, 0.9, 0.9);
682 
683  // Fit the smoothed distribution
684  if (smoothed) {
685  auto drawOptions = "pe";
686 
687  auto scale = 1. / pull_smt[name]->Integral("width");
688  pull_smt[name]->Scale(scale);
689  pull_smt[name]->Draw(drawOptions);
690  pull_smt[name]->Write();
691 
692  legend->AddEntry(pull_smt[name], "smoothed", "pe");
693 
694  if (fitSmoothed) {
695  pull_smt[name]->Fit("gaus", "q");
696  TF1* gauss = pull_smt[name]->GetFunction("gaus");
697  gauss->SetLineColorAlpha(kBlack, 0.5);
698 
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);
703 
704  legend->AddEntry(gauss, mu_fit_info.Data(), "l");
705  legend->AddEntry(pull_prt[name], su_fit_info.Data(), "");
706  }
707  }
708 
709  if (filtered) {
710  auto drawOptions = smoothed ? "pe same" : "pe";
711 
712  auto scale = 1. / pull_flt[name]->Integral("width");
713  pull_flt[name]->Scale(scale);
714  pull_flt[name]->Draw(drawOptions);
715  pull_flt[name]->Write();
716 
717  legend->AddEntry(pull_flt[name], "filtered", "pe");
718 
719  if (fitFiltered) {
720  pull_flt[name]->Fit("gaus", "q", "same");
721  TF1* gauss = pull_flt[name]->GetFunction("gaus");
722  gauss->SetLineColorAlpha(kBlue, 0.5);
723 
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);
728 
729  legend->AddEntry(gauss, mu_fit_info.Data(), "l");
730  legend->AddEntry(pull_prt[name], su_fit_info.Data(), "");
731  }
732  }
733 
734  if (predicted) {
735  auto drawOptions = (smoothed or filtered) ? "pe same" : "pe";
736 
737  auto scale = 1. / pull_prt[name]->Integral("width");
738  pull_prt[name]->Scale(scale);
739  pull_prt[name]->Draw(drawOptions);
740  pull_prt[name]->Write();
741 
742  legend->AddEntry(pull_prt[name], "predicted", "pe");
743 
744  if (fitPredicted) {
745  pull_prt[name]->Fit("gaus", "q", "same");
746  TF1* gauss = pull_prt[name]->GetFunction("gaus");
747  gauss->SetLineColorAlpha(kRed, 0.5);
748 
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);
753 
754  legend->AddEntry(gauss, mu_fit_info.Data(), "l");
755  legend->AddEntry(pull_prt[name], su_fit_info.Data(), "");
756  }
757  }
758 
759  // Reference standard normal pdf
760  auto ref = new TF1("ref", "ROOT::Math::normal_pdf(x,1,0)", -5, 5);
761  ref->SetLineColor(kGreen);
762  ref->Draw("same");
763  legend->AddEntry(ref, "#mu = 0 #sigma = 1", "l");
764 
765  legend->SetBorderSize(0);
766  legend->SetFillStyle(0);
767  legend->SetTextFont(42);
768  legend->Draw();
769  }
770 
771  // Save the Canvases as pictures
772  if (not saveAs.empty()) {
773  pulls->SaveAs((vlID + std::string("pulls.") + saveAs).c_str());
774  }
775  }
776  }
777 
778  if (out != nullptr) {
779  out->Close();
780  }
781  return 0;
782 }