Belle II Software  release-08-01-10
TOPValidationAlgorithm.cc
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 
9 #include <top/calibration/TOPValidationAlgorithm.h>
10 #include <top/utilities/Chi2MinimumFinder1D.h>
11 
12 #include <TFile.h>
13 #include <TH1F.h>
14 #include <TH1D.h>
15 #include <TH2F.h>
16 
17 #include <string>
18 #include <vector>
19 #include <map>
20 #include <set>
21 #include <cmath>
22 
23 using namespace std;
24 
25 namespace Belle2 {
30  namespace TOP {
31 
32  TOPValidationAlgorithm::TOPValidationAlgorithm():
33  CalibrationAlgorithm("TOPValidationCollector")
34  {
35  setDescription("Calibration algorithm for automatic validation of the calibration");
36  }
37 
38 
40  {
41  // open output root file
42 
43  auto* file = TFile::Open("TOPCalValidation.root", "recreate");
44 
45  // write-out the histograms that are ready
46 
47  B2INFO("TOPValidationAlgorithm: get basic histograms");
48 
49  double scaleFactor = 1; // scaling factor for module and common T0 uncertainties
50  auto h1 = getObjectPtr<TH1F>("moduleT0_pulls");
51  if (h1) {
52  if (h1->GetEntries() > 1) scaleFactor = h1->GetRMS();
53  h1->Write();
54  }
55 
56  std::vector<string> slotNames;
57  for (unsigned slot = 1; slot <= c_numModules; slot++) {
58  string slotName = "slot";
59  if (slot < 10) slotName += "0";
60  slotName += to_string(slot);
61  slotNames.push_back(slotName);
62  }
63 
64  for (unsigned slot = 0; slot < c_numModules; slot++) {
65  string name = "hits_" + slotNames[slot];
66  auto h = getObjectPtr<TH2F>(name);
67  if (h) h->Write();
68  }
69 
70  // channelT0 residuals: determine Chi2 minima and save the results into histograms
71 
72  B2INFO("TOPValidationAlgorithm: determine channel T0 residuals");
73 
74  for (unsigned slot = 0; slot < c_numModules; slot++) {
75  string name = "chi2_" + slotNames[slot];
76  auto h = getObjectPtr<TH2F>(name);
77  if (not h) continue;
78 
79  int nx = h->GetNbinsX();
80  string nameOut = "channelT0_" + slotNames[slot];
81  string title = "ChannelT0, " + slotNames[slot] + "; channel; residuals [ns]";
82  auto* h0 = new TH1F(nameOut.c_str(), title.c_str(), nx, 0, nx);
83 
84  for (int ibin = 1; ibin <= nx; ibin++) {
85  Chi2MinimumFinder1D finder(std::shared_ptr<TH1D>(h->ProjectionY("py", ibin, ibin)));
86  const auto& minimum = finder.getMinimum();
87  if (not minimum.valid) continue;
88  h0->SetBinContent(ibin, minimum.position);
89  h0->SetBinError(ibin, minimum.error);
90  }
91  }
92 
93  // moduleT0, commonT0 and others: get input tree and set branch addresses
94 
95  B2INFO("TOPValidationAlgorithm: get input tree");
96 
97  auto inputTree = getObjectPtr<TTree>("tree");
98  if (not inputTree) {
99  B2ERROR("TOPValidationAlgorithm: input tree not found");
100  file->Write();
101  file->Close();
102  return c_Failure;
103  }
104  inputTree->SetBranchAddress("expNo", &m_treeEntry.expNo);
105  inputTree->SetBranchAddress("runNo", &m_treeEntry.runNo);
106  inputTree->SetBranchAddress("numTracks", &m_treeEntry.numTracks);
107  inputTree->SetBranchAddress("commonT0", &m_treeEntry.commonT0);
108  inputTree->SetBranchAddress("commonT0Err", &m_treeEntry.commonT0Err);
109  inputTree->SetBranchAddress("moduleT0", &m_treeEntry.moduleT0);
110  inputTree->SetBranchAddress("moduleT0Err", &m_treeEntry.moduleT0Err);
111  inputTree->SetBranchAddress("numTBCalibrated", &m_treeEntry.numTBCalibrated);
112  inputTree->SetBranchAddress("numT0Calibrated", &m_treeEntry.numT0Calibrated);
113  inputTree->SetBranchAddress("numActive", &m_treeEntry.numActive);
114  inputTree->SetBranchAddress("numActiveCalibrated", &m_treeEntry.numActiveCalibrated);
115  inputTree->SetBranchAddress("thrEffi", &m_treeEntry.thrEffi);
116  inputTree->SetBranchAddress("asicShifts", &m_treeEntry.asicShifts);
117  inputTree->SetBranchAddress("svdOffset", &m_treeEntry.svdOffset);
118  inputTree->SetBranchAddress("svdSigma", &m_treeEntry.svdSigma);
119  inputTree->SetBranchAddress("cdcOffset", &m_treeEntry.cdcOffset);
120  inputTree->SetBranchAddress("cdcSigma", &m_treeEntry.cdcSigma);
121  inputTree->SetBranchAddress("fillPatternOffset", &m_treeEntry.fillPatternOffset);
122  inputTree->SetBranchAddress("fillPatternFraction", &m_treeEntry.fillPatternFraction);
123 
124  // sort input tree entries according to exp/run
125 
126  B2INFO("TOPValidationAlgorithm: sort input tree entries");
127 
128  std::multimap<int, ValidationTreeStruct> sortedEntries;
129  std::set<int> sortedRuns;
130  for (int iev = 0; iev < inputTree->GetEntries(); iev++) {
131  inputTree->GetEntry(iev);
132  int expRun = (m_treeEntry.expNo << 16) + m_treeEntry.runNo;
133  sortedEntries.emplace(expRun, m_treeEntry);
134  sortedRuns.insert(expRun);
135  }
136 
137  // merge input entries of the same exp/run and rescale the errors
138 
139  B2INFO("TOPValidationAlgorithm: merge input tree entries");
140 
141  std::vector<ValidationTreeStruct> mergedEntries;
142  for (auto expRun : sortedRuns) {
143  ValidationTreeStruct mergedEntry;
144  const auto range = sortedEntries.equal_range(expRun);
145  for (auto it = range.first; it != range.second; ++it) {
146  mergedEntry.merge(it->second);
147  }
148  mergedEntry.rescaleErrors(scaleFactor);
149  mergedEntries.push_back(mergedEntry);
150  }
151 
152  // make histograms vs. run index
153 
154  int nx = mergedEntries.size();
155  if (nx == 0) {
156  B2ERROR("TOPValidationAlgorithm: input tree is empty");
157  file->Write();
158  file->Close();
159  return c_Failure;
160  }
161 
162  file->cd();
163 
164  std::set<int> sortedExps;
165  for (auto& entry : mergedEntries) sortedExps.insert(entry.expNo);
166  std::string titleIndex = "Experiment";
167  if (sortedExps.size() > 1) titleIndex += "s";
168  for (auto expNo : sortedExps) titleIndex += " " + to_string(expNo) + ",";
169  titleIndex.pop_back();
170  auto* h_index = new TH1F("runIndex", (titleIndex + "; run index; run number").c_str(), nx, 0, nx);
171  for (int i = 0; i < nx; i++) {
172  h_index->SetBinContent(i + 1, mergedEntries[i].runNo);
173  }
174 
175  auto* h_numTracks = new TH1F("numTracks", "Number of tracks; run index; number of tracks", nx, 0, nx);
176  for (int i = 0; i < nx; i++) {
177  h_numTracks->SetBinContent(i + 1, mergedEntries[i].numTracks);
178  }
179 
180  auto* h_numMerged = new TH1F("numMerged", "Number of merged entries; run index; number of merged", nx, 0, nx);
181  for (int i = 0; i < nx; i++) {
182  h_numMerged->SetBinContent(i + 1, mergedEntries[i].numMerged);
183  }
184 
185  auto* h_commonT0 = new TH1F("commonT0", "commonT0; run index; residuals [ns]", nx, 0, nx);
186  for (int i = 0; i < nx; i++) {
187  h_commonT0->SetBinContent(i + 1, mergedEntries[i].commonT0);
188  h_commonT0->SetBinError(i + 1, mergedEntries[i].commonT0Err);
189  }
190 
191  for (unsigned slot = 0; slot < c_numModules; slot++) {
192  string name = "moduleT0_" + slotNames[slot];
193  string title = "moduleT0, " + slotNames[slot] + "; run index; residuals [ns]";
194  auto* h = new TH1F(name.c_str(), title.c_str(), nx, 0, nx);
195  for (int i = 0; i < nx; i++) {
196  h->SetBinContent(i + 1, mergedEntries[i].moduleT0[slot]);
197  h->SetBinError(i + 1, mergedEntries[i].moduleT0Err[slot]);
198  }
199  }
200 
201  for (unsigned slot = 0; slot < c_numModules; slot++) {
202  string name = "numTBCalibrated_" + slotNames[slot];
203  string title = "Time base calibrated, " + slotNames[slot] + "; run index; fraction";
204  auto* h = new TH1F(name.c_str(), title.c_str(), nx, 0, nx);
205  for (int i = 0; i < nx; i++) {
206  h->SetBinContent(i + 1, mergedEntries[i].numTBCalibrated[slot] / 512.0);
207  }
208  }
209 
210  for (unsigned slot = 0; slot < c_numModules; slot++) {
211  string name = "numT0Calibrated_" + slotNames[slot];
212  string title = "channel T0 calibrated, " + slotNames[slot] + "; run index; fraction";
213  auto* h = new TH1F(name.c_str(), title.c_str(), nx, 0, nx);
214  for (int i = 0; i < nx; i++) {
215  h->SetBinContent(i + 1, mergedEntries[i].numT0Calibrated[slot] / 512.0);
216  }
217  }
218 
219  for (unsigned slot = 0; slot < c_numModules; slot++) {
220  string name = "numActive_" + slotNames[slot];
221  string title = "Active, " + slotNames[slot] + "; run index; fraction";
222  auto* h = new TH1F(name.c_str(), title.c_str(), nx, 0, nx);
223  for (int i = 0; i < nx; i++) {
224  h->SetBinContent(i + 1, mergedEntries[i].numActive[slot] / 512.0);
225  }
226  }
227 
228  for (unsigned slot = 0; slot < c_numModules; slot++) {
229  string name = "numActiveCalibrated_" + slotNames[slot];
230  string title = "Active and calibrated, " + slotNames[slot] + "; run index; fraction";
231  auto* h = new TH1F(name.c_str(), title.c_str(), nx, 0, nx);
232  for (int i = 0; i < nx; i++) {
233  h->SetBinContent(i + 1, mergedEntries[i].numActiveCalibrated[slot] / 512.0);
234  }
235  }
236 
237  for (unsigned slot = 0; slot < c_numModules; slot++) {
238  string name = "thrEffi_" + slotNames[slot];
239  string title = "Threshold efficiency, " + slotNames[slot] + "; run index; efficiency";
240  auto* h = new TH1F(name.c_str(), title.c_str(), nx, 0, nx);
241  for (int i = 0; i < nx; i++) {
242  h->SetBinContent(i + 1, mergedEntries[i].thrEffi[slot]);
243  }
244  }
245 
246  for (unsigned carrier = 0; carrier < 4; carrier++) {
247  string name = "asicShifts_" + to_string(carrier);
248  string title = "BS13d, carrier " + to_string(carrier) + "; run index; shift [ns]";
249  auto* h = new TH1F(name.c_str(), title.c_str(), nx, 0, nx);
250  for (int i = 0; i < nx; i++) {
251  auto y = mergedEntries[i].asicShifts[carrier];
252  if (isnan(y)) {
253  h->SetBinError(i + 1, 0);
254  } else {
255  h->SetBinContent(i + 1, y);
256  h->SetBinError(i + 1, 0.001);
257  }
258  }
259  }
260 
261  auto* h_svdOffset = new TH1F("svdOffset", "Event T0 offset w.r.t TOP; run index; offset [ns]", nx, 0, nx);
262  auto* h_svdSigma = new TH1F("svdSigma", "Event T0 resolution; run index; sigma [ns]", nx, 0, nx);
263  for (int i = 0; i < nx; i++) {
264  h_svdOffset->SetBinContent(i + 1, mergedEntries[i].svdOffset);
265  h_svdOffset->SetBinError(i + 1, mergedEntries[i].svdSigma / 1000);
266  h_svdSigma->SetBinContent(i + 1, mergedEntries[i].svdSigma);
267  h_svdSigma->SetBinError(i + 1, 0.001);
268  }
269 
270  auto* h_cdcOffset = new TH1F("cdcOffset", "Event T0 offset w.r.t TOP; run index; offset [ns]", nx, 0, nx);
271  auto* h_cdcSigma = new TH1F("cdcSigma", "Event T0 resolution; run index; sigma [ns]", nx, 0, nx);
272  for (int i = 0; i < nx; i++) {
273  h_cdcOffset->SetBinContent(i + 1, mergedEntries[i].cdcOffset);
274  h_cdcOffset->SetBinError(i + 1, mergedEntries[i].cdcSigma / 1000);
275  h_cdcSigma->SetBinContent(i + 1, mergedEntries[i].cdcSigma);
276  h_cdcSigma->SetBinError(i + 1, 0.001);
277  }
278 
279  auto* h_fillPatternOffset = new TH1F("fillPatternOffset", "Fill pattern offset; run index; offset [RF clk]", nx, 0, nx);
280  auto* h_fillPatternFract = new TH1F("fillPatternFract",
281  "Fraction of buckets matched with filled ones; run index; fraction", nx, 0, nx);
282  for (int i = 0; i < nx; i++) {
283  if (isnan(mergedEntries[i].fillPatternOffset)) {
284  h_fillPatternOffset->SetBinError(i + 1, 0);
285  } else {
286  h_fillPatternOffset->SetBinContent(i + 1, mergedEntries[i].fillPatternOffset);
287  h_fillPatternOffset->SetBinError(i + 1, 0.01);
288  }
289  h_fillPatternFract->SetBinContent(i + 1, mergedEntries[i].fillPatternFraction);
290  }
291 
292  // write-out the results
293 
294  B2INFO("TOPValidationAlgorithm: write the results");
295 
296  file->Write();
297  file->Close();
298 
299  return c_OK;
300  }
301 
302  } // end namespace TOP
304 } // end namespace Belle2
Base class for calibration algorithms.
void setDescription(const std::string &description)
Set algorithm description (in constructor)
EResult
The result of calibration.
@ c_OK
Finished successfuly =0 in Python.
@ c_Failure
Failed =3 in Python.
Minimum finder using tabulated chi^2 values in one dimension.
virtual EResult calibrate() final
algorithm implementation
ValidationTreeStruct m_treeEntry
input tree entry
Abstract base class for different kinds of events.
Calibration validation tree structure.
int numTBCalibrated[c_numModules]
number of timebase calibrated channels, index = slot - 1
void rescaleErrors(double scaleFactor)
Rescale errors.
int numTracks
number of selected tracks
int numActive[c_numModules]
number of active channels, index = slot - 1
float svdOffset
SVD event T0 offset.
float svdSigma
SVD event T0 resolution.
float fillPatternFraction
fraction of reconstructed buckets matched with filled ones
void merge(const ValidationTreeStruct &other)
Merge two structures.
float fillPatternOffset
fill pattern offset
float commonT0Err
common T0 uncertainty (not scaled)
float thrEffi[c_numModules]
threshold efficiency: average over active calibrated channels, index = slot - 1
float moduleT0[c_numModules]
module T0 residuals, index = slot - 1
int numActiveCalibrated[c_numModules]
number of active calibrated channels, index = slot - 1
float cdcSigma
CDC event T0 resolution.
float cdcOffset
CDC event T0 offset.
float moduleT0Err[c_numModules]
module T0 uncertainties (not scaled), index = slot - 1
float asicShifts[4]
carrier shifts of BS13d, index = carrier number
int numT0Calibrated[c_numModules]
number of channel T0 calibrated channels, index = slot - 1