Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
trackSummaryAnalysis.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file trackSummaryAnalysis.C
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2019-2021 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 <bitset>
10 #include <fstream>
11 #include <iostream>
12 #include <limits>
13 #include <map>
14 #include <string>
15 #include <vector>
16 
17 #include <TCanvas.h>
18 #include <TChain.h>
19 #include <TF1.h>
20 #include <TFile.h>
21 #include <TH1F.h>
22 #include <TH2F.h>
23 #include <TH3F.h>
24 #include <TLegend.h>
25 #include <TMath.h>
26 #include <TStyle.h>
27 #include <TTree.h>
28 #include <TVectorF.h>
29 
30 #include "CommonUtils.h"
31 #include "TreeReader.h"
32 
33 using namespace ROOT;
34 
53 int trackSummaryAnalysis(
54  const std::vector<std::string>& inFiles, const std::string& inTree,
55  const std::string& outFile, const std::string& inConfig = "",
56  const std::string& outConfig = "", unsigned long nEntries = 0,
57  unsigned int nPeakEntries = 0, float pullRange = 6.,
58  unsigned int nHistBins = 61, unsigned int nPhiBins = 10,
59  const std::array<float, 2>& phiRange = {-M_PI, M_PI},
60  unsigned int nEtaBins = 10, const std::array<float, 2>& etaRange = {-3, 3},
61  const std::vector<double>& ptBorders =
62  {0., std::numeric_limits<double>::infinity()},
63  const std::bitset<7>& residualPulls = std::bitset<7>{"1111111"},
64  const std::bitset<5>& auxiliary = std::bitset<5>{"11111"}) {
65  // Load the tree chain
66  TChain* treeChain = new TChain(inTree.c_str());
67  for (const auto& inFile : inFiles) {
68  treeChain->Add(inFile.c_str());
69  // Open root file written by RootTrajectoryWriter
70  std::cout << "*** Adding file: " << inFile << std::endl;
71  }
72 
73  if (treeChain->GetEntries() == 0) {
74  std::cout << "xxx No entries found ... bailing out." << std::endl;
75  return -1;
76  }
77 
78  TCanvas* rangeCanvas =
79  new TCanvas("rangeCanvas", "Range Estimation", 100, 100, 620, 400);
80 
81  TrackSummaryReader tracks(treeChain, true);
82 
83  // Loop over the entries of the keys
84  unsigned long entries = estimateEntries(*tracks.tree, nEntries);
85  unsigned long peakEntries = estimateEntries(*tracks.tree, nPeakEntries);
86 
87  // Temporary counter
88  unsigned int histBarcode = 0;
89 
90  // Deduction
91  unsigned int nPtBins = ptBorders.size() - 1;
92 
93  // One time initialization of the residual and pull handles
94  using ResidualPulls = std::vector<ResidualPullHandle>;
95  ResidualPulls baseResidualPulls = {};
96 
97  if (residualPulls.test(0)) {
98  // A standard d0Handle
99  ResidualPullHandle d0Handle;
100  d0Handle.tag = "d0";
101  d0Handle.residualStr = "d_{0}^{rec} - d_{0}^{true}";
102  d0Handle.residualUnit = "[mm]";
103  d0Handle.errorStr = "#sigma(d_{0})";
104  d0Handle.rangeDrawStr = "eLOC0_fit-t_d0";
105  d0Handle.rangeMaxStr = "(1000,-10,10)";
106  d0Handle.value = ResidualAccessor{tracks.eLOC0_fit, tracks.t_d0};
107  d0Handle.error = DirectAccessor<float>{tracks.err_eLOC0_fit};
108  d0Handle.accept = AcceptAll{};
109  // Push it
110  baseResidualPulls.push_back(d0Handle);
111  }
112 
113  if (residualPulls.test(1)) {
114  // A standard z0Handle
115  ResidualPullHandle z0Handle;
116  z0Handle.tag = "z0";
117  z0Handle.residualStr = "z_{0}^{rec} - z_{0}^{true}";
118  z0Handle.residualUnit = "[mm]";
119  z0Handle.errorStr = "#sigma(z_{0})";
120  z0Handle.rangeDrawStr = "eLOC1_fit-t_z0";
121  z0Handle.rangeMaxStr = "(1000,-10,10)";
122  z0Handle.value = ResidualAccessor{tracks.eLOC1_fit, tracks.t_z0};
123  z0Handle.error = DirectAccessor<float>{tracks.err_eLOC1_fit};
124  z0Handle.accept = AcceptAll{};
125  // Push it
126  baseResidualPulls.push_back(z0Handle);
127  }
128 
129  if (residualPulls.test(2)) {
130  // A standard phi0Handle
131  ResidualPullHandle phi0Handle;
132  phi0Handle.tag = "phi0";
133  phi0Handle.residualStr = "#phi_{0}^{rec} - #phi_{0}^{true}";
134  phi0Handle.errorStr = "#sigma(phi_{0})";
135  phi0Handle.rangeDrawStr = "ePHI_fit-t_phi";
136  phi0Handle.rangeMaxStr = "(1000,-0.01,0.01)";
137  phi0Handle.value = ResidualAccessor{tracks.ePHI_fit, tracks.t_phi};
138  phi0Handle.error = DirectAccessor<float>{tracks.err_ePHI_fit};
139  phi0Handle.accept = AcceptAll{};
140  // Push it
141  baseResidualPulls.push_back(phi0Handle);
142  }
143 
144  if (residualPulls.test(3)) {
145  // A standard theta0Handle
146  ResidualPullHandle theta0Handle;
147  theta0Handle.tag = "theta0";
148  theta0Handle.residualStr = "#theta_{0}^{rec} - #theta_{0}^{true}";
149  theta0Handle.errorStr = "#sigma(theta_{0})";
150  theta0Handle.rangeDrawStr = "eTHETA_fit-t_theta";
151  theta0Handle.rangeMaxStr = "(1000,-0.01,0.01)";
152  theta0Handle.value = ResidualAccessor{tracks.eTHETA_fit, tracks.t_theta};
153  theta0Handle.error = DirectAccessor<float>{tracks.err_eTHETA_fit};
154  theta0Handle.accept = AcceptAll{};
155  // Push it
156  baseResidualPulls.push_back(theta0Handle);
157  }
158 
159  if (residualPulls.test(4)) {
160  // The standard qop Handle
161  ResidualPullHandle qopHandle;
162  qopHandle.tag = "qop";
163  qopHandle.residualStr = "q/p^{rec} - q/p^{true}";
164  qopHandle.residualUnit = "[GeV^{-1}]";
165  qopHandle.errorStr = "#sigma(q/p)";
166  qopHandle.rangeDrawStr = "eQOP_fit-t_charge/t_p";
167  qopHandle.rangeMaxStr = "(1000,-0.1,0.1)";
168  qopHandle.value =
169  QopResidualAccessor{tracks.eQOP_fit, tracks.t_charge, tracks.t_p};
170  qopHandle.error = DirectAccessor<float>{tracks.err_eQOP_fit};
171  qopHandle.accept = AcceptAll{};
172  // Push it
173  baseResidualPulls.push_back(qopHandle);
174  }
175 
176  if (residualPulls.test(5)) {
177  // The pt measurement
178  ResidualPullHandle tHandle;
179  tHandle.tag = "t";
180  tHandle.residualStr = "t^{rec} - t^{true}";
181  tHandle.residualUnit = "[ns]";
182  tHandle.errorStr = "#sigma(t})";
183  tHandle.rangeDrawStr = "eT_fit - t_time";
184  tHandle.rangeMaxStr = "(1000,-10.,10.)";
185  tHandle.value = ResidualAccessor{tracks.eT_fit, tracks.t_time};
186  tHandle.error = DirectAccessor<float>{tracks.err_eT_fit};
187  tHandle.accept = AcceptAll{};
188  // Push it
189  baseResidualPulls.push_back(tHandle);
190  }
191 
192  if (residualPulls.test(6)) {
193  // The pt measurement
194  ResidualPullHandle ptHandle;
195  ptHandle.tag = "pt";
196  ptHandle.residualStr = "p_{T}^{rec} - p_{T}^{true}";
197  ptHandle.residualUnit = "[GeV]";
198  ptHandle.errorStr = "#sigma(p_{T})";
199  ptHandle.rangeDrawStr = "1./abs(eQOP_fit) * sin(eTHETA_fit) - (1./t_pT)";
200  ptHandle.rangeMaxStr = "(1000,-10.,10.)";
201  ptHandle.value =
202  PtResidualAccessor{tracks.eQOP_fit, tracks.eTHETA_fit, tracks.t_pT};
203  ptHandle.error = PtErrorAccessor{tracks.eQOP_fit, tracks.err_eQOP_fit,
204  tracks.eTHETA_fit, tracks.err_eTHETA_fit};
205  ptHandle.accept = AcceptAll{};
206  // Push it
207  baseResidualPulls.push_back(ptHandle);
208  }
209 
210  using Auxiliaries = std::vector<SingleHandle>;
211  Auxiliaries baseAuxilaries;
212 
213  if (auxiliary.test(0)) {
214  // Chi2/ndf
215  SingleHandle chi2ndf;
216  chi2ndf.tag = "chi2ndf";
217  chi2ndf.label = "#Chi^{2}/ndf";
218  chi2ndf.bins = nHistBins;
219  chi2ndf.range = {0., 5.};
220  chi2ndf.value =
221  DivisionAccessor<float, unsigned int>{tracks.chi2Sum, tracks.NDF};
222  chi2ndf.accept = AcceptAll{};
223  // Push it
224  baseAuxilaries.push_back(chi2ndf);
225  }
226 
227  if (auxiliary.test(1)) {
228  // Measurements
229  SingleHandle measurements;
230  measurements.tag = "measurements";
231  measurements.label = "#(measurements)";
232  measurements.rangeDrawStr = "nMeasurements";
233  measurements.value = DirectAccessor<unsigned int>{tracks.nMeasurements};
234  measurements.accept = AcceptAll{};
235  estimateIntegerRange(measurements, *rangeCanvas, *tracks.tree, peakEntries,
236  250, 5, ++histBarcode);
237  // Push it
238  baseAuxilaries.push_back(measurements);
239  }
240 
241  if (auxiliary.test(2)) {
242  // Holes
243  SingleHandle holes;
244  holes.tag = "holes";
245  holes.label = "#(holes)";
246  holes.rangeDrawStr = "nHoles";
247  holes.value = DirectAccessor<unsigned int>{tracks.nHoles};
248  holes.accept = AcceptAll{};
249  estimateIntegerRange(holes, *rangeCanvas, *tracks.tree, peakEntries, 10, 2,
250  ++histBarcode);
251  // Push it
252  baseAuxilaries.push_back(holes);
253  }
254 
255  if (auxiliary.test(3)) {
256  // Holes
257  SingleHandle outliers;
258  outliers.tag = "outliers";
259  outliers.label = "#(outliers)";
260  outliers.rangeDrawStr = "nOutliers";
261  outliers.value = DirectAccessor<unsigned int>{tracks.nOutliers};
262  outliers.accept = AcceptAll{};
263  estimateIntegerRange(outliers, *rangeCanvas, *tracks.tree, peakEntries, 10,
264  2, ++histBarcode);
265  // Push it
266  baseAuxilaries.push_back(outliers);
267  }
268 
269  if (auxiliary.test(4)) {
270  // Holes
271  SingleHandle shared;
272  shared.tag = "shared";
273  shared.label = "#(shared)";
274  shared.rangeDrawStr = "nSharedHits";
275  shared.value = DirectAccessor<unsigned int>{tracks.nSharedHits};
276  shared.accept = AcceptAll{};
277  estimateIntegerRange(shared, *rangeCanvas, *tracks.tree, peakEntries, 10, 2,
278  ++histBarcode);
279  // Push it
280  baseAuxilaries.push_back(shared);
281  }
282 
283  // Preparation phase for handles
284 #ifdef NLOHMANN_AVAILABLE
285  nlohmann::json handle_configs;
286  if (not inConfig.empty()) {
287  std::ifstream ifs(inConfig.c_str());
288  handle_configs = nlohmann::json::parse(ifs);
289  }
290 
299  auto handleRange = [&](ResidualPullHandle& handle, const TString& handleTag, unsigned int peakE) -> void {
300  bool rangeDetermined = false;
301  if (not inConfig.empty()) {
302  if (handle_configs.contains((handleTag).Data())) {
303  auto handle_config = handle_configs[(handleTag).Data()];
304  handle.range = handle_config["range"].get<std::array<float, 2>>();
305  rangeDetermined = true;
306  }
307  }
308  if (not rangeDetermined) {
309  estimateResiudalRange(handle, *rangeCanvas, *tracks.tree, peakE,
310  ++histBarcode);
311  }
312 
313  if (not outConfig.empty()) {
314  nlohmann::json range_config;
315  range_config["range"] = handle.range;
316  handle_configs[(handleTag).Data()] = range_config;
317  }
318  };
319 #else
320 
321 
322 
323 
324 
325  auto handleRange = [&](ResidualPullHandle& handle, const TString& handleTag,
326  unsigned long peakEntries) -> void {
327  estimateResiudalRange(handle, *rangeCanvas, *tracks.tree, peakEntries,
328  ++histBarcode);
329  };
330 #endif
331 
332  // Full Range handles - they accept all tracks
333  ResidualPulls fullResidualPulls = baseResidualPulls;
334  // Range Estimation and booking histogram
335  for (auto& fHandle : fullResidualPulls) {
336  // The full tag
337  TString handleTag(fHandle.tag + std::string("_all"));
338  handleRange(fHandle, handleTag, peakEntries);
339  // Book histograms
340  bookHistograms(fHandle, pullRange, nHistBins, ++histBarcode);
341  // The Histogram names
342  TString residualN = TString("res_") + handleTag;
343  TString pullN = TString("pull_") + handleTag;
344  // Style and name
345  setHistStyle(fHandle.residualHist);
346  fHandle.residualHist->SetName(residualN);
347  setHistStyle(fHandle.pullHist);
348  fHandle.pullHist->SetName(pullN);
349  }
350 
351  // Regional/Binned handles
352  using ResidualPullsVector = std::vector<ResidualPulls>;
353  using ResidualPullsMatrix = std::vector<ResidualPullsVector>;
354 
355  // Eta-Pt residual/pull handles
356  ResidualPullsVector ptResidualPulls =
357  ResidualPullsVector(nPtBins, baseResidualPulls);
358  ResidualPullsMatrix etaPtResidualPulls =
359  ResidualPullsMatrix(nEtaBins, ptResidualPulls);
360 
361  // Eta-Phi residual/pull handles
362  ResidualPullsVector phiResidualPulls =
363  ResidualPullsVector(nPhiBins, baseResidualPulls);
364  ResidualPullsMatrix etaPhiResidualPulls =
365  ResidualPullsMatrix(nEtaBins, phiResidualPulls);
366 
367  Auxiliaries fullAuxiliaries = baseAuxilaries;
368 
369  // Histogram booking for full range auxiliaries
370  for (auto& fAuxiliary : fullAuxiliaries) {
371  fAuxiliary.hist =
372  new TH1F(fAuxiliary.tag.c_str(), fAuxiliary.tag.c_str(),
373  fAuxiliary.bins, fAuxiliary.range[0], fAuxiliary.range[1]);
374  fAuxiliary.hist->GetXaxis()->SetTitle(fAuxiliary.label.c_str());
375  fAuxiliary.hist->GetYaxis()->SetTitle("Entries");
376  setHistStyle(fAuxiliary.hist);
377  }
378 
379  using AuxiliariesVector = std::vector<Auxiliaries>;
380  using AuxiliariesMatrix = std::vector<AuxiliariesVector>;
381 
382  // Eta-Pt auxiliaries
383  AuxiliariesVector ptAuxiliaries = AuxiliariesVector(nPtBins, baseAuxilaries);
384  AuxiliariesMatrix etaPtAuxiliaries =
385  AuxiliariesMatrix(nEtaBins, ptAuxiliaries);
386 
387  // Eta-Phi auxiliaries
388  AuxiliariesVector phiAuxiliaries =
389  AuxiliariesVector(nPhiBins, baseAuxilaries);
390  AuxiliariesMatrix etaPhiAuxiliaries =
391  AuxiliariesMatrix(nEtaBins, phiAuxiliaries);
392 
393  // Loop over the binned handles & fill the acceptors
394  float phiStep = (phiRange[1] - phiRange[0]) / nPhiBins;
395  float etaStep = (etaRange[1] - etaRange[0]) / nEtaBins;
396 
397  TVectorF phiVals(nPhiBins + 1);
398  TVectorF etaVals(nEtaBins + 1);
399  TVectorF ptVals(nPtBins + 1);
400 
401  phiVals[0] = phiRange[0];
402  etaVals[0] = etaRange[0];
403  ptVals[0] = ptBorders[0];
404 
405 #ifdef BOOST_AVAILABLE
406  std::cout << "*** Handle Preparation: " << std::endl;
407  progress_display handle_preparation_progress(
408  nPhiBins * nEtaBins * nPtBins * baseResidualPulls.size());
409 #endif
410 
411  std::string land = " && ";
412 
414  for (unsigned int iphi = 0; iphi < nPhiBins; ++iphi) {
415  // Prepare the phi range for this batch
416  float phiMin = phiRange[0] + iphi * phiStep;
417  float phiMax = phiRange[0] + (iphi + 1) * phiStep;
418  phiVals[iphi + 1] = phiMax;
419 
420  // Acceptance range
421  AcceptRange phiBin{tracks.t_phi, {phiMin, phiMax}};
422  // Name tag
423  TString phiTag = "_phi";
424  phiTag += iphi;
425  // Range cut string
426  std::string phiCut = "t_phi >= ";
427  phiCut += std::to_string(phiMin);
428  phiCut += land;
429  phiCut += std::string("t_phi < ");
430  phiCut += std::to_string(phiMax);
431 
432  for (unsigned int ieta = 0; ieta < nEtaBins; ++ieta) {
433  // Prepare the eta ragne for this batch
434  float etaMin = etaRange[0] + ieta * etaStep;
435  float etaMax = etaRange[0] + (ieta + 1) * etaStep;
436  etaVals[ieta + 1] = etaMax;
437  // Acceptance range
438  AcceptRange etaBin{tracks.t_eta, {etaMin, etaMax}};
439  // Name tag
440  TString etaTag = "_eta";
441  etaTag += ieta;
442  // Range cut string
443  std::string etaCut = "t_eta >= ";
444  etaCut += std::to_string(etaMin);
445  etaCut += land;
446  etaCut += std::string("t_eta < ");
447  etaCut += std::to_string(etaMax);
448 
449  // Combined eta/phi acceptance
450  AcceptCombination etaPhiBin{etaBin, phiBin};
451  TString etaPhiTag = etaTag + phiTag;
452 
453  for (unsigned int ipt = 0; ipt < nPtBins; ++ipt) {
454  // Acceptance range for this pT bin
455  float ptMin = static_cast<float>(ptBorders[ipt]);
456  float ptMax = static_cast<float>(ptBorders[ipt + 1]);
457  AcceptRange ptBin{tracks.t_pT, {ptMin, ptMax}};
458 
459  float upperPtBorder =
460  ptBorders[ipt + 1] == std::numeric_limits<double>::infinity()
461  ? 100.
462  : ptBorders[ipt + 1];
463  ptVals[ipt + 1] = upperPtBorder;
464  // Name tag
465  TString ptTag = "_pt";
466  ptTag += ipt;
467 
468  // Range cut string
469  std::string ptCut = "t_pT >= ";
470  ptCut += std::to_string(ptMin);
471  ptCut += land;
472  ptCut += std::string("t_pT < ");
473  ptCut += std::to_string(ptMax);
474 
475  // Combined eta/pt acceptance
476  AcceptCombination etaPtBin{etaBin, ptBin};
477 
478  for (unsigned int iresp = 0; iresp < baseResidualPulls.size();
479  ++iresp) {
480  // Eta-Pt handles -- restrict for iphi == 0
481  if (iphi == 0) {
482  auto& etaPtHandle = etaPtResidualPulls[ieta][ipt][iresp];
483  // Create the handle tag
484  TString handleTag(etaPtHandle.tag + etaTag + ptTag);
485  // Accept range and range cut
486  etaPtHandle.accept = etaPtBin;
487  etaPtHandle.rangeCutStr = ptCut + land + etaCut;
488  handleRange(etaPtHandle, handleTag, peakEntries);
489  bookHistograms(etaPtHandle, pullRange, nHistBins, ++histBarcode);
490 
491  // Set name and style/
492  TString residualN = TString("res_") + handleTag;
493  TString pullN = TString("pull_") + handleTag;
494  // Range histogram does not exist when read in from configuration
495  if (etaPtHandle.rangeHist != nullptr) {
496  TString rangeN = TString("range_") + handleTag;
497  setHistStyle(etaPtHandle.rangeHist);
498  etaPtHandle.rangeHist->SetName(rangeN);
499  }
500  setHistStyle(etaPtHandle.residualHist);
501  etaPtHandle.residualHist->SetName(residualN);
502  setHistStyle(etaPtHandle.pullHist);
503  etaPtHandle.pullHist->SetName(pullN);
504  }
505  // Eta-Phi handles --- restrice for ipt == 0
506  if (ipt == 0) {
507  auto& etaPhiHandle = etaPhiResidualPulls[ieta][iphi][iresp];
508  // Create the handle tag
509  TString handleTag(etaPhiHandle.tag + etaTag + phiTag);
510  etaPhiHandle.accept = etaPhiBin;
511  handleRange(etaPhiHandle, handleTag, peakEntries);
512  bookHistograms(etaPhiHandle, pullRange, nHistBins, ++histBarcode);
513 
514  // Set name and style
515  TString residualN = TString("res_") + handleTag;
516  TString pullN = TString("pull_") + handleTag;
517 
518  // Range histogram does not exist when read in from configuration
519  if (etaPhiHandle.rangeHist != nullptr) {
520  TString rangeN = TString("range_") + handleTag;
521  setHistStyle(etaPhiHandle.rangeHist);
522  etaPhiHandle.rangeHist->SetName(rangeN);
523  }
524  setHistStyle(etaPhiHandle.residualHist);
525  etaPhiHandle.residualHist->SetName(residualN);
526  setHistStyle(etaPhiHandle.pullHist);
527  etaPhiHandle.pullHist->SetName(pullN);
528  }
529 
530 #ifdef BOOST_AVAILABLE
531  ++handle_preparation_progress;
532 #endif
533  }
534 
535  // Auxiliary plots
536  for (unsigned int iaux = 0; iaux < baseAuxilaries.size(); ++iaux) {
537  // Eta-Pt handles - restrict to iphi == 0
538  if (iphi == 0) {
539  auto& etaPtAux = etaPtAuxiliaries[ieta][ipt][iaux];
540  etaPtAux.accept = etaPtBin;
541  TString handleTag(etaPtAux.tag + etaTag + ptTag);
542  etaPtAux.hist =
543  new TH1F(handleTag.Data(), etaPtAux.tag.c_str(), etaPtAux.bins,
544  etaPtAux.range[0], etaPtAux.range[1]);
545  etaPtAux.hist->GetXaxis()->SetTitle(etaPtAux.label.c_str());
546  etaPtAux.hist->GetYaxis()->SetTitle("Entries");
547  setHistStyle(etaPtAux.hist);
548  }
549 
550  // Eta-Phi handles - restrict to ipt == 0
551  if (ipt == 0) {
552  auto& etaPhiAux = etaPhiAuxiliaries[ieta][iphi][iaux];
553  etaPhiAux.accept = etaPhiBin;
554  TString handleTag(etaPhiAux.tag + etaTag + phiTag);
555  etaPhiAux.hist = new TH1F(handleTag.Data(), etaPhiAux.tag.c_str(),
556  etaPhiAux.bins, etaPhiAux.range[0],
557  etaPhiAux.range[1]);
558  etaPhiAux.hist->GetXaxis()->SetTitle(etaPhiAux.label.c_str());
559  etaPhiAux.hist->GetYaxis()->SetTitle("Entries");
560  setHistStyle(etaPhiAux.hist);
561  }
562  }
563  }
564  }
565  }
566 
567 #ifdef NLOHMANN_AVAILABLE
568  if (not outConfig.empty()) {
569  std::ofstream config_out;
570  config_out.open(outConfig.c_str());
571  config_out << handle_configs.dump(4);
572  }
573 #endif
574 
575 #ifdef BOOST_AVAILABLE
576  std::cout << "*** Event Loop: " << std::endl;
577  progress_display event_loop_progress(entries);
578 #endif
579 
580  for (unsigned long ie = 0; ie < entries; ++ie) {
581 #ifdef BOOST_AVAILABLE
582  ++event_loop_progress;
583 #endif
584 
585  // Make sure you have the entry
586  tracks.tree->GetEntry(ie);
587  size_t nTracks = tracks.hasFittedParams->size();
588  for (size_t it = 0; it < nTracks; ++it) {
589  if (tracks.hasFittedParams->at(it)) {
590  // Residual handlesL
591  // Full range handles
592  for (auto& fHandle : fullResidualPulls) {
593  fHandle.fill(it);
594  }
595 
596  // Eta-Pt handles
597  for (auto& etaBatch : etaPtResidualPulls) {
598  for (auto& ptBatch : etaBatch) {
599  for (auto& bHandle : ptBatch) {
600  bHandle.fill(it);
601  }
602  }
603  }
604 
605  // Eta-Phi handles
606  for (auto& etaBatch : etaPhiResidualPulls) {
607  for (auto& phiBatch : etaBatch) {
608  for (auto& bHandle : phiBatch) {
609  bHandle.fill(it);
610  }
611  }
612  }
613  }
614 
615  // Auxiliary handles:
616  // Full range handles
617  for (auto& fAuxiliary : fullAuxiliaries) {
618  fAuxiliary.fill(it);
619  }
620 
621  // Eta-Pt handles
622  for (auto& etaBatch : etaPtAuxiliaries) {
623  for (auto& ptBatch : etaBatch) {
624  for (auto& bHandle : ptBatch) {
625  bHandle.fill(it);
626  }
627  }
628  }
629 
630  // Eta-Phi handles
631  for (auto& etaBatch : etaPhiAuxiliaries) {
632  for (auto& phiBatch : etaBatch) {
633  for (auto& bHandle : phiBatch) {
634  bHandle.fill(it);
635  }
636  }
637  }
638  }
639  }
640 
641  // The output file section
642  auto output = TFile::Open(outFile.c_str(), "recreate");
643  output->cd();
644 
645  // Full range handles : residual and pulls
646  for (auto& fHandle : fullResidualPulls) {
647  if (fHandle.rangeHist != nullptr) {
648  fHandle.rangeHist->Write();
649  }
650  fHandle.residualHist->Write();
651  fHandle.pullHist->Write();
652  }
653 
654  // Full range handles : auxiliaries
655  for (auto& fAuxiliary : fullAuxiliaries) {
656  fAuxiliary.hist->SetName(fAuxiliary.hist->GetName() + TString("_all"));
657  fAuxiliary.hist->Write();
658  }
659 
660  struct SummaryHistograms {
661  TH2F* fillStats = nullptr;
662 
663  std::vector<TH2F*> residualRMS;
664  std::vector<TH2F*> residualMean;
665  std::vector<TH2F*> pullSigma;
666  std::vector<TH2F*> pullMean;
667 
668  std::vector<TH2F*> auxiliaries;
669  };
670 
681  auto analyseBins = [&](ResidualPullsMatrix& residualPullsMatrix,
682  AuxiliariesMatrix& auxiliaryMatrix,
683  const TString& matrixTag, const TVectorF& outerBorders,
684  const TVectorF& innerBorders,
685  const TString& fXTitle = "#eta",
686  const TString& sXTitle = "#phi") -> void {
687  // The summary histogram set
688  SummaryHistograms summary;
689 
690  // 2D handles ---------------------------
691  unsigned int nOuterBins = outerBorders.GetNrows() - 1;
692  auto outerValues = outerBorders.GetMatrixArray();
693  unsigned int nInnerBins = innerBorders.GetNrows() - 1;
694  auto innerValues = innerBorders.GetMatrixArray();
695 
696  TString statN = TString("entries") + matrixTag;
697  summary.fillStats =
698  new TH2F(statN, "", nOuterBins, outerValues, nInnerBins, innerValues);
699 
700 #ifdef BOOST_AVAILABLE
701  progress_display analysis_progress(nOuterBins * nInnerBins);
702 #endif
703 
704  // Prepare by looping over the base bHandles - residuals
705  for (auto& bHandle : baseResidualPulls) {
706  // Create a unique handle tag
707  TString handleTag = TString(bHandle.tag) + matrixTag;
708  // ... and names
709  TString residualRMSN = TString("res_rms_") + handleTag;
710  TString residualMeanN = TString("res_mean_") + handleTag;
711  TString pullSigmaN = TString("pull_sigma_") + handleTag;
712  TString pullMeanN = TString("pull_mean_") + handleTag;
713 
714  TH2F* residualRMS = new TH2F(residualRMSN, "", nOuterBins, outerValues,
715  nInnerBins, innerValues);
716  TH2F* residualMean = new TH2F(residualMeanN, "", nOuterBins, outerValues,
717  nInnerBins, innerValues);
718  TH2F* pullSigma = new TH2F(pullSigmaN, "", nOuterBins, outerValues,
719  nInnerBins, innerValues);
720  TH2F* pullMean = new TH2F(pullMeanN, "", nOuterBins, outerValues,
721  nInnerBins, innerValues);
722  // Booked & ready
723  summary.residualRMS.push_back(residualRMS);
724  summary.residualMean.push_back(residualMean);
725  summary.pullSigma.push_back(pullSigma);
726  summary.pullMean.push_back(pullMean);
727  }
728 
729  // Prepare by looping over the base handles - auxiliaries
730  for (auto& aHandle : baseAuxilaries) {
731  // Create a unique handle tag
732  TString auxiliaryTag = TString(aHandle.tag) + matrixTag;
733  TH2F* auxHist = new TH2F(auxiliaryTag, auxiliaryTag, nOuterBins,
734  outerValues, nInnerBins, innerValues);
735  summary.auxiliaries.push_back(auxHist);
736  }
737 
738  unsigned int io = 0;
739  for (auto& outerBatch : residualPullsMatrix) {
740  unsigned int ii = 0;
741  for (auto& innerBatch : outerBatch) {
742  // residual/pull loop
743  unsigned int iresp = 0;
744  for (auto& bHandle : innerBatch) {
745  // Range estimates / could be empty if read from configuration
746  if (bHandle.rangeHist != nullptr) {
747  bHandle.rangeHist->Write();
748  }
749  // Fill the stats
750  if (iresp == 0) {
751  summary.fillStats->SetBinContent(io + 1, ii + 1, bHandle.accepted);
752  }
753 
754  // Residuals
755  // Get RMS/ RMSError
756  float rrms = bHandle.residualHist->GetRMS();
757  float rrerr = bHandle.residualHist->GetRMSError();
758  float rmean = bHandle.residualHist->GetMean();
759  float rmerr = bHandle.residualHist->GetMeanError();
760  summary.residualRMS[iresp]->SetBinContent(io + 1, ii + 1, rrms);
761  summary.residualRMS[iresp]->SetBinError(io + 1, ii + 1, rrerr);
762  summary.residualMean[iresp]->SetBinContent(io + 1, ii + 1, rmean);
763  summary.residualMean[iresp]->SetBinError(io + 1, ii + 1, rmerr);
764  bHandle.residualHist->Write();
765  // Pulls
766  bHandle.pullHist->Fit("gaus", "q");
767  TF1* gauss = bHandle.pullHist->GetFunction("gaus");
768  if (gauss != nullptr) {
769  float pmu = gauss->GetParameter(1);
770  float pmerr = gauss->GetParError(1);
771  float psigma = gauss->GetParameter(2);
772  float pserr = gauss->GetParError(2);
773  summary.pullSigma[iresp]->SetBinContent(io + 1, ii + 1, psigma);
774  summary.pullSigma[iresp]->SetBinError(io + 1, ii + 1, pserr);
775  summary.pullMean[iresp]->SetBinContent(io + 1, ii + 1, pmu);
776  summary.pullMean[iresp]->SetBinError(io + 1, ii + 1, pmerr);
777  }
778  bHandle.pullHist->Write();
779  ++iresp;
780  }
781 
782  // auxilaiary loop
783  auto auxiliaryBatch = auxiliaryMatrix[io][ii];
784  unsigned int iaux = 0;
785  for (auto& aHandle : auxiliaryBatch) {
786  float value = aHandle.hist->GetMean();
787  float error = aHandle.hist->GetMeanError();
788  summary.auxiliaries[iaux]->SetBinContent(io + 1, ii + 1, value);
789  summary.auxiliaries[iaux]->SetBinError(io + 1, ii + 1, error);
790  aHandle.hist->Write();
791 
792  ++iaux;
793  }
794 #ifdef BOOST_AVAILABLE
795  ++analysis_progress;
796 #endif
797  ++ii;
798  }
799  ++io;
800  }
801 
809  auto writeProjections = [](const TH2F& h2, const TString& fXTitleP = "#eta",
810  const TString& fYTitleP = "sigma",
811  const TString& sXTitleP = "#phi",
812  const TString& sYTitleP = "sigma") -> void {
813  const TString& fTag = "_pX";
814  const TString& sTag = "_pY";
815 
816  unsigned int nBinsX = h2.GetXaxis()->GetNbins();
817  unsigned int nBinsY = h2.GetYaxis()->GetNbins();
818  if (fTag != "") {
819  TH1D* pX =
820  dynamic_cast<TH1D*>(h2.ProjectionX((h2.GetName() + fTag).Data()));
821  setHistStyle(pX);
822  if (pX != nullptr) {
823  pX->GetXaxis()->SetTitle(fXTitleP.Data());
824  pX->GetYaxis()->SetTitle(fYTitleP.Data());
825  pX->Write();
826  }
827  // Bin-wise projections
828  for (unsigned int iy = 1; iy <= nBinsY; ++iy) {
829  pX = dynamic_cast<TH1D*>(h2.ProjectionX(
830  (h2.GetName() + fTag + sTag + (iy - 1)).Data(), iy, iy));
831  setHistStyle(pX);
832  if (pX != nullptr) {
833  pX->GetXaxis()->SetTitle(fXTitleP.Data());
834  pX->GetYaxis()->SetTitle(fYTitleP.Data());
835  pX->Write();
836  }
837  }
838  }
839  if (sTag != "") {
840  TH1D* pY =
841  dynamic_cast<TH1D*>(h2.ProjectionY((h2.GetName() + sTag).Data()));
842  setHistStyle(pY);
843  if (pY != nullptr) {
844  pY->GetXaxis()->SetTitle(sXTitleP.Data());
845  pY->GetYaxis()->SetTitle(sYTitleP.Data());
846  pY->Write();
847  }
848  // Bin-wise projections
849  for (unsigned int ix = 1; ix <= nBinsX; ++ix) {
850  pY = dynamic_cast<TH1D*>(h2.ProjectionY(
851  (h2.GetName() + sTag + fTag + (ix - 1)).Data(), ix, ix));
852  setHistStyle(pY);
853  if (pY != nullptr) {
854  pY->GetXaxis()->SetTitle(sXTitleP.Data());
855  pY->GetYaxis()->SetTitle(sYTitleP.Data());
856  pY->Write();
857  }
858  }
859  }
860  };
861 
862  setHistStyle(summary.fillStats);
863  summary.fillStats->Write();
864 
865  // Write mapped residual/pull histograms and their projections
866  for (unsigned int iresp = 0; iresp < baseResidualPulls.size(); ++iresp) {
867  // Get the handle for writing out
868  auto bHandle = baseResidualPulls[iresp];
869 
870  TString rrms = TString("RMS[") + bHandle.residualStr + TString("] ") +
871  bHandle.residualUnit;
872  TString rmu = TString("#mu[") + bHandle.residualStr + TString("] ") +
873  bHandle.residualUnit;
874 
875  TString psigma = TString("#sigma[(") + bHandle.residualStr +
876  TString(")/") + bHandle.errorStr + TString("]");
877  TString pmu = TString("#mu[(") + bHandle.residualStr + TString(")/") +
878  bHandle.errorStr + TString("]");
879 
880  // 2D map histograms
881  setHistStyle(summary.residualRMS[iresp]);
882  summary.residualRMS[iresp]->GetXaxis()->SetTitle(fXTitle);
883  summary.residualRMS[iresp]->GetYaxis()->SetTitle(sXTitle);
884  summary.residualRMS[iresp]->GetZaxis()->SetTitle(rrms);
885  summary.residualRMS[iresp]->Write();
886 
887  setHistStyle(summary.residualMean[iresp]);
888  summary.residualMean[iresp]->GetXaxis()->SetTitle(fXTitle);
889  summary.residualMean[iresp]->GetYaxis()->SetTitle(sXTitle);
890  summary.residualMean[iresp]->GetZaxis()->SetTitle(rmu);
891  summary.residualMean[iresp]->Write();
892 
893  setHistStyle(summary.pullSigma[iresp]);
894  adaptColorPalette(summary.pullSigma[iresp], 0., 4., 1., 0.1, 104);
895  summary.pullSigma[iresp]->GetXaxis()->SetTitle(fXTitle);
896  summary.pullSigma[iresp]->GetYaxis()->SetTitle(sXTitle);
897  summary.pullSigma[iresp]->GetZaxis()->SetRangeUser(0., 4.);
898  summary.pullSigma[iresp]->GetZaxis()->SetTitle(psigma);
899  summary.pullSigma[iresp]->Write();
900 
901  setHistStyle(summary.pullMean[iresp]);
902  adaptColorPalette(summary.pullMean[iresp], -1., 1., 0., 0.1, 104);
903  summary.pullMean[iresp]->GetXaxis()->SetTitle(fXTitle);
904  summary.pullMean[iresp]->GetYaxis()->SetTitle(sXTitle);
905  summary.pullMean[iresp]->GetZaxis()->SetRangeUser(-1., 1.);
906  summary.pullMean[iresp]->GetZaxis()->SetTitle(pmu);
907  summary.pullMean[iresp]->Write();
908 
909  // Write the projection histograms
910  writeProjections(*summary.residualRMS[iresp], fXTitle, rrms, sXTitle,
911  rrms);
912  writeProjections(*summary.residualMean[iresp], fXTitle, rmu, sXTitle,
913  rmu);
914  writeProjections(*summary.pullSigma[iresp], fXTitle, psigma, sXTitle,
915  psigma);
916  writeProjections(*summary.pullMean[iresp], fXTitle, pmu, sXTitle, pmu);
917  }
918 
919  // Write mapped auxiliary histograms and their projections
920  for (unsigned int iaux = 0; iaux < baseAuxilaries.size(); ++iaux) {
921  setHistStyle(summary.auxiliaries[iaux]);
922  summary.auxiliaries[iaux]->GetXaxis()->SetTitle(fXTitle);
923  summary.auxiliaries[iaux]->GetYaxis()->SetTitle(sXTitle);
924  summary.auxiliaries[iaux]->GetZaxis()->SetTitle(
925  baseAuxilaries[iaux].label.c_str());
926  summary.auxiliaries[iaux]->Write();
927 
928  writeProjections(*summary.auxiliaries[iaux], fXTitle,
929  baseAuxilaries[iaux].label, sXTitle,
930  baseAuxilaries[iaux].label);
931  }
932 
933  return;
934  };
935 
936 // The handle matrices
937 #ifdef BOOST_AVAILABLE
938  std::cout << "*** Bin/Projection Analysis: " << std::endl;
939 #endif
940  analyseBins(etaPtResidualPulls, etaPtAuxiliaries, TString("_eta_pt"), etaVals,
941  ptVals, "#eta", "p_{T} [GeV]");
942  analyseBins(etaPhiResidualPulls, etaPhiAuxiliaries, TString("_eta_phi"),
943  etaVals, phiVals, "#eta", "#phi");
944 
945  output->Close();
946 
947  return 1;
948 }