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