Belle II Software development
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
23using namespace std;
24
25namespace Belle2 {
30 namespace TOP {
31
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
void setDescription(const std::string &description)
Set algorithm description (in constructor)
EResult
The result of calibration.
@ c_OK
Finished successfully =0 in Python.
CalibrationAlgorithm(const std::string &collectorModuleName)
Constructor - sets the prefix for collected objects (won't be accesses until execute(....
Minimum finder using tabulated chi^2 values in one dimension.
virtual EResult calibrate() final
algorithm implementation
ValidationTreeStruct m_treeEntry
input tree entry
std::shared_ptr< T > getObjectPtr(const std::string &name, const std::vector< Calibration::ExpRun > &requestedRuns)
Get calibration data object by name and list of runs, the Merge function will be called to generate t...
Abstract base class for different kinds of events.
STL namespace.
Calibration validation tree structure.
void rescaleErrors(double scaleFactor)
Rescale errors.
void merge(const ValidationTreeStruct &other)
Merge two structures.