Belle II Software development
TOPLocalCalFitter.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/TOPLocalCalFitter.h>
10
11// C++ STL
12#include <algorithm>
13#include <cmath>
14#include <iomanip>
15#include <limits>
16
17// ROOT
18#include <TF1.h>
19#include <TFile.h>
20#include <TH2F.h>
21#include <TMath.h>
22#include <TROOT.h>
23#include <TTree.h>
24
25// Belle II
26#include <framework/logging/Logger.h>
27#include <top/dbobjects/TOPCalChannelT0.h>
28
29void Belle2::TOP::TOPLocalCalFitter::fitChannel(short iSlot, short iChannel, TH1* h_profile)
30{
31 fitChannel(iSlot, iChannel, h_profile, /*inBins=*/false, /*frac=*/0.0);
32}
33
34void Belle2::TOP::TOPLocalCalFitter::buildChannelMaps()
35{
36 if (!m_treeTTS) {
37 B2ERROR("TOPLocalCalFitter::buildChannelMaps called with null m_treeTTS.");
38 return;
39 }
40
41 // The TTS tree is indexed as (channel + 512 * slot), 0-based for both.
42 for (short slot = 0; slot < 16; ++slot) {
43 for (short ch = 0; ch < 512; ++ch) {
44 const Long64_t idx = static_cast<Long64_t>(ch) + 512LL * slot;
45 m_treeTTS->GetEntry(idx);
46 m_rowOf[slot][ch] = static_cast<short>(m_pixelRow);
47 m_colOf[slot][ch] = static_cast<short>(m_pixelCol);
48 }
49 }
50 m_hasChannelMaps = true;
51}
52
53// Destructor
54Belle2::TOP::TOPLocalCalFitter::~TOPLocalCalFitter()
55{
56 // Close & delete input files
57 if (m_inputTTS) {
58 m_inputTTS->Close();
59 delete m_inputTTS;
60 m_inputTTS = nullptr;
61 }
62 if (m_inputConstraints) {
63 m_inputConstraints->Close();
64 delete m_inputConstraints;
65 m_inputConstraints = nullptr;
66 }
67
68 // Delete trees if still around (ROOT does not auto-delete in-memory TTrees)
69 if (m_fitTree) {
70 delete m_fitTree;
71 m_fitTree = nullptr;
72 }
73 if (m_timewalkTree) {
74 delete m_timewalkTree;
75 m_timewalkTree = nullptr;
76 }
77 if (m_crosstalkTree) {
78 delete m_crosstalkTree;
79 m_crosstalkTree = nullptr;
80 }
81 if (m_fitTree_noXtalk) {
82 delete m_fitTree_noXtalk;
83 m_fitTree_noXtalk = nullptr;
84 }
85
86 // Close & delete output file last
87 if (m_histFile) {
88 m_histFile->Close();
89 delete m_histFile;
90 m_histFile = nullptr;
91 }
92}
93
94static double crystalball_function(double x, double alpha, double n, double sigma, double mean)
95{
96 // evaluate the crystal ball function
97 if (sigma < 0.) return 0.;
98 double z = (x - mean) / sigma;
99 if (alpha < 0) z = -z;
100 double abs_alpha = std::abs(alpha);
101 if (z > - abs_alpha)
102 return std::exp(- 0.5 * z * z);
103 else {
104 double nDivAlpha = n / abs_alpha;
105 double AA = std::exp(-0.5 * abs_alpha * abs_alpha);
106 double B = nDivAlpha - abs_alpha;
107 double arg = nDivAlpha / (B - z);
108 return AA * std::pow(arg, n);
109 }
110}
111
112static double crystalball_pdf(double x, double alpha, double n, double sigma, double mean)
113{
114 // evaluation of the PDF ( is defined only for n >1)
115 if (sigma < 0.) return 0.;
116 if (n <= 1) return std::numeric_limits<double>::quiet_NaN(); // pdf is not normalized for n <=1
117 double abs_alpha = std::abs(alpha);
118 double C = n / abs_alpha * 1. / (n - 1.) * std::exp(-alpha * alpha / 2.);
119 double D = std::sqrt(M_PI / 2.) * (1. + erf(abs_alpha / std::sqrt(2.)));
120 double N = 1. / (sigma * (C + D));
121 return N * crystalball_function(x, alpha, n, sigma, mean);
122}
123
124// Two gaussians to model the TTS of the MCP-PMTs.
125static double TTSPDF(double x, double time, double deltaT, double sigma1, double sigma2, double f1)
126{
127 return f1 * TMath::Gaus(x, time, sigma1, kTRUE) + (1 - f1) * TMath::Gaus(x, time + deltaT, sigma2, kTRUE) ;
128}
129
130// Full PDF to fit the laser, made of:
131// 2 TTSPDF
132// 1 crystal ball PDF for the extra peak at +1 ns we don't understand
133// 1 gaussian to help modelling the tail
134// cppcheck-suppress constParameterCallback
135static double laserPDF(double* x, double* p)
136{
137 // Define parameters
138 double time = p[0];
139 double sigma = p[1];
140 double fraction = p[2];
141 double deltaTLaser = p[3];
142 double sigmaRatio = p[4];
143 double deltaTTS = p[5];
144 double f1 = p[6];
145 double norm = p[7];
146 double deltaTExtra = p[8];
147 double sigmaExtra = p[9];
148 double normExtra = p[10];
149 double deltaTBkg = p[11];
150 double sigmaBkg = p[12];
151 double bkg = p[13];
152 double alpha = p[14];
153 double n = p[15];
154
155 // Define function
156 double mainPeak = fraction * TTSPDF(x[0], time, deltaTTS, sigma, TMath::Sqrt(sigma * sigma + sigmaRatio * sigmaRatio), f1);
157 double secondaryPeak = (1. - fraction) * TTSPDF(x[0], time + deltaTLaser, deltaTTS, sigma,
158 TMath::Sqrt(sigma * sigma + sigmaRatio * sigmaRatio), f1);
159 double extraPeak = crystalball_pdf(x[0], alpha, n, sigmaExtra, time + deltaTExtra);
160 double background = TMath::Gaus(x[0], time + deltaTBkg + deltaTTS, sigmaBkg, kTRUE);
161
162 return norm * (mainPeak + secondaryPeak) + normExtra * extraPeak + bkg * background;
163}
164
165Belle2::TOP::TOPLocalCalFitter::TOPLocalCalFitter(): CalibrationAlgorithm("TOPLaserCalibratorCollector")
166{
167 setDescription(
168 "Perform the fit of the laser and pulser runs"
169 );
170
171}
172
173void Belle2::TOP::TOPLocalCalFitter::loadMCInfoTrees()
174{
175 m_inputTTS = TFile::Open(m_TTSData.c_str());
176 m_inputConstraints = TFile::Open(m_fitConstraints.c_str());
177
178 B2INFO("Getting the TTS parameters from " << m_TTSData);
179 m_inputTTS->cd();
180 m_inputTTS->GetObject("tree", m_treeTTS);
181 m_treeTTS->SetBranchAddress("mean2", &m_mean2);
182 m_treeTTS->SetBranchAddress("sigma1", &m_sigma1);
183 m_treeTTS->SetBranchAddress("sigma2", &m_sigma2);
184 m_treeTTS->SetBranchAddress("fraction1", &m_f1);
185 m_treeTTS->SetBranchAddress("fraction2", &m_f2);
186 m_treeTTS->SetBranchAddress("pixelRow", &m_pixelRow);
187 m_treeTTS->SetBranchAddress("pixelCol", &m_pixelCol);
188
189 buildChannelMaps(); // build the maps rowOf[slot][channel], colOf[slot][channel]
190
191 if (m_fitterMode == "MC")
192 std::cout << "Running in MC mode, not constraints will be set" << std::endl;
193 else {
194 B2INFO("Getting the laser fit parameters from " << m_fitConstraints);
195 m_inputConstraints->cd();
196 m_inputConstraints->GetObject("fitTree", m_treeConstraints);
197 m_treeConstraints->SetBranchAddress("peakTime", &m_peakTimeConstraints);
198 m_treeConstraints->SetBranchAddress("deltaT", &m_deltaTConstraints);
199 m_treeConstraints->SetBranchAddress("fraction", &m_fractionConstraints);
200 if (m_fitterMode == "monitoring") {
201 m_treeConstraints->SetBranchAddress("timeExtra", &m_timeExtraConstraints);
202 m_treeConstraints->SetBranchAddress("sigmaExtra", &m_sigmaExtraConstraints);
203 m_treeConstraints->SetBranchAddress("alphaExtra", &m_alphaExtraConstraints);
204 m_treeConstraints->SetBranchAddress("nExtra", &m_nExtraConstraints);
205 m_treeConstraints->SetBranchAddress("timeBackground", &m_timeBackgroundConstraints);
206 m_treeConstraints->SetBranchAddress("sigmaBackground", &m_sigmaBackgroundConstraints);
207 }
208 }
209 return;
210}
211
212void Belle2::TOP::TOPLocalCalFitter::setupOutputTreeAndFile()
213{
214 m_histFile = new TFile(m_output.c_str(), "recreate");
215 m_histFile->cd();
216 m_fitTree = new TTree("fitTree", "fitTree");
217 m_fitTree->Branch<short>("channel", &m_channel);
218 m_fitTree->Branch<short>("slot", &m_slot);
219 m_fitTree->Branch<short>("row", &m_row);
220 m_fitTree->Branch<short>("col", &m_col);
221 m_fitTree->Branch<short>("asic", &m_asic);
222 m_fitTree->Branch<short>("asicChannel", &m_asicChannel);
223 m_fitTree->Branch<short>("boardstack", &m_boardstack);
224 m_fitTree->Branch<float>("peakTime", &m_peakTime);
225 m_fitTree->Branch<float>("peakTimeErr", &m_peakTimeErr);
226 m_fitTree->Branch<float>("deltaT", &m_deltaT);
227 m_fitTree->Branch<float>("deltaTErr", &m_deltaTErr);
228 m_fitTree->Branch<float>("sigma", &m_sigma);
229 m_fitTree->Branch<float>("sigmaErr", &m_sigmaErr);
230 m_fitTree->Branch<float>("fraction", &m_fraction);
231 m_fitTree->Branch<float>("fractionErr", &m_fractionErr);
232 m_fitTree->Branch<float>("yieldLaser", &m_yieldLaser);
233 m_fitTree->Branch<float>("yieldLaserErr", &m_yieldLaserErr);
234 m_fitTree->Branch<float>("timeExtra", &m_timeExtra);
235 m_fitTree->Branch<float>("sigmaExtra", &m_sigmaExtra);
236 m_fitTree->Branch<float>("nExtra", &m_nExtra);
237 m_fitTree->Branch<float>("alphaExtra", &m_alphaExtra);
238 m_fitTree->Branch<float>("yieldLaserExtra", &m_yieldLaserExtra);
239 m_fitTree->Branch<float>("timeBackground", &m_timeBackground);
240 m_fitTree->Branch<float>("sigmaBackground", &m_sigmaBackground);
241 m_fitTree->Branch<float>("yieldLaserBackground", &m_yieldLaserBackground);
242 m_fitTree->Branch<float>("fractionMC", &m_fractionMC);
243 m_fitTree->Branch<float>("deltaTMC", &m_deltaTMC);
244 m_fitTree->Branch<float>("peakTimeMC", &m_peakTimeMC);
245 m_fitTree->Branch<float>("firstPulserTime", &m_firstPulserTime);
246 m_fitTree->Branch<float>("firstPulserSigma", &m_firstPulserSigma);
247 m_fitTree->Branch<float>("secondPulserTime", &m_secondPulserTime);
248 m_fitTree->Branch<float>("secondPulserSigma", &m_secondPulserSigma);
249 m_fitTree->Branch<short>("fitStatus", &m_fitStatus);
250 m_fitTree->Branch<double>("width", &m_width);
251 m_fitTree->Branch<double>("amplitude", &m_amplitude);
252 m_fitTree->Branch<float>("chi2", &m_chi2);
253 m_fitTree->Branch<float>("rms", &m_rms);
254
255
256 if (m_isFitInAmplitudeBins) {
257 m_timewalkTree = new TTree("timewalkTree", "timewalkTree");
258 m_timewalkTree->Branch<float>("binLowerEdge", &m_binLowerEdge);
259 m_timewalkTree->Branch<float>("binUpperEdge", &m_binUpperEdge);
260 m_timewalkTree->Branch<short>("channel", &m_channel);
261 m_timewalkTree->Branch<short>("slot", &m_slot);
262 m_timewalkTree->Branch<short>("row", &m_row);
263 m_timewalkTree->Branch<short>("col", &m_col);
264 m_timewalkTree->Branch<short>("asic", &m_asic);
265 m_timewalkTree->Branch<short>("asicChannel", &m_asicChannel);
266 m_timewalkTree->Branch<short>("boardstack", &m_boardstack);
267 m_timewalkTree->Branch<float>("histoIntegral", &m_histoIntegral);
268 m_timewalkTree->Branch<float>("peakTime", &m_peakTime);
269 m_timewalkTree->Branch<float>("peakTimeErr", &m_peakTimeErr);
270 m_timewalkTree->Branch<float>("deltaT", &m_deltaT);
271 m_timewalkTree->Branch<float>("deltaTErr", &m_deltaTErr);
272 m_timewalkTree->Branch<float>("sigma", &m_sigma);
273 m_timewalkTree->Branch<float>("sigmaErr", &m_sigmaErr);
274 m_timewalkTree->Branch<float>("fraction", &m_fraction);
275 m_timewalkTree->Branch<float>("fractionErr", &m_fractionErr);
276 m_timewalkTree->Branch<float>("yieldLaser", &m_yieldLaser);
277 m_timewalkTree->Branch<float>("yieldLaserErr", &m_yieldLaserErr);
278 m_timewalkTree->Branch<float>("timeExtra", &m_timeExtra);
279 m_timewalkTree->Branch<float>("sigmaExtra", &m_sigmaExtra);
280 m_timewalkTree->Branch<float>("nExtra", &m_nExtra);
281 m_timewalkTree->Branch<float>("alphaExtra", &m_alphaExtra);
282 m_timewalkTree->Branch<float>("yieldLaserExtra", &m_yieldLaserExtra);
283 m_timewalkTree->Branch<float>("timeBackground", &m_timeBackground);
284 m_timewalkTree->Branch<float>("sigmaBackground", &m_sigmaBackground);
285 m_timewalkTree->Branch<float>("yieldLaserBackground", &m_yieldLaserBackground);
286 m_timewalkTree->Branch<float>("fractionMC", &m_fractionMC);
287 m_timewalkTree->Branch<float>("deltaTMC", &m_deltaTMC);
288 m_timewalkTree->Branch<float>("peakTimeMC", &m_peakTimeMC);
289 m_timewalkTree->Branch<float>("firstPulserTime", &m_firstPulserTime);
290 m_timewalkTree->Branch<float>("firstPulserSigma", &m_firstPulserSigma);
291 m_timewalkTree->Branch<float>("secondPulserTime", &m_secondPulserTime);
292 m_timewalkTree->Branch<float>("secondPulserSigma", &m_secondPulserSigma);
293 m_timewalkTree->Branch<short>("fitStatus", &m_fitStatus);
294 m_timewalkTree->Branch<double>("width", &m_width);
295 m_timewalkTree->Branch<double>("amplitude", &m_amplitude);
296 m_timewalkTree->Branch<float>("chi2", &m_chi2);
297 m_timewalkTree->Branch<float>("rms", &m_rms);
298 }
299
300 if (m_detectCrosstalk) {
301
302 // Create tree that stores candidate crosstalk events
303 m_crosstalkTree = new TTree("crosstalkTree", "Tree containing candidate crosstalks");
304 m_crosstalkTree->Branch<short>("sl0", &m_sl0);
305 m_crosstalkTree->Branch<short>("sl1", &m_sl1); // slot numbers (they will be the same, including them for checks)
306 m_crosstalkTree->Branch<short>("ch0", &m_ch0);
307 m_crosstalkTree->Branch<short>("ch1", &m_ch1); // channel numbers
308 m_crosstalkTree->Branch<float>("ht0", &m_ht0);
309 m_crosstalkTree->Branch<float>("ht1", &m_ht1); // hit times
310 m_crosstalkTree->Branch<float>("a0", &m_a0);
311 m_crosstalkTree->Branch<float>("a1", &m_a1); // amplitudes
312 m_crosstalkTree->Branch<float>("w0", &m_w0);
313 m_crosstalkTree->Branch<float>("w1", &m_w1); // widths
314 m_crosstalkTree->Branch<float>("q0", &m_q0);
315 m_crosstalkTree->Branch<float>("q1", &m_q1); // integrated charges
316 m_crosstalkTree->Branch<float>("f_q0", &m_f_q0); // fraction of charge on channel 0
317
318 // Create tree that stores fit results of hits without associated crosstalk
319 // Unlike the "vanilla" fitTree, this doesn't contain the results of the fits to the calibration pulses
320 m_fitTree_noXtalk = new TTree("fitTreeNoXTalk", "Fits to channels with no detected crosstalk");
321 m_fitTree_noXtalk->Branch<short>("channel", &m_channel);
322 m_fitTree_noXtalk->Branch<short>("slot", &m_slot);
323 m_fitTree_noXtalk->Branch<short>("row", &m_row);
324 m_fitTree_noXtalk->Branch<short>("col", &m_col);
325 m_fitTree_noXtalk->Branch<short>("asic", &m_asic);
326 m_fitTree_noXtalk->Branch<short>("asicChannel", &m_asicChannel);
327 m_fitTree_noXtalk->Branch<short>("boardstack", &m_boardstack);
328 m_fitTree_noXtalk->Branch<float>("peakTime", &m_peakTime);
329 m_fitTree_noXtalk->Branch<float>("peakTimeErr", &m_peakTimeErr);
330 m_fitTree_noXtalk->Branch<float>("deltaT", &m_deltaT);
331 m_fitTree_noXtalk->Branch<float>("deltaTErr", &m_deltaTErr);
332 m_fitTree_noXtalk->Branch<float>("sigma", &m_sigma);
333 m_fitTree_noXtalk->Branch<float>("sigmaErr", &m_sigmaErr);
334 m_fitTree_noXtalk->Branch<float>("fraction", &m_fraction);
335 m_fitTree_noXtalk->Branch<float>("fractionErr", &m_fractionErr);
336 m_fitTree_noXtalk->Branch<float>("yieldLaser", &m_yieldLaser);
337 m_fitTree_noXtalk->Branch<float>("yieldLaserErr", &m_yieldLaserErr);
338 m_fitTree_noXtalk->Branch<float>("timeExtra", &m_timeExtra);
339 m_fitTree_noXtalk->Branch<float>("sigmaExtra", &m_sigmaExtra);
340 m_fitTree_noXtalk->Branch<float>("nExtra", &m_nExtra);
341 m_fitTree_noXtalk->Branch<float>("alphaExtra", &m_alphaExtra);
342 m_fitTree_noXtalk->Branch<float>("yieldLaserExtra", &m_yieldLaserExtra);
343 m_fitTree_noXtalk->Branch<float>("timeBackground", &m_timeBackground);
344 m_fitTree_noXtalk->Branch<float>("sigmaBackground", &m_sigmaBackground);
345 m_fitTree_noXtalk->Branch<float>("yieldLaserBackground", &m_yieldLaserBackground);
346 m_fitTree_noXtalk->Branch<float>("fractionMC", &m_fractionMC);
347 m_fitTree_noXtalk->Branch<float>("deltaTMC", &m_deltaTMC);
348 m_fitTree_noXtalk->Branch<float>("peakTimeMC", &m_peakTimeMC);
349 m_fitTree_noXtalk->Branch<float>("firstPulserTime", &m_firstPulserTime);
350 m_fitTree_noXtalk->Branch<float>("firstPulserSigma", &m_firstPulserSigma);
351 m_fitTree_noXtalk->Branch<float>("secondPulserTime", &m_secondPulserTime);
352 m_fitTree_noXtalk->Branch<float>("secondPulserSigma", &m_secondPulserSigma);
353 m_fitTree_noXtalk->Branch<short>("fitStatus", &m_fitStatus);
354 m_fitTree_noXtalk->Branch<double>("width", &m_width);
355 m_fitTree_noXtalk->Branch<double>("amplitude", &m_amplitude);
356 m_fitTree_noXtalk->Branch<float>("chi2", &m_chi2);
357 m_fitTree_noXtalk->Branch<float>("rms", &m_rms);
358
359 }
360
361 return;
362}
363
364void Belle2::TOP::TOPLocalCalFitter::fitChannel(short iSlot, short iChannel, TH1* h_profile, bool inBins, double frac)
365{
366 // loads the TTS infos and the fit constraint for the given channel and slot
367 if (m_fitterMode == "monitoring")
368 m_treeConstraints->GetEntry(iChannel + 512 * iSlot);
369 else if (m_fitterMode == "calibration") // The MC-based constraint file has only slot 1 at the moment
370 m_treeConstraints->GetEntry(iChannel);
371
372 m_treeTTS->GetEntry(iChannel + 512 * iSlot);
373 // finds the maximum of the hit timing histogram and adjust the histogram range around it (3 ns window)
374 double maxpos = h_profile->GetBinCenter(h_profile->GetMaximumBin());
375 h_profile->GetXaxis()->SetRangeUser(maxpos - 1, maxpos + 2.);
376
377 // gets the histogram integral to give a starting value to the fitter
378 double integral = h_profile->Integral();
379
380 // creates the fit function
381 TF1 laser = TF1("laser", laserPDF, maxpos - 1, maxpos + 2., 16);
382
383 // par[0] = peakTime
384 laser.SetParameter(0, maxpos);
385 laser.SetParLimits(0, maxpos - 0.06, maxpos + 0.06);
386
387 // par[1] = sigma
388 laser.SetParameter(1, 0.1);
389 laser.SetParLimits(1, 0.05, 0.25);
390 if (m_fitterMode == "MC") {
391 laser.SetParameter(1, 0.02);
392 laser.SetParLimits(1, 0., 0.04);
393 }
394
395 // par[2] = fraction of the main peak respect to the total
396 laser.SetParameter(2, m_fractionConstraints);
397 laser.SetParLimits(2, 0.5, 1.);
398 if (inBins) {
399 laser.FixParameter(2, frac);
400 }
401
402 // par[3]= time difference between the main and secondary path. fixed to the MC value
403 laser.FixParameter(3, m_deltaTConstraints);
404
405 // This is an hack: in some channels the MC sees one peak only, while in the data there are clearly
406 // two well distinguished peaks. This will disappear if we'll ever get a better laser simulation.
407 if (m_deltaTConstraints > -0.001) {
408 laser.SetParameter(3, -0.3);
409 laser.SetParLimits(3, -0.4, -0.2);
410 }
411
412 // par[4] is the quadratic difference of the sigmas of the two TTS gaussians (tail - core)
413 laser.FixParameter(4, TMath::Sqrt(m_sigma2 * m_sigma2 - m_sigma1 * m_sigma1));
414 // par[5] is the position of the second TTS gaussian w/ respect to the first one
415 laser.FixParameter(5, m_mean2);
416 // par[6] is the relative contribution of the second TTS gaussian
417 laser.FixParameter(6, m_f1);
418 if (m_fitterMode == "MC")
419 laser.FixParameter(6, 0);
420
421 // par[7] is the PDF normalization, = integral*bin width
422 const double binw = h_profile->GetXaxis()->GetBinWidth(1);
423 laser.SetParameter(7, integral * binw);
424 laser.SetParLimits(7, 0.2 * integral * binw, 2.*integral * binw);
425
426 // par[8-10] are the relative position, the sigma and the integral of the extra peak
427 laser.SetParameter(8, 1.);
428 laser.SetParLimits(8, 0.3, 2.);
429 laser.SetParameter(9, 0.2);
430 laser.SetParLimits(9, 0.08, 1.);
431 laser.SetParameter(10, 0.1 * integral * binw);
432 laser.SetParLimits(10, 0., 0.2 * integral * binw);
433 // par[14-15] are the tail parameters of the crystal ball function used to describe the extra peak
434 laser.SetParameter(14, -2.);
435 laser.SetParameter(15, 2.);
436 laser.SetParLimits(15, 1.01, 20.);
437
438 // par[11-13] are relative position, sigma and integral of the broad gaussian added to better describe the tail at high times
439 laser.SetParameter(11, 1.);
440 laser.SetParLimits(11, 0.1, 5.);
441 laser.SetParameter(12, 0.8);
442 laser.SetParLimits(12, 0., 5.);
443 laser.SetParameter(13, 0.01 * integral * binw);
444 laser.SetParLimits(13, 0., 0.2 * integral * binw);
445
446 // if it's a monitoring fit, fix a buch more parameters.
447 if (m_fitterMode == "monitoring") {
448 laser.FixParameter(2, m_fractionConstraints);
449 laser.FixParameter(3, m_deltaTConstraints);
450 laser.FixParameter(8, m_timeExtraConstraints);
451 laser.FixParameter(9, m_sigmaExtraConstraints);
452 laser.FixParameter(14, m_alphaExtraConstraints);
453 laser.FixParameter(15, m_nExtraConstraints);
454 laser.FixParameter(11, m_timeBackgroundConstraints);
455 laser.FixParameter(12, m_sigmaBackgroundConstraints);
456 }
457
458 // if it's a MC fit, fix a buch more parameters.
459 if (m_fitterMode == "MC") {
460 laser.SetParameter(2, 0.8);
461 laser.SetParLimits(2, 0., 1.);
462 laser.SetParameter(3, -0.1);
463 laser.SetParLimits(3, -0.4, -0.);
464 // The following are just random reasonable number, only to pin-point the tail components to some value and remove them form the fit
465 laser.FixParameter(8, 0);
466 laser.FixParameter(9, 0.1);
467 laser.FixParameter(14, -2.);
468 laser.FixParameter(15, 2);
469 laser.FixParameter(11, 1.);
470 laser.FixParameter(12, 0.1);
471 laser.FixParameter(13, 0.);
472 laser.FixParameter(10, 0.);
473 }
474
475 // make the plot of the fit function nice setting 2000 sampling points
476 laser.SetNpx(2000);
477
478 // do the fit!
479 h_profile->Fit("laser", "R L Q");
480
481 // Add by hand the different fit components to the histogram, mostly for debugging/presentation purposes
482 TF1* peak1 = new TF1("peak1", laserPDF, maxpos - 1, maxpos + 2., 16);
483 TF1* peak2 = new TF1("peak2", laserPDF, maxpos - 1, maxpos + 2., 16);
484 TF1* extra = new TF1("extra", laserPDF, maxpos - 1, maxpos + 2., 16);
485 TF1* background = new TF1("background", laserPDF, maxpos - 1, maxpos + 2., 16);
486 for (int iPar = 0; iPar < 16; iPar++) {
487 peak1->FixParameter(iPar, laser.GetParameter(iPar));
488 peak2->FixParameter(iPar, laser.GetParameter(iPar));
489 extra->FixParameter(iPar, laser.GetParameter(iPar));
490 background->FixParameter(iPar, laser.GetParameter(iPar));
491 }
492 peak1->FixParameter(2, 0.);
493 peak1->FixParameter(7, (1 - laser.GetParameter(2))*laser.GetParameter(7));
494 peak1->FixParameter(10, 0.);
495 peak1->FixParameter(13, 0.);
496 peak2->FixParameter(2, 1.);
497 peak2->FixParameter(7, laser.GetParameter(2)*laser.GetParameter(7));
498 peak2->FixParameter(10, 0.);
499 peak2->FixParameter(13, 0.);
500 extra->FixParameter(7, 0.);
501 extra->FixParameter(13, 0.);
502 background->FixParameter(7, 0.);
503 background->FixParameter(10, 0.);
504
505 h_profile->GetListOfFunctions()->Add(peak1);
506 h_profile->GetListOfFunctions()->Add(peak2);
507 h_profile->GetListOfFunctions()->Add(extra);
508 h_profile->GetListOfFunctions()->Add(background);
509
510 // save the results in the variables linked to the tree branches
511 m_channel = iChannel;
512 setHardwareIdentifiers(m_channel);
513 m_row = rowOf(iSlot, iChannel);
514 m_col = colOf(iSlot, iChannel);
515 m_slot = iSlot + 1;
516 m_peakTime = laser.GetParameter(0);
517 m_peakTimeErr = laser.GetParError(0);
518 m_deltaT = laser.GetParameter(3);
519 m_deltaTErr = laser.GetParError(3);
520 m_sigma = laser.GetParameter(1);
521 m_sigmaErr = laser.GetParError(1);
522 m_fraction = laser.GetParameter(2);
523 m_fractionErr = laser.GetParError(2);
524 m_yieldLaser = laser.GetParameter(7) / binw;
525 m_yieldLaserErr = laser.GetParError(7) / binw;
526 m_timeExtra = laser.GetParameter(8);
527 m_sigmaExtra = laser.GetParameter(9);
528 m_yieldLaserExtra = laser.GetParameter(10) / binw;
529 m_alphaExtra = laser.GetParameter(14);
530 m_nExtra = laser.GetParameter(15);
531 m_timeBackground = laser.GetParameter(11);
532 m_sigmaBackground = laser.GetParameter(12);
533 m_yieldLaserBackground = laser.GetParameter(13) / binw;
534 m_chi2 = laser.GetChisquare() / laser.GetNDF();
535
536 // copy some MC information to the output tree
537 m_fractionMC = m_fractionConstraints;
538 m_deltaTMC = m_deltaTConstraints;
539 m_peakTimeMC = m_peakTimeConstraints;
540
541 return;
542}
543
544void Belle2::TOP::TOPLocalCalFitter::fitPulser(TH1* h_profileFirstPulser, TH1* h_profileSecondPulser)
545{
546 float maxpos = h_profileFirstPulser->GetBinCenter(h_profileFirstPulser->GetMaximumBin());
547 h_profileFirstPulser->GetXaxis()->SetRangeUser(maxpos - 1, maxpos + 1.);
548 if (h_profileFirstPulser->Integral() > 1000) {
549 TF1 pulser1 = TF1("pulser1", "[0]*TMath::Gaus(x, [1], [2], kTRUE)", maxpos - 1, maxpos + 1.);
550 pulser1.SetParameter(0, 1.);
551 pulser1.SetParameter(1, maxpos);
552 pulser1.SetParameter(2, 0.05);
553 h_profileFirstPulser->Fit("pulser1", "R Q");
554 m_firstPulserTime = pulser1.GetParameter(1);
555 m_firstPulserSigma = pulser1.GetParameter(2);
556 h_profileFirstPulser->Write();
557 } else {
558 m_firstPulserTime = -999;
559 m_firstPulserSigma = -999;
560 }
561
562 maxpos = h_profileSecondPulser->GetBinCenter(h_profileSecondPulser->GetMaximumBin());
563 h_profileSecondPulser->GetXaxis()->SetRangeUser(maxpos - 1, maxpos + 1.);
564 if (h_profileSecondPulser->Integral() > 1000) {
565 TF1 pulser2 = TF1("pulser2", "[0]*TMath::Gaus(x, [1], [2], kTRUE)", maxpos - 1, maxpos + 1.);
566 pulser2.SetParameter(0, 1.);
567 pulser2.SetParameter(1, maxpos);
568 pulser2.SetParameter(2, 0.05);
569 h_profileSecondPulser->Fit("pulser2", "R Q");
570 m_secondPulserTime = pulser2.GetParameter(1);
571 m_secondPulserSigma = pulser2.GetParameter(2);
572 h_profileSecondPulser->Write();
573 } else {
574 m_secondPulserTime = -999;
575 m_secondPulserSigma = -999;
576 }
577 return;
578}
579
580void Belle2::TOP::TOPLocalCalFitter::determineFitStatus()
581{
582 if (m_chi2 < 4 && m_sigma < 0.2 && m_yieldLaser > 1000) {
583 m_fitStatus = 0;
584 } else {
585 m_fitStatus = 1;
586 }
587 return;
588}
589
590void Belle2::TOP::TOPLocalCalFitter::calculateChannelT0()
591{
592 Long64_t nEntries = m_fitTree->GetEntries();
593 if (nEntries != 8192) {
594 B2ERROR("fitTree does not contain an entry with a fit result for each channel. Found " << nEntries <<
595 " instead of 8192. Perhaps you tried to run the commonT0 calculation before finishing the fitting?");
596 return;
597 }
598
599 // Create and fill the TOPCalChannelT0 object
600 auto* channelT0 = new TOPCalChannelT0();
601 short nCal[16] = {0};
602 for (Long64_t i = 0; i < nEntries; i++) {
603 m_fitTree->GetEntry(i);
604 channelT0->setT0(m_slot, m_channel, m_peakTime - m_peakTimeMC, m_peakTimeErr);
605 if (m_fitStatus == 0) {
606 nCal[m_slot - 1]++;
607 } else {
608 channelT0->setUnusable(m_slot, m_channel);
609 }
610 }
611
612 // Normalize the constants
613 channelT0->suppressAverage();
614
615 // create the localDB
616 saveCalibration(channelT0);
617
618 short nCalTot = 0;
619 B2INFO("Summary: ");
620 for (int iSlot = 1; iSlot < 17; iSlot++) {
621 B2INFO("--> Number of calibrated channels on Slot " << iSlot << " : " << nCal[iSlot - 1] << "/512");
622 B2INFO("--> Cal on ch 1, 256 and 511: " << channelT0->getT0(iSlot, 0) << ", " << channelT0->getT0(iSlot,
623 257) << ", " << channelT0->getT0(iSlot, 511));
624 nCalTot += nCal[iSlot - 1];
625 }
626
627 B2RESULT("Channel T0 calibration constants imported to database, calibrated channels: " << nCalTot << "/ 8192");
628
629 // Loop again on the output tree to save the constants there too, adding two more branches.
630 TBranch* channelT0Branch = m_fitTree->Branch<float>("channelT0", &m_channelT0);
631 TBranch* channelT0ErrBranch = m_fitTree->Branch<float>("channelT0Err", &m_channelT0Err);
632
633 for (int i = 0; i < nEntries; i++) {
634 m_fitTree->GetEntry(i);
635 m_channelT0 = channelT0->getT0(m_slot, m_channel);
636 m_channelT0Err = channelT0->getT0Error(m_slot, m_channel);
637 channelT0Branch->Fill();
638 channelT0ErrBranch->Fill();
639 }
640
641 return;
642
643}
644
645// Main fitter function
646Belle2::CalibrationAlgorithm::EResult Belle2::TOP::TOPLocalCalFitter::calibrate()
647{
648
649 gROOT->SetBatch();
650
651 // Load MC constraints
652 loadMCInfoTrees();
653
654 // Prepare output
655 setupOutputTreeAndFile();
656
657 // Load the tree with the hits (output of TOPLaserCalibratorCollector)
658 auto hitTree = getObjectPtr<TTree>("hitTree");
659 int event;
660 float amplitude, width, hitTime;
661 short channel, slot; //, row, col;
662 bool refTimeValid;
663 hitTree->SetBranchAddress("event", &event);
664 hitTree->SetBranchAddress("amplitude", &amplitude);
665 hitTree->SetBranchAddress("width", &width);
666 hitTree->SetBranchAddress("hitTime", &hitTime);
667 hitTree->SetBranchAddress("channel", &channel);
668 //hitTree->SetBranchAddress("row", &row);
669 //hitTree->SetBranchAddress("col", &col);
670 hitTree->SetBranchAddress("slot", &slot);
671 hitTree->SetBranchAddress("refTimeValid", &refTimeValid);
672
673 // Prepare histogram to save interesting features of each channel
674 TH2F* h_hitTime = new TH2F("h_hitTime", " ", 512 * 16, 0., 512 * 16, 22000, -70, 40.); // 5 ps bins
675 TH2F* h_amplitude2D = new TH2F("h_amplitude", " ", 512 * 16, 0., 512 * 16, 600, 0, 2200.);
676 TH2F* h_width2D = new TH2F("h_width", " ", 512 * 16, 0., 512 * 16, 1000, 0, 2.);
677
678 // Prepare vector of hitTime vs channel histograms for fits in amplitude bins
679 // (attempt to speed things up looping over the hitTree only once).
680
681 std::vector<TH2F*> h_hitTimeLaserHistos = {};
682 for (int iLowerEdge = 0; iLowerEdge < (int)m_binEdges.size() - 1; iLowerEdge++) {
683 TH2F* h_hitTimeLaser = new TH2F(("h_hitTimeLaser_" + std::to_string(iLowerEdge + 1)).c_str(), " ",
684 512 * 16, 0., 512 * 16, 14000, -70, 0.); // 5 ps bins
685 h_hitTimeLaserHistos.push_back(h_hitTimeLaser);
686 }
687
688 // Per-event accumulator for crosstalk logic
689 struct Hit {
690 short slot;
691 short ch;
692 float t; // hitTime
693 float a; // amplitude
694 float w; // width
695 float q; // charge proxy = conv * a * w
696 };
697 std::vector<Hit> evtHits;
698 //evtHits.reserve(64); // typical multiplicity in laser runs?
699
700 // Track current event id we are accumulating
701 int prev_evt = std::numeric_limits<int>::min();
702
703 // Conversion factor to compute approximate integrated charge (ADC·ns)
704 const float conv = 1.f / (0.3989f * 2.35f);
705
706 // Crosstalk tunables (consider moving to members with setters)
707 const float fracMin = 0.25f; // f_q0 < fracMin or > 1-fracMin => crosstalk
708 const float dtMax = 0.30f; // ns, time-coincidence window for pair
709 const float epsQ = 1e-6f; // guard against zero-sum charges
710
711 // Prepares histogram to store features of each channel (if no crosstalk detected)
712 // which will be fit later
713 TH2F* h_hitTime_noXtalk = new TH2F("h_hitTime_noXtalk", " ", 512 * 16, 0., 512 * 16, 22000, -70, 40.); // 5 ps bins
714 TH2F* h_amplitude2D_noXtalk = new TH2F("h_amplitude2D_noXtalk", " ", 512 * 16, 0., 512 * 16, 600, 0, 2200.);
715 TH2F* h_width2D_noXtalk = new TH2F("h_width2D_noXtalk", " ", 512 * 16, 0., 512 * 16, 1000, 0, 2.);
716
717 // Get number of entries in hitTree
718 Long64_t nhits = hitTree->GetEntries();
719 const Long64_t step = std::max<Long64_t>(1, nhits / 100); // for progress bar
720
721 // Loop on all hits to retrieve information
722 for (Long64_t i = 0; i < nhits; i++) {
723
724 // Print percentage of completion
725 if (i % step == 0) {
726 std::cout << "Processing hit " << i << " of " << nhits << " ("
727 << std::setprecision(3) << (100. * i) / nhits << " %)" << std::endl;
728 }
729
730 // Process entry
731 hitTree->GetEntry(i);
732
733 // Fill hit time histograms for each bin of pulse heigth (if activated)
734 if (m_isFitInAmplitudeBins) {
735 auto it = std::lower_bound(m_binEdges.cbegin(), m_binEdges.cend(), amplitude); // std::vector iterator
736 int iLowerEdge = std::distance(m_binEdges.cbegin(), it) - 1;
737 if (iLowerEdge >= 0 && iLowerEdge < static_cast<int>(m_binEdges.size()) - 1 && refTimeValid)
738 h_hitTimeLaserHistos[iLowerEdge]->Fill(channel + (slot - 1) * 512, hitTime);
739 }
740
741 // Check if pulse is at least 80 ADC and has valid reference time (suppress noise)
742 if (amplitude > 80. && refTimeValid) {
743
744 // Fill the hitTime vs channel histogram
745 h_hitTime->Fill(channel + (slot - 1) * 512, hitTime);
746
747 // If entry is found with -65 < hitTime < -10 (laser pulse), fill amplitude and width histograms
748 if ((hitTime > -65) && (hitTime < -10)) { // use logical &&
749
750 // Fill amplitude and width histograms
751 h_amplitude2D->Fill(channel + (slot - 1) * 512, amplitude);
752 h_width2D->Fill(channel + (slot - 1) * 512, width);
753
754 // Crosstalk finder algorithm (only on laser pulses)
755 if (m_detectCrosstalk) {
756
757 const int curr_evt = event;
758
759 // If this entry belongs to a NEW event, process the accumulated previous event
760 if (curr_evt != prev_evt && !evtHits.empty()) {
761
762 // ---- step 1: find crosstalk pairs and mark involved hits ----
763 std::vector<char> isXtalk(evtHits.size(), 0);
764 for (size_t ii = 0; ii + 1 < evtHits.size(); ++ii) {
765 const auto& hi = evtHits[ii];
766 for (size_t jj = ii + 1; jj < evtHits.size(); ++jj) {
767 const auto& hj = evtHits[jj];
768
769 // same slot and neighboring pixels
770 if (hi.slot != hj.slot) continue;
771 if (!areNeighbors(hi.slot - 1, hi.ch, hj.ch)) continue;
772
773 // near-coincident in time (laser)
774 if (std::fabs(hi.t - hj.t) > dtMax) continue;
775
776 // robust fraction of shared charge
777 const float qsum = hi.q + hj.q;
778 if (qsum <= epsQ) continue;
779 const float f_q0 = hi.q / qsum;
780
781 if (f_q0 < fracMin || f_q0 > (1.f - fracMin)) {
782 isXtalk[ii] = 1;
783 isXtalk[jj] = 1;
784
785 // save diagnostics once per pair
786 if (m_crosstalkTree) {
787 m_sl0 = hi.slot; m_sl1 = hj.slot;
788 m_ch0 = hi.ch; m_ch1 = hj.ch;
789 m_ht0 = hi.t; m_ht1 = hj.t;
790 m_a0 = hi.a; m_a1 = hj.a;
791 m_w0 = hi.w; m_w1 = hj.w;
792 m_q0 = hi.q; m_q1 = hj.q;
793 m_f_q0 = f_q0;
794 m_crosstalkTree->Fill();
795 }
796 }
797 }
798 }
799
800 // ---- step 2: fill "no-crosstalk" histograms exactly once per clean hit ----
801 for (size_t kk = 0; kk < evtHits.size(); ++kk) {
802 if (isXtalk[kk]) continue;
803 const auto& h = evtHits[kk];
804 const int gch = h.ch + (h.slot - 1) * 512;
805 h_hitTime_noXtalk->Fill(gch, h.t);
806 h_amplitude2D_noXtalk->Fill(gch, h.a);
807 h_width2D_noXtalk->Fill(gch, h.w);
808 }
809
810 // clear for next event
811 evtHits.clear();
812 }
813 // accumulate the current hit for the (possibly new) event
814 evtHits.push_back(Hit{slot, channel, hitTime, amplitude, width, conv* amplitude * width});
815
816 // update current event id
817 prev_evt = curr_evt;
818 }
819 }
820 }
821 }
822
823 // Final flush for the last accumulated event (if any)
824 if (m_detectCrosstalk && !evtHits.empty()) {
825 std::vector<char> isXtalk(evtHits.size(), 0);
826 for (size_t ll = 0; ll + 1 < evtHits.size(); ++ll) {
827 const auto& hi = evtHits[ll];
828 for (size_t mm = ll + 1; mm < evtHits.size(); ++mm) {
829 const auto& hj = evtHits[mm];
830 if (hi.slot != hj.slot) continue;
831 if (!areNeighbors(hi.slot - 1, hi.ch, hj.ch)) continue;
832 if (std::fabs(hi.t - hj.t) > dtMax) continue;
833 const float qsum = hi.q + hj.q;
834 if (qsum <= epsQ) continue;
835 const float f_q0 = hi.q / qsum;
836 if (f_q0 < fracMin || f_q0 > (1.f - fracMin)) {
837 isXtalk[ll] = 1;
838 isXtalk[mm] = 1;
839 if (m_crosstalkTree) {
840 m_sl0 = hi.slot; m_sl1 = hj.slot;
841 m_ch0 = hi.ch; m_ch1 = hj.ch;
842 m_ht0 = hi.t; m_ht1 = hj.t;
843 m_a0 = hi.a; m_a1 = hj.a;
844 m_w0 = hi.w; m_w1 = hj.w;
845 m_q0 = hi.q; m_q1 = hj.q;
846 m_f_q0 = f_q0;
847 m_crosstalkTree->Fill();
848 }
849 }
850 }
851 }
852 for (size_t nn = 0; nn < evtHits.size(); ++nn) {
853 if (isXtalk[nn]) continue;
854 const auto& h = evtHits[nn];
855 const int gch = h.ch + (h.slot - 1) * 512;
856 h_hitTime_noXtalk->Fill(gch, h.t);
857 h_amplitude2D_noXtalk->Fill(gch, h.a);
858 h_width2D_noXtalk->Fill(gch, h.w);
859 }
860 evtHits.clear();
861 }
862
863 // Save tree filled with candidate crosstalks => useful for further studies
864 if (m_detectCrosstalk) {
865 std::cout << "Writing crosstalkTree (candidate crosstalk channel pairs) to output file" << std::endl;
866 m_histFile->cd();
867 m_crosstalkTree->Write();
868 }
869
870 // Write hitTime histograms to file
871 m_histFile->cd();
872 h_hitTime->Write();
873
874 // After filling hitTime, amplitude, width vs. channel histograms,
875 // loop on each slot and channel to perform channelT0 fits on hitTime profiles
876 for (short iSlot = 0; iSlot < 16; iSlot++) {
877 std::cout << "fitting slot " << iSlot + 1 << std::endl;
878 for (short iChannel = 0; iChannel < 512; iChannel++) {
879
880 // Project to 1-d hitTime distribution
881 TH1D* h_profile = h_hitTime->ProjectionY(
882 ("profile_" + std::to_string(iSlot + 1) + "_" + std::to_string(iChannel)).c_str(),
883 iSlot * 512 + iChannel + 1, iSlot * 512 + iChannel + 1
884 );
885
886 // Set hitTime range based on fitter mode
887 if (m_fitterMode == "MC")
888 //h_profile->GetXaxis()->SetRangeUser(-10, -10);
889 h_profile->GetXaxis()->SetRangeUser(-65, -1);
890 else // if you will even change the limits, make sure not to include the h_hitTime overflow bins in this range
891 h_profile->GetXaxis()->SetRangeUser(-65, -5);
892
893 // Run fit and determine status
894 fitChannel(iSlot, iChannel, h_profile);
895 determineFitStatus();
896
897 // Now let's fit the pulser
898 TH1D* h_profileFirstPulser = h_hitTime->ProjectionY(
899 ("profileFirstPulser_" + std::to_string(iSlot + 1) + "_" + std::to_string(
900 iChannel)).c_str(), iSlot * 512 + iChannel + 1, iSlot * 512 + iChannel + 1
901 );
902 TH1D* h_profileSecondPulser = h_hitTime->ProjectionY(
903 ("profileSecondPulser_" + std::to_string(iSlot + 1) + "_" + std::to_string(
904 iChannel)).c_str(), iSlot * 512 + iChannel + 1, iSlot * 512 + iChannel + 1
905 );
906 h_profileFirstPulser->GetXaxis()->SetRangeUser(-10, 10);
907 h_profileSecondPulser->GetXaxis()->SetRangeUser(10, 40);
908 fitPulser(h_profileFirstPulser, h_profileSecondPulser);
909
910
911 // Get pulse heigth [ADC] and width [ns]
912 TH1D* h_amplitude = h_amplitude2D->ProjectionY(
913 ("AmpProfile_" + std::to_string(iSlot + 1) + "_" + std::to_string(iChannel)).c_str(),
914 iSlot * 512 + iChannel + 1, iSlot * 512 + iChannel + 1
915 );
916 TH1D* h_width = h_width2D->ProjectionY(
917 ("WidthProfile_" + std::to_string(iSlot + 1) + "_" + std::to_string(iChannel)).c_str(),
918 iSlot * 512 + iChannel + 1, iSlot * 512 + iChannel + 1
919 );
920
921 // Set values for pulse amplitude and width (mean)
922 Double_t q = 0.5; // quantile for median
923 h_amplitude->GetQuantiles(1, &m_amplitude, &q);
924 h_width->GetQuantiles(1, &m_width, &q);
925
926 m_fitTree->Fill();
927 h_profile->Write();
928 h_profileFirstPulser->Write();
929 h_profileSecondPulser->Write();
930 h_amplitude->Write();
931 h_width->Write();
932
933 // Free memory: TH2D::ProjectionY allocates memory dynamically
934 delete h_profile;
935 delete h_profileFirstPulser;
936 delete h_profileSecondPulser;
937 delete h_amplitude;
938 delete h_width;
939
940 }
941
942 // Write hitTime histogram to output tree
943 h_hitTime->Write();
944
945 }
946
947 // Compute calibration constants (w/ error) and add corresponding branches to output fitTree
948 calculateChannelT0();
949
950 // Write fit results to tree
951 m_fitTree->Write();
952
953 // ChannelT0 fits in bins of pulse heigth
954 if (m_isFitInAmplitudeBins) {
955
956 std::cout << "Fitting in bins of pulse heigth" << std::endl;
957
958 for (short iSlot = 0; iSlot < 16; iSlot++) {
959 std::cout << " Fitting slot " << iSlot + 1 << std::endl;
960 for (short iChannel = 0; iChannel < 512; iChannel++) {
961
962 // The fraction parameter should not depend on amplitude ==> let's fix it
963 // to the value we get from the fit integrated on all amplitudes
964 TH1D* h_profile_full = h_hitTime->ProjectionY(
965 ("profile_" + std::to_string(iSlot + 1) + "_" + std::to_string(iChannel)).c_str(),
966 iSlot * 512 + iChannel + 1, iSlot * 512 + iChannel + 1
967 );
968 fitChannel(iSlot, iChannel, h_profile_full);
969 float ff = m_fraction;
970 // Fit done, free the memory
971 delete h_profile_full;
972
973 // Loop on the amplitude bins
974 for (int iLowerEdge = 0; iLowerEdge < (int)m_binEdges.size() - 1; iLowerEdge++) {
975
976 // Get current bin edges
977 m_binLowerEdge = m_binEdges[iLowerEdge];
978 m_binUpperEdge = m_binEdges[iLowerEdge + 1];
979 std::cout << "Fitting the amplitude interval (" << m_binLowerEdge << ", " << m_binUpperEdge << " )" << std::endl;
980
981 // Get profile for current amplitude bin
982 TH1D* h_profile = h_hitTimeLaserHistos[iLowerEdge]->ProjectionY(
983 ("profile_" + std::to_string(iSlot + 1) + "_" + std::to_string(
984 iChannel) + "_" + std::to_string(iLowerEdge)).c_str(),
985 iSlot * 512 + iChannel + 1, iSlot * 512 + iChannel + 1
986 );
987 // Set range of fit based on fitter mode
988 if (m_fitterMode == "MC")
989 h_profile->GetXaxis()->SetRangeUser(-10, -10);
990 else // if you will even change it, make sure not to include the h_hitTime overflow bins in this range
991 h_profile->GetXaxis()->SetRangeUser(-65, -5);
992
993
994 // Fit the hitTime distribution
995 fitChannel(iSlot, iChannel, h_profile, true, ff);
996 m_histoIntegral = h_profile->Integral();
997 determineFitStatus();
998
999 // Get amplitude and width profiles in order to get median amp and width of the channel
1000 TH1D* h_amplitude = h_amplitude2D->ProjectionY(
1001 ("AmpProfile_" + std::to_string(iSlot + 1) +
1002 "_" + std::to_string(iChannel) +
1003 "_" + std::to_string(iLowerEdge)).c_str(),
1004 iSlot * 512 + iChannel + 1, iSlot * 512 + iChannel + 1
1005 );
1006 TH1D* h_width = h_width2D->ProjectionY(
1007 ("WidthProfile_" + std::to_string(iSlot + 1) +
1008 "_" + std::to_string(iChannel) +
1009 "_" + std::to_string(iLowerEdge)).c_str(),
1010 iSlot * 512 + iChannel + 1, iSlot * 512 + iChannel + 1
1011 );
1012
1013 // Measure the *median* amplitude and widths => less sensitive to outliers
1014 Double_t q = 0.5; // quantile for median
1015 h_amplitude->GetQuantiles(1, &m_amplitude, &q);
1016 h_width->GetQuantiles(1, &m_width, &q);
1017
1018 // Fill tree and write histograms to output file
1019 m_timewalkTree->Fill();
1020 h_profile->Write();
1021 h_amplitude->Write();
1022 h_width->Write();
1023
1024 // Try to avoid memory leaks
1025 delete h_profile;
1026 delete h_amplitude;
1027 delete h_width;
1028
1029 }
1030 }
1031 }
1032 m_timewalkTree->Write();
1033 }
1034
1035 // Run fits on channels that didn't have cross-talk
1036 if (m_detectCrosstalk) {
1037 std::cout << "Fitting channels with no crosstalk detected" << std::endl;
1038 for (short iSlot = 0; iSlot < 16; iSlot++) {
1039 std::cout << " Fitting slot " << iSlot + 1 << std::endl;
1040 for (short iChannel = 0; iChannel < 512; iChannel++) {
1041 // Project to 1-d hitTime distribution
1042 TH1D* h_profile = h_hitTime_noXtalk->ProjectionY(
1043 ("profile_" + std::to_string(iSlot + 1) + "_" + std::to_string(iChannel)).c_str(),
1044 iSlot * 512 + iChannel + 1, iSlot * 512 + iChannel + 1
1045 );
1046
1047 // Set hitTime range based on fitter mode
1048 if (m_fitterMode == "MC")
1049 //h_profile->GetXaxis()->SetRangeUser(-10, -10);
1050 h_profile->GetXaxis()->SetRangeUser(-65, -1);
1051 else // if you will even change the limits, make sure not to include the h_hitTime overflow bins in this range
1052 h_profile->GetXaxis()->SetRangeUser(-65, -5);
1053
1054 // Run fit and determine status
1055 fitChannel(iSlot, iChannel, h_profile);
1056 determineFitStatus();
1057
1058 // Get pulse heigth [ADC] and width [ns]
1059 TH1D* h_amplitude = h_amplitude2D_noXtalk->ProjectionY(
1060 ("AmpProfile_" + std::to_string(iSlot + 1) + "_" + std::to_string(iChannel)).c_str(),
1061 iSlot * 512 + iChannel + 1, iSlot * 512 + iChannel + 1
1062 );
1063 TH1D* h_width = h_width2D_noXtalk->ProjectionY(
1064 ("WidthProfile_" + std::to_string(iSlot + 1) + "_" + std::to_string(iChannel)).c_str(),
1065 iSlot * 512 + iChannel + 1, iSlot * 512 + iChannel + 1
1066 );
1067
1068 // Set values for pulse amplitude and width (mean)
1069 Double_t q = 0.5; // quantile for median
1070 h_amplitude->GetQuantiles(1, &m_amplitude, &q);
1071 h_width->GetQuantiles(1, &m_width, &q);
1072
1073 m_fitTree_noXtalk->Fill();
1074 h_profile->Write();
1075 h_amplitude->Write();
1076 h_width->Write();
1077
1078 // Free memory: TH2D::ProjectionY allocates memory dynamically
1079 delete h_profile;
1080 delete h_amplitude;
1081 delete h_width;
1082 }
1083 }
1084 // Write fitTree with no crosstalk (after computing calibration constants)
1085 calculateChannelT0();
1086 m_fitTree_noXtalk->Write();
1087 }
1088
1089 m_histFile->Close();
1090
1091 return c_OK;
1092}
Belle2::CalibrationAlgorithm::saveCalibration
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
Definition CalibrationAlgorithm.cc:297
Belle2::CalibrationAlgorithm::setDescription
void setDescription(const std::string &description)
Set algorithm description (in constructor)
Definition CalibrationAlgorithm.h:321
Belle2::CalibrationAlgorithm::EResult
EResult
The result of calibration.
Definition CalibrationAlgorithm.h:40
Belle2::CalibrationAlgorithm::c_OK
@ c_OK
Finished successfully =0 in Python.
Definition CalibrationAlgorithm.h:41
Belle2::CalibrationAlgorithm::CalibrationAlgorithm
CalibrationAlgorithm(const std::string &collectorModuleName)
Constructor - sets the prefix for collected objects (won't be accesses until execute(....
Definition CalibrationAlgorithm.h:140
Belle2::TOPCalChannelT0
Channel T0 calibration constants for all 512 channels of 16 modules.
Definition TOPCalChannelT0.h:23
Belle2::TOP::TOPLocalCalFitter::m_fitStatus
short m_fitStatus
Fit quality flag, propagated to the constants.
Definition TOPLocalCalFitter.h:255
Belle2::TOP::TOPLocalCalFitter::colOf
short colOf(short slot, short ch) const noexcept
Column index for (slot,channel), or -1 if out of bounds.
Definition TOPLocalCalFitter.h:271
Belle2::TOP::TOPLocalCalFitter::m_asicChannel
short m_asicChannel
ASIC channel number (0-7)
Definition TOPLocalCalFitter.h:206
Belle2::TOP::TOPLocalCalFitter::m_crosstalkTree
TTree * m_crosstalkTree
Output tree for crosstalk candidates.
Definition TOPLocalCalFitter.h:159
Belle2::TOP::TOPLocalCalFitter::m_yieldLaserErr
float m_yieldLaserErr
Statistical error on yield.
Definition TOPLocalCalFitter.h:220
Belle2::TOP::TOPLocalCalFitter::m_binLowerEdge
float m_binLowerEdge
Lower edge of the amplitude bin in which this fit is performed.
Definition TOPLocalCalFitter.h:199
Belle2::TOP::TOPLocalCalFitter::m_nExtraConstraints
float m_nExtraConstraints
parameter n of the tail of the extra peak
Definition TOPLocalCalFitter.h:193
Belle2::TOP::TOPLocalCalFitter::m_chi2
float m_chi2
Reduced chi2 of the fit.
Definition TOPLocalCalFitter.h:238
Belle2::TOP::TOPLocalCalFitter::m_timeExtraConstraints
float m_timeExtraConstraints
Position of the gaussian used to describe the extra peak on the timing distribution tail.
Definition TOPLocalCalFitter.h:190
Belle2::TOP::TOPLocalCalFitter::fitChannel
void fitChannel(short slot, short channel, TH1 *h)
Fits the laser light on one channel.
Definition TOPLocalCalFitter.cc:29
Belle2::TOP::TOPLocalCalFitter::m_f1
float m_f1
Fraction of the first gaussian on the TTS parametrization.
Definition TOPLocalCalFitter.h:181
Belle2::TOP::TOPLocalCalFitter::m_fraction
float m_fraction
Fraction of events in the secondary peak.
Definition TOPLocalCalFitter.h:212
Belle2::TOP::TOPLocalCalFitter::m_sigmaExtra
float m_sigmaExtra
Gaussian sigma of the extra peak in the timing tail.
Definition TOPLocalCalFitter.h:223
Belle2::TOP::TOPLocalCalFitter::m_secondPulserSigma
float m_secondPulserSigma
Time resolution from the fit of the first electronic pulse, from a Gaussian fit.
Definition TOPLocalCalFitter.h:253
Belle2::TOP::TOPLocalCalFitter::m_yieldLaserBackground
float m_yieldLaserBackground
Integral of the background gaussian.
Definition TOPLocalCalFitter.h:230
Belle2::TOP::TOPLocalCalFitter::m_sigma
float m_sigma
Gaussian time resolution, fitted.
Definition TOPLocalCalFitter.h:211
Belle2::TOP::TOPLocalCalFitter::m_peakTimeMC
float m_peakTimeMC
Time of the main peak in the MC simulation, i.e.
Definition TOPLocalCalFitter.h:234
Belle2::TOP::TOPLocalCalFitter::m_output
std::string m_output
Name of the output file.
Definition TOPLocalCalFitter.h:138
Belle2::TOP::TOPLocalCalFitter::m_q0
float m_q0
Integrated charge for channel 0 in pair.
Definition TOPLocalCalFitter.h:173
Belle2::TOP::TOPLocalCalFitter::m_timeBackground
float m_timeBackground
Position of the gaussian used to describe the background, w/ respect to peakTime.
Definition TOPLocalCalFitter.h:228
Belle2::TOP::TOPLocalCalFitter::loadMCInfoTrees
void loadMCInfoTrees()
loads the TTS parameters and the MC truth info
Definition TOPLocalCalFitter.cc:173
Belle2::TOP::TOPLocalCalFitter::m_deltaTMC
float m_deltaTMC
Time difference between the main peak and the secondary peak in the MC simulation.
Definition TOPLocalCalFitter.h:233
Belle2::TOP::TOPLocalCalFitter::m_peakTimeErr
float m_peakTimeErr
Statistical error on peakTime.
Definition TOPLocalCalFitter.h:216
Belle2::TOP::TOPLocalCalFitter::m_channel
short m_channel
Channel number (0-511)
Definition TOPLocalCalFitter.h:201
Belle2::TOP::TOPLocalCalFitter::m_channelT0Err
float m_channelT0Err
Statistical error on channelT0.
Definition TOPLocalCalFitter.h:246
Belle2::TOP::TOPLocalCalFitter::m_detectCrosstalk
bool m_detectCrosstalk
Enables the crosstalk detection algorithm.
Definition TOPLocalCalFitter.h:157
Belle2::TOP::TOPLocalCalFitter::m_TTSData
std::string m_TTSData
File with the TTS parametrization.
Definition TOPLocalCalFitter.h:141
Belle2::TOP::TOPLocalCalFitter::m_binEdges
std::vector< float > m_binEdges
Amplitude bins.
Definition TOPLocalCalFitter.h:145
Belle2::TOP::TOPLocalCalFitter::m_sigmaBackground
float m_sigmaBackground
Sigma of the gaussian used to describe the background.
Definition TOPLocalCalFitter.h:229
Belle2::TOP::TOPLocalCalFitter::m_peakTime
float m_peakTime
Fitted time of the main (i.e.
Definition TOPLocalCalFitter.h:208
Belle2::TOP::TOPLocalCalFitter::m_fractionErr
float m_fractionErr
Statistical error on fraction.
Definition TOPLocalCalFitter.h:219
Belle2::TOP::TOPLocalCalFitter::m_isFitInAmplitudeBins
bool m_isFitInAmplitudeBins
Enables the fit in amplitude bins.
Definition TOPLocalCalFitter.h:144
Belle2::TOP::TOPLocalCalFitter::setHardwareIdentifiers
void setHardwareIdentifiers(short channel)
Set the hardware identifiers corresponding to a TOP channel.
Definition TOPLocalCalFitter.h:284
Belle2::TOP::TOPLocalCalFitter::m_alphaExtra
float m_alphaExtra
alpha parameter of the tail of the extra peak.
Definition TOPLocalCalFitter.h:225
Belle2::TOP::TOPLocalCalFitter::m_a0
float m_a0
Amplitude for channel 0 in pair.
Definition TOPLocalCalFitter.h:169
Belle2::TOP::TOPLocalCalFitter::m_rms
float m_rms
RMS of the histogram used for the fit.
Definition TOPLocalCalFitter.h:239
Belle2::TOP::TOPLocalCalFitter::m_mean2
float m_mean2
Position of the second gaussian of the TTS parametrization with respect to the first one.
Definition TOPLocalCalFitter.h:178
Belle2::TOP::TOPLocalCalFitter::m_fitConstraints
std::string m_fitConstraints
File with the Fit constraints.
Definition TOPLocalCalFitter.h:139
Belle2::TOP::TOPLocalCalFitter::m_hasChannelMaps
bool m_hasChannelMaps
Flag indicating if channel->(row,col) maps have been built.
Definition TOPLocalCalFitter.h:263
Belle2::TOP::TOPLocalCalFitter::m_timewalkTree
TTree * m_timewalkTree
Output of the fitter.
Definition TOPLocalCalFitter.h:153
Belle2::TOP::TOPLocalCalFitter::m_w0
float m_w0
Width for channel 0 in pair.
Definition TOPLocalCalFitter.h:171
Belle2::TOP::TOPLocalCalFitter::m_yieldLaser
float m_yieldLaser
Total number of laser hits from the fitting function integral.
Definition TOPLocalCalFitter.h:213
Belle2::TOP::TOPLocalCalFitter::areNeighbors
bool areNeighbors(short slot, short a, short b, int drMax=1, int dcMax=1) const noexcept
Return true if channels a and b are neighbors on the same slot in row/col space.
Definition TOPLocalCalFitter.h:298
Belle2::TOP::TOPLocalCalFitter::m_inputTTS
TFile * m_inputTTS
File containing m_treeTTS.
Definition TOPLocalCalFitter.h:146
Belle2::TOP::TOPLocalCalFitter::m_nExtra
float m_nExtra
parameter n of the tail of the extra peak
Definition TOPLocalCalFitter.h:226
Belle2::TOP::TOPLocalCalFitter::m_ht1
float m_ht1
Hit time for channel 1 in pair.
Definition TOPLocalCalFitter.h:168
Belle2::TOP::TOPLocalCalFitter::m_alphaExtraConstraints
float m_alphaExtraConstraints
alpha parameter of the tail of the extra peak.
Definition TOPLocalCalFitter.h:192
Belle2::TOP::TOPLocalCalFitter::m_channelT0
float m_channelT0
Raw, channelT0 calibration, defined as peakTime-peakTimeMC.
Definition TOPLocalCalFitter.h:241
Belle2::TOP::TOPLocalCalFitter::m_ch1
short m_ch1
Channel number (0-511)
Definition TOPLocalCalFitter.h:166
Belle2::TOP::TOPLocalCalFitter::m_sigmaErr
float m_sigmaErr
Statistical error on sigma.
Definition TOPLocalCalFitter.h:218
Belle2::TOP::TOPLocalCalFitter::m_binUpperEdge
float m_binUpperEdge
Upper edge of the amplitude bin in which this fit is performed.
Definition TOPLocalCalFitter.h:200
Belle2::TOP::TOPLocalCalFitter::m_f_q0
float m_f_q0
Fraction of charge on channel 0 in pair.
Definition TOPLocalCalFitter.h:175
Belle2::TOP::TOPLocalCalFitter::m_histFile
TFile * m_histFile
Output of the fitter.
Definition TOPLocalCalFitter.h:151
Belle2::TOP::TOPLocalCalFitter::m_treeConstraints
TTree * m_treeConstraints
Input to the fitter.
Definition TOPLocalCalFitter.h:149
Belle2::TOP::TOPLocalCalFitter::m_sigmaExtraConstraints
float m_sigmaExtraConstraints
Width of the gaussian used to describe the extra peak on the timing distribution tail.
Definition TOPLocalCalFitter.h:191
Belle2::TOP::TOPLocalCalFitter::m_w1
float m_w1
Width for channel 1 in pair.
Definition TOPLocalCalFitter.h:172
Belle2::TOP::TOPLocalCalFitter::m_firstPulserSigma
float m_firstPulserSigma
Time resolution from the fit of the first electronic pulse, from a Gaussian fit.
Definition TOPLocalCalFitter.h:249
Belle2::TOP::TOPLocalCalFitter::m_ht0
float m_ht0
Hit time for channel 0 in pair.
Definition TOPLocalCalFitter.h:167
Belle2::TOP::TOPLocalCalFitter::m_colOf
std::array< std::array< short, 512 >, 16 > m_colOf
Column index for (slot,channel), or -1 if out of bounds.
Definition TOPLocalCalFitter.h:262
Belle2::TOP::TOPLocalCalFitter::rowOf
short rowOf(short slot, short ch) const noexcept
Row index for (slot,channel), or -1 if out of bounds.
Definition TOPLocalCalFitter.h:266
Belle2::TOP::TOPLocalCalFitter::m_q1
float m_q1
Integrated charge for channel 1 in pair.
Definition TOPLocalCalFitter.h:174
Belle2::TOP::TOPLocalCalFitter::m_deltaT
float m_deltaT
Time difference between the main peak and the secondary peak.
Definition TOPLocalCalFitter.h:209
Belle2::TOP::TOPLocalCalFitter::m_fitTree
TTree * m_fitTree
Output of the fitter.
Definition TOPLocalCalFitter.h:152
Belle2::TOP::TOPLocalCalFitter::m_f2
float m_f2
Fraction of the second gaussian on the TTS parametrization.
Definition TOPLocalCalFitter.h:182
Belle2::TOP::TOPLocalCalFitter::determineFitStatus
void determineFitStatus()
determines if the constant obtained by the fit are good or not
Definition TOPLocalCalFitter.cc:580
Belle2::TOP::TOPLocalCalFitter::buildChannelMaps
void buildChannelMaps()
Build (row,col) lookup tables from the TTS tree; call after opening m_treeTTS.
Definition TOPLocalCalFitter.cc:34
Belle2::TOP::TOPLocalCalFitter::m_secondPulserTime
float m_secondPulserTime
Average time of the second electronic pulse respect to the reference pulse, from a gaussian fit.
Definition TOPLocalCalFitter.h:251
Belle2::TOP::TOPLocalCalFitter::m_boardstack
short m_boardstack
Boardstack number (0-3)
Definition TOPLocalCalFitter.h:207
Belle2::TOP::TOPLocalCalFitter::m_peakTimeConstraints
float m_peakTimeConstraints
Time of the main laser peak in the MC simulation (aka MC correction)
Definition TOPLocalCalFitter.h:187
Belle2::TOP::TOPLocalCalFitter::m_sigma1
float m_sigma1
Width of the first gaussian on the TTS parametrization.
Definition TOPLocalCalFitter.h:179
Belle2::TOP::TOPLocalCalFitter::m_fractionMC
float m_fractionMC
Fraction of events in the secondary peak form the MC simulation.
Definition TOPLocalCalFitter.h:232
Belle2::TOP::TOPLocalCalFitter::m_sigmaBackgroundConstraints
float m_sigmaBackgroundConstraints
Sigma of the gaussian used to describe the background.
Definition TOPLocalCalFitter.h:195
Belle2::TOP::TOPLocalCalFitter::m_treeTTS
TTree * m_treeTTS
Input to the fitter.
Definition TOPLocalCalFitter.h:148
Belle2::TOP::TOPLocalCalFitter::setupOutputTreeAndFile
void setupOutputTreeAndFile()
prepares the output tree
Definition TOPLocalCalFitter.cc:212
Belle2::TOP::TOPLocalCalFitter::~TOPLocalCalFitter
~TOPLocalCalFitter() override
Destructor.
Definition TOPLocalCalFitter.cc:54
Belle2::TOP::TOPLocalCalFitter::m_fitTree_noXtalk
TTree * m_fitTree_noXtalk
Output tree for non-crosstalk candidates.
Definition TOPLocalCalFitter.h:160
Belle2::TOP::TOPLocalCalFitter::m_timeBackgroundConstraints
float m_timeBackgroundConstraints
Position of the gaussian used to describe the background, w/ respect to peakTime.
Definition TOPLocalCalFitter.h:194
Belle2::TOP::TOPLocalCalFitter::m_ch0
short m_ch0
Channel number (0-511)
Definition TOPLocalCalFitter.h:165
Belle2::TOP::TOPLocalCalFitter::m_inputConstraints
TFile * m_inputConstraints
File containing m_treeConstraints.
Definition TOPLocalCalFitter.h:147
Belle2::TOP::TOPLocalCalFitter::m_deltaTErr
float m_deltaTErr
Statistical error on deltaT.
Definition TOPLocalCalFitter.h:217
Belle2::TOP::TOPLocalCalFitter::calculateChannelT0
void calculateChannelT0()
Calculates the commonT0 calibration after the fits have been done.
Definition TOPLocalCalFitter.cc:590
Belle2::TOP::TOPLocalCalFitter::m_rowOf
std::array< std::array< short, 512 >, 16 > m_rowOf
Row index for (slot,channel), or -1 if out of bounds.
Definition TOPLocalCalFitter.h:261
Belle2::TOP::TOPLocalCalFitter::m_firstPulserTime
float m_firstPulserTime
Average time of the first electronic pulse respect to the reference pulse, from a Gaussian fit.
Definition TOPLocalCalFitter.h:248
Belle2::TOP::TOPLocalCalFitter::m_histoIntegral
float m_histoIntegral
Integral of the fitted histogram.
Definition TOPLocalCalFitter.h:214
Belle2::TOP::TOPLocalCalFitter::m_yieldLaserExtra
float m_yieldLaserExtra
Integral of the extra peak.
Definition TOPLocalCalFitter.h:224
Belle2::TOP::TOPLocalCalFitter::m_sigma2
float m_sigma2
Width of the second gaussian on the TTS parametrization.
Definition TOPLocalCalFitter.h:180
Belle2::TOP::TOPLocalCalFitter::m_a1
float m_a1
Amplitude for channel 1 in pair.
Definition TOPLocalCalFitter.h:170
Belle2::TOP::TOPLocalCalFitter::fitPulser
void fitPulser(TH1 *, TH1 *)
Fits the two pulsers.
Definition TOPLocalCalFitter.cc:544
Belle2::TOP::TOPLocalCalFitter::calibrate
EResult calibrate() override
Runs the algorithm on events.
Definition TOPLocalCalFitter.cc:646
Belle2::TOP::TOPLocalCalFitter::m_timeExtra
float m_timeExtra
Position of the extra peak seen in the timing tail, w/ respect to peakTime.
Definition TOPLocalCalFitter.h:222
Belle2::TOP::TOPLocalCalFitter::m_fractionConstraints
float m_fractionConstraints
Fraction of the main peak.
Definition TOPLocalCalFitter.h:189
Belle2::TOP::TOPLocalCalFitter::m_deltaTConstraints
float m_deltaTConstraints
Distance between the main and the secondary laser peak.
Definition TOPLocalCalFitter.h:188
Belle2::CalibrationAlgorithm::getObjectPtr
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...
Definition CalibrationAlgorithm.h:414
Belle2::Hit
Structure to hold some of the calpulse data.
Definition TOPTimeBaseCalibratorModule.h:30