Belle II Software  release-08-01-10
laserResolutionTools.py
1 #!/usr/bin/env python3
2 # -*- coding: utf-8 -*-
3 
4 
11 
12 # ------------------------------------------------------------------------
13 # Tools to fit the laser timing distribution and plot the results.
14 # This code is intended for studying the TOP time resolution
15 # ------------------------------------------------------------------------
16 
17 
18 import basf2 as b2
19 from ROOT import TF1, TFile, TH1F, TSpectrum, TTree
20 from array import array
21 
22 
23 def fitLaserResolution(
24  mcCorrectionsFile='/group/belle2/group/detector/TOP/calibration/MCreferences/t0MC.root',
25  dataFile='',
26  outputFile='laserResolutionResults.root',
27  pdfType='cb',
28  maxPeaks=2,
29  saveFits=False,
30  includeBackscatter=True,
31  useSinglePDF=True,
32  useMCPeaks=True):
33  """
34  Function that implements the fit of the laser time distribution correcting for the propagation time in the prism.
35  The inputs to the fitter are hit time distributions per channel, in the form of TH2F.
36 
37  Two files containing the containing the hit timing VS channels distribution must be provided:
38 
39  * The first file is used to estimate the light propagation correction (aka MC correction),
40  the number of peaks and the separation bewteen these peaks.
41  Since the correction is the same for all the slots, the corresponding TH2F can have only 512 bins on the x-axis.
42  A default file is available, where the hit time is taken from
43  the TOPSimHits, and the laser jitter has been set to zero in the simulation.
44  * The second file contains the distribution obtained from the data do be used to estimate the resolution.
45  In this case the histrogram must have
46  (512 x 16) bins on the x-axis, to store the time distribution for each channel of the detector.
47 
48  Parameters:
49  mcCorrectionsFile (str): File with the TH2F containing the TOPSImHits time distribution.
50  Used to calculate the MC corrections
51  dataFile (str): File with the TH2F containing the TOPDigit time distribution on which the resolutions
52  must be evaluated
53  outputFile (str): rootFile containing the outupt trees and plots
54  pdfType (str): PDF used in the data fit: cb for a CrystalBall, gaus for a Gaussian or
55  gausExpo for an exponentially modified Gaussian
56  maxPeaks (int): maximum number of peaks that can be used in the data fit. Note that the
57  actual value is decided channel-by-channel
58  according to the MC simulation, this parameter is only a maximum cap.
59  In the current implementation only 1 or 2 peaks are supported, therefore maxPeaks can only go up to 2
60  saveFits (bool): set to True to save the histrogram and the fit for eavey single channel
61  includeBackscatter (bool): set to True to add an extra peak in the fit to account for the
62  MCP backscattering
63  useMCPeaks (bool): ste to True to fix the light path time difference from the MC. Set to False
64  to let this parameter floating in the fit
65  Returns:
66  void
67  """
68 
69  if maxPeaks > 2 or maxPeaks < 1:
70  b2.B2FATAL('Usupported value for the maximum number for peaks (maxPeaks = ' +
71  str(maxPeaks) + '). Please set maxPeak o either 1 or 2.')
72  if pdfType != 'cb' and pdfType != 'gaus' and pdfType != 'gausExpo':
73  b2.B2FATAL(
74  'Unknown PDF type ' +
75  pdfType +
76  ". The possible options are cb for the Crystal Ball,"
77  "gaus for the Gaussian and gausExpo for exponentially modified Gaussian")
78 
79  # First loops over the MC file to get the time propagation corrections
80  # Only the two highest peaks are stored
81 
82  # List with the position fo the MC peaks (2 per channel). This is where the
83  # MC peak positions are stored to be used later to correct the data
84  mcPeaks = [[-1. for second in range(2)] for first in range(512)]
85 
86  # Opens the file containg the MC distribution and loads it
87  tFileMC = TFile(mcCorrectionsFile)
88  histoMC = tFileMC.Get('LaserTimingVSChannelOneSlot')
89 
90  # Loop over the channels of a slot (i.e the bins on the x axis of histoMC)
91  for kCh in range(512):
92  # Gets a one-channel wide slice and adjusts the range
93  timeProjectionMC = histoMC.ProjectionY('projectionMC', kCh + 1, kCh + 1)
94  timeProjectionMC.GetXaxis().SetRangeUser(0., 1.)
95  # Initializes TSpectrum
96  spectrum = TSpectrum()
97  # By default dentifies a peak only if it is 1-sigma apart from another peak
98  # and has an aplitud of at least 10% of the largest peak.
99  numPeaks = spectrum.Search(timeProjectionMC, 1., 'nobackground', 0.1)
100  # This is a bit ugl way to store the peak locations and amplitudes
101  # in a list. At this stage the maxima are no granted to be sorted
102  peaksX_ = spectrum.GetPositionX()
103  peaksY_ = spectrum.GetPositionY()
104  peaksX = [peaksX_[i] for i in range(numPeaks)]
105  peaksY = [peaksY_[j] for j in range(numPeaks)]
106 
107  # Sorts the peaks starting form the largest one and stores the first
108  # two in the mcPeaks matrix
109  for iMax in range(numPeaks):
110  if iMax < 2:
111  maxPosition = int(peaksY.index(max(peaksY)))
112  peaksY[maxPosition] = -1. # get rid of the maximum
113  mcPeaks[kCh][iMax] = peaksX[maxPosition]
114 
115  tFileMC.Close()
116 
117  # End of the MC corrections part
118 
119  # Opens the file with the data and loads the timing histogram
120  tFileData = TFile(dataFile)
121  histoData = tFileData.Get('LaserTimingVSChannel')
122 
123  # Prepares the output file and tree
124  # tFileOut = TFile(outputFile, "RECREATE")
125  outTree = TTree("tree", "tree")
126 
127  # Slot ID (1-based)
128  slotID = array('i', [0])
129  # Hardware channle number (0-based)
130  hwChannel = array('i', [0])
131  # Number of peaks used in the fit
132  nPeaks = array('i', [0])
133  # Degrees of freedom of the fit
134  ndf = array('f', [0])
135  # Chi2 of the fit
136  chi2 = array('f', [0.])
137  # Sigma of the Gaussian core
138  sigma = array('f', [0.])
139  # Statistical error on sigma
140  errSigma = array('f', [0.])
141  # Time of the primary peak (maximum of the fit fuction)
142  peakTime = array('f', [0.])
143  # Statistical error on manPeak
144  errPeakTime = array('f', [0.])
145 
146  # Time of the spurious peak
147  extraTime = array('f', [0.])
148  # Statistical error on extraTime
149  errExtraTime = array('f', [0.])
150 
151  # scale factor for Sigma of the Gaussian core
152  sigmaExtra = array('f', [0.])
153  # Statistical error on sigma scale factor
154  errSigmaExtra = array('f', [0.])
155 
156  # Alpha parameter of the CB functions
157  alpha = array('f', [0.])
158  # Statistical error on alpha
159  errAlpha = array('f', [0.])
160  # n parameter of the CB functions
161  n = array('f', [0.])
162  # statistical error on n
163  errN = array('f', [0.])
164  # lambda parameter of the gauss exponiential
165  lam = array('f', [0.])
166  # Statistical error on lam
167  errLam = array('f', [0.])
168  # Normalization of the primary peak
169  norm1 = array('f', [0.])
170  # Statistical error on norm1
171  errNorm1 = array('f', [0.])
172  # Normalization of the secondary peak
173  norm2 = array('f', [0.])
174  # Statistical error on norm2
175  errNorm2 = array('f', [0.])
176 
177  # Normalization of the spurious peak
178  normExtra = array('f', [0.])
179  # Statistical error on normExtra
180  errNormExtra = array('f', [0.])
181 
182  # Separation between primary and secondary peak from the MC
183  deltaMC = array('f', [0.])
184  # Laser time corrected for the propagation time in the prism
185  correctedTime = array('f', [0.])
186  # Time of the late peak in the MC
187  latePeak = array('f', [0.])
188  # Time of the earlier peak in the MC
189  earlyPeak = array('f', [0.])
190 
191  # Sets the tree branches
192  outTree.Branch('slotID', slotID, 'slotID/I')
193  outTree.Branch('hwChannel', hwChannel, 'hwChannel/I')
194  outTree.Branch('ndf', ndf, 'ndf/F')
195  outTree.Branch('chi2', chi2, 'chi2/F')
196  outTree.Branch('nPeaks', nPeaks, 'nPeaks/I')
197  outTree.Branch('sigma', sigma, 'sigma/F')
198  outTree.Branch('errSigma', errSigma, 'errSigma/F')
199  outTree.Branch('peakTime', peakTime, 'peakTime/F')
200  outTree.Branch('errPeakTime', errPeakTime, 'errPeakTime/F')
201  outTree.Branch('extraTime', extraTime, 'extraTime/F')
202  outTree.Branch('errExtraTime', errExtraTime, 'errExtraTime/F')
203  outTree.Branch('alpha', alpha, 'alpha/F')
204  outTree.Branch('errAlpha', errAlpha, 'errAlpha/F')
205  outTree.Branch('n', n, 'n/F')
206  outTree.Branch('errN', errN, 'errN/F')
207  outTree.Branch('lam', lam, 'lam/F')
208  outTree.Branch('errLam', errLam, 'errLam/F')
209  outTree.Branch('norm1', norm1, 'norm1/F')
210  outTree.Branch('errNorm1', errNorm1, 'errNorm1/F')
211  outTree.Branch('norm2', norm2, 'norm2/F')
212  outTree.Branch('errNorm2', errNorm2, 'errNorm2/F')
213  outTree.Branch('normExtra', normExtra, 'normExtra/F')
214  outTree.Branch('errNormExtra', errNormExtra, 'errNormExtra/F')
215  outTree.Branch('sigmaExtra', sigmaExtra, 'sigmaExtra/F')
216  outTree.Branch('errSigmaExtra', errSigmaExtra, 'errSigmaExtra/F')
217  outTree.Branch('deltaMC', deltaMC, 'deltaMC/F')
218  outTree.Branch('correctedTime', correctedTime, 'correctedTime/F')
219  outTree.Branch('latePeak', latePeak, 'latePeak/F')
220  outTree.Branch('earlyPeak', earlyPeak, 'earlyPeak/F')
221 
222  # Main loop_: takes the histogram, slices it and fits each slide according to the number of peaks found in the MC and the
223  # pdf chosen by te user
224  for iSlot in range(0, 16):
225  print('Fitting slot ' + str(iSlot))
226  for k in range(0, 512):
227  timeProjection = histoData.ProjectionY('projection_' + str(iSlot + 1) + '_' +
228  str(k), k + iSlot * 512 + 1, k + iSlot * 512 + 1)
229  maximum = timeProjection.GetXaxis().GetBinCenter(timeProjection.GetMaximumBin())
230  if timeProjection.Integral() < 10:
231  continue
232  timeProjection.GetXaxis().SetRangeUser(maximum - 1., maximum + 2.)
233 
234  slotID[0] = iSlot + 1
235  hwChannel[0] = k
236 
237  # Saves the MC information in the tree variables
238  if(mcPeaks[k][1] > 0):
239  # By convention the first component of the PDF represents the later peak,
240  # which in the MC appears to be also the main one. However this may change with
241  # different MC simulations, therefore we have to check which one is
242  # the first peak and which one is the second peak
243  if mcPeaks[k][1] > mcPeaks[k][0]:
244  earlyPeak[0] = mcPeaks[k][0]
245  latePeak[0] = mcPeaks[k][1]
246  else:
247  earlyPeak[0] = mcPeaks[k][1]
248  latePeak[0] = mcPeaks[k][0]
249  deltaMC[0] = earlyPeak[0] - latePeak[0]
250  nPeaks[0] = 2
251  else:
252  nPeaks[0] = 1
253  deltaMC[0] = 0.
254  latePeak[0] = mcPeaks[k][0]
255  earlyPeak[0] = 0.
256 
257  # fit with the CB
258  if pdfType == 'cb':
259  cb = TF1(
260  "cb",
261  "[0]*ROOT:: Math: : crystalball_function(x, [2], [3], [4], [5]) + "
262  "[1]*[0]*ROOT::Math::crystalball_function(x, [2], [3], [4], [5]+[6]) + "
263  "[7]*[0]*ROOT:: Math::crystalball_function(x, [2], [3], [9]*[4], [5]+[8]) ",
264  maximum -
265  1.,
266  maximum +
267  2.)
268  # main peak
269  cb.SetParameter(0, 10.) # norm 1
270  cb.SetParameter(2, -1.) # alpha
271  cb.SetParLimits(2, -30., -0.01) # alpha
272  cb.SetParameter(3, 3.) # n
273  cb.SetParLimits(3, 0.01, 6.) # n must be on the right side
274  cb.SetParameter(4, 0.05) # sigma
275  cb.SetParLimits(4, 0.03, 0.200) # range for sigma
276  cb.SetParameter(5, maximum) # maximum of the latest peak
277  cb.SetParLimits(5, maximum - 0.1, maximum + 0.1) # maximum of the latest peak
278  # backscatter peak
279  cb.SetParameter(7, 0.01) # yield of the spurious peak
280  cb.SetParLimits(7, 0., .1) # yield of the spurious peak
281  cb.SetParameter(8, 1.) # delay of the spurious peak (~ 1 ns)
282  cb.SetParLimits(8, 0.2, 2.) # delay of the spurious peak (~ 1 ns)
283  cb.SetParameter(9, 1.5) # scale factro for sigma of the spurious peak (~2)
284  cb.SetParLimits(9, 1.1, 3.) # scale factro for sigma of the spurious peak (~2)
285  if not includeBackscatter:
286  cb.FixParameter(7, 0.) # yield of the spurious peak
287  cb.FixParameter(8, 1.) # delay of the spurious peak (~ 1 ns)
288  cb.FixParameter(9, 1.5) # scale factro for sigma of the spurious peak (~2)
289 
290  if(mcPeaks[k][1] > 0 and not useSinglePDF):
291  cb.SetParameter(1, 0.3) # norm 2
292  cb.SetParLimits(1, 0., 100.) # norm 1
293  if useMCPeaks:
294  cb.FixParameter(6, deltaMC[0])
295  else:
296  cb.SetParameter(6, deltaMC[0])
297  cb.SetParLimits(6, deltaMC[0] - 0.2, -0.001)
298  else:
299  cb.FixParameter(1, 0.)
300  cb.FixParameter(6, 0.)
301 
302  cb.SetNpx(1000)
303  timeProjection.Fit('cb', 'RQS')
304  timeProjection.Fit('cb', 'RQS')
305  result = timeProjection.Fit('cb', 'RQS')
306 
307  # Dumps the fit results into the tree variables
308  chi2[0] = result.Chi2()
309  ndf[0] = result.Ndf()
310  alpha[0] = cb.GetParameter(2)
311  errAlpha[0] = cb.GetParError(2)
312  n[0] = cb.GetParameter(3)
313  errN[0] = cb.GetParError(3)
314  lam[0] = 0.
315  errLam[0] = 0.
316  peakTime[0] = cb.GetParameter(5)
317  errPeakTime[0] = cb.GetParError(5)
318  sigma[0] = cb.GetParameter(4)
319  errSigma[0] = cb.GetParError(4)
320  norm1[0] = cb.GetParameter(0)
321  errNorm1[0] = cb.GetParError(0)
322  norm2[0] = cb.GetParameter(1)
323  errNorm2[0] = cb.GetParError(1)
324  normExtra[0] = cb.GetParameter(7)
325  errNormExtra[0] = cb.GetParError(7)
326  extraTime[0] = cb.GetParameter(8)
327  errExtraTime[0] = cb.GetParError(8)
328  sigmaExtra[0] = cb.GetParameter(9)
329  errSigmaExtra[0] = cb.GetParError(9)
330 
331  # models each peak as an exponentially-modified gaussian
332  if pdfType == 'gausExpo':
333 
334  gausExpo = TF1(
335  "gausExpo",
336  "[0]*(([4]/2.)*TMath::Exp(([4]/2.)*(2*[2] + [4]*[3]*[3]- 2.*x)))*"
337  "ROOT::Math::erfc(([2] + [4]*[3]*[3] - x)/(TMath::Sqrt(2.)*[3])) + "
338  "[1]*[0]*(([4]/2.)*TMath::Exp(([4]/2.)*(2*([2]+[5])+[4]*[3]*[3]-2.*x)))*"
339  "ROOT::Math::erfc((([2]+[5]) + [4]*[3]*[3] - x)/(TMath::Sqrt(2.)*[3])) +"
340  "[6]*[0]*(([4]/2.)*TMath::Exp(([4]/2.)*(2*([2]+[7])+[4]*[8]*[3]*[8]*[3]-2.*x)))*"
341  "ROOT::Math::erfc((([2]+[7]) + [4]*[8]*[3]*[8]*[3] - x)/(TMath::Sqrt(2.)*[8]*[3]))",
342  maximum -
343  1.,
344  maximum +
345  2.)
346 
347  # main peak
348  gausExpo.SetParameter(0, 1000.) # norm 1
349  gausExpo.SetParameter(2, maximum) # maximum of the latest peak
350  gausExpo.SetParLimits(2, maximum - 0.1, maximum + 0.1) # maximum of the latest peak
351  gausExpo.SetParameter(3, 0.05) # sigma
352  gausExpo.SetParLimits(3, 0.03, 0.200) # range for sigma
353  gausExpo.SetParameter(4, 1.2) # lambda
354  # backscatter peak
355  gausExpo.SetParameter(6, 0.01) # norm extra
356  gausExpo.SetParLimits(6, 0., 0.1) # norm extra
357  gausExpo.SetParameter(7, 1.) # delay of the spurious peak (~ 1 ns)
358  gausExpo.SetParLimits(7, 0.2, 2.) # delay of the spurious peak (~ 1 ns)
359  gausExpo.SetParameter(8, 1.5) # scale factro for sigma of the spurious peak (~2)
360  gausExpo.SetParLimits(8, 1.1, 3.) # scale factro for sigma of the spurious peak (~2)
361  if not includeBackscatter:
362  gausExpo.FixParameter(6, 0.) # yield of the spurious peak
363  gausExpo.FixParameter(7, 1.) # delay of the spurious peak (~ 1 ns)
364  gausExpo.FixParameter(8, 1.5) # scale factor for sigma of the spurious peak (~2)
365 
366  if(mcPeaks[k][1] > 0 and not useSinglePDF):
367  gausExpo.SetParameter(1, 0.3) # norm 2
368  gausExpo.SetParLimits(1, 0., 100.) # norm 1
369  if useMCPeaks:
370  gausExpo.FixParameter(5, deltaMC[0])
371  else:
372  gausExpo.SetParameter(5, deltaMC[0])
373  gausExpo.SetParLimits(5, deltaMC[0] - 0.2, -0.001)
374  else:
375  gausExpo.FixParameter(1, 0.)
376  gausExpo.FixParameter(5, 0.)
377 
378  # Perform the fit
379  gausExpo.SetNpx(1000)
380  timeProjection.Fit('gausExpo', 'RQS')
381  timeProjection.Fit('gausExpo', 'RQS')
382  result = timeProjection.Fit('gausExpo', 'RQS')
383 
384  # save the results
385  chi2[0] = result.Chi2()
386  ndf[0] = result.Ndf()
387  alpha[0] = 0.
388  errAlpha[0] = 0.
389  n[0] = 0.
390  errN[0] = 0.
391  lam[0] = gausExpo.GetParameter(4)
392  errLam[0] = gausExpo.GetParError(4)
393  peakTime[0] = gausExpo.GetParameter(2)
394  errPeakTime[0] = gausExpo.GetParError(2)
395  sigma[0] = gausExpo.GetParameter(3)
396  errSigma[0] = gausExpo.GetParError(3)
397  norm1[0] = gausExpo.GetParameter(0)
398  errNorm1[0] = gausExpo.GetParError(0)
399  norm2[0] = gausExpo.GetParameter(1)
400  errNorm2[0] = gausExpo.GetParError(1)
401  normExtra[0] = gausExpo.GetParameter(6)
402  errNormExtra[0] = gausExpo.GetParError(6)
403  extraTime[0] = gausExpo.GetParameter(7)
404  errExtraTime[0] = gausExpo.GetParError(7)
405  sigmaExtra[0] = gausExpo.GetParameter(8)
406  errSigmaExtra[0] = gausExpo.GetParError(8)
407 
408  # fit using a simple gaussian
409  if pdfType == 'gaus':
410  gaussian = TF1(
411  "gaussian",
412  "[0]*TMath::Gaus([2], [3], kTRUE)) + "
413  "[0]*[1]*TMath::Gaus([2]+[4], [3], kTRUE)) +"
414  "[0]*[5]*TMath::Gaus([2]+[6], [3]*[7], kTRUE))",
415  maximum -
416  1.,
417  maximum +
418  2.)
419  # main peak
420  gaussian.SetParameter(0, 10.) # norm 1
421  gaussian.SetParameter(2, maximum) # maximum of the latest peak
422  gaussian.SetParLimits(2, maximum - 0.1, maximum + 0.1) # maximum of the latest peak
423  gaussian.SetParameter(3, 0.05) # sigma
424  gaussian.SetParLimits(3, 0.03, 0.200) # range for sigma
425  # backscatter peak
426  gaussian.SetParameter(5, 0.01) # yield of the spurious peak
427  gaussian.SetParLimits(5, 0., .1) # yield of the spurious peak
428  gaussian.SetParameter(6, 1.) # delay of the spurious peak (~ 1 ns)
429  gaussian.SetParLimits(6, 0.2, 2.) # delay of the spurious peak (~ 1 ns)
430  gaussian.SetParameter(7, 1.5) # scale factro for sigma of the spurious peak (~2)
431  gaussian.SetParLimits(7, 1.1, 3.) # scale factro for sigma of the spurious peak (~2)
432  if not includeBackscatter:
433  gaussian.FixParameter(5, 0.) # yield of the spurious peak
434  gaussian.FixParameter(6, 1.) # delay of the spurious peak (~ 1 ns)
435  gaussian.FixParameter(7, 1.5) # scale factro for sigma of the spurious peak (~2)
436 
437  if(mcPeaks[k][1] > 0 and not useSinglePDF):
438  gaussian.SetParameter(1, 0.3) # norm 2
439  gaussian.SetParLimits(1, 0., 100.) # norm 1
440  if useMCPeaks:
441  gaussian.FixParameter(4, deltaMC[0])
442  else:
443  gaussian.SetParameter(4, deltaMC[0])
444  gaussian.SetParLimits(4, deltaMC[0] - 0.2, -0.001)
445  else:
446  gaussian.FixParameter(1, 0.)
447  gaussian.FixParameter(4, 0.)
448 
449  gaussian.SetNpx(1000)
450  timeProjection.Fit('gaussian', 'RQS')
451  timeProjection.Fit('gaussian', 'RQS')
452  result = timeProjection.Fit('gaussian', 'RQS')
453 
454  # Dumps the fit results into the tree variables
455  chi2[0] = result.Chi2()
456  ndf[0] = result.Ndf()
457  alpha[0] = 0.
458  errAlpha[0] = 0.
459  n[0] = 0.
460  errN[0] = 0.
461  lam[0] = 0.
462  errLam[0] = 0.
463  peakTime[0] = gaussian.GetParameter(2)
464  errPeakTime[0] = gaussian.GetParError(2)
465  sigma[0] = gaussian.GetParameter(3)
466  errSigma[0] = gaussian.GetParError(3)
467  norm1[0] = gaussian.GetParameter(0)
468  errNorm1[0] = gaussian.GetParError(0)
469  norm2[0] = gaussian.GetParameter(1)
470  errNorm2[0] = gaussian.GetParError(1)
471  normExtra[0] = gaussian.GetParameter(5)
472  errNormExtra[0] = gaussian.GetParError(5)
473  extraTime[0] = gaussian.GetParameter(6)
474  errExtraTime[0] = gaussian.GetParError(6)
475  sigmaExtra[0] = gaussian.GetParameter(7)
476  errSigmaExtra[0] = gaussian.GetParError(7)
477 
478  # calculate the MC-corrected time sobtracting the timing of the main peak in the MC
479  correctedTime[0] = peakTime[0] - mcPeaks[k][0]
480 
481  outTree.Fill()
482 
483  if saveFits:
484  timeProjection.Write()
485 
486  histoData.Write()
487  outTree.Write()
488 
489 
490 def plotLaserResolution(fitResultsFile='laserResolutionResults.root'):
491  """
492  Function that creates some standard plots from the tree containg the flaset fit results
493  """
494 
495  # Opens the file with the data and loads the timing histogram
496  tFileData = TFile(fitResultsFile)
497  treeData = tFileData.Get('tree')
498  histoData = tFileData.Get('LaserTimingVSChannel')
499 
500  # Findes the average value of correctedTime to optimize the hisotgram axis ranges
501  h_test = TH1F("h_test", "a test histogram to determine the histogram boundaries", 1000, 0., 100.)
502  treeData.Draw("correctedTime>>h_test")
503  meanCorrTime = h_test.GetMean()
504 
505  # laser peak time in each slot, after the MC correction (10 ps binnig)
506  h_localT0Resolution = [TH1F("h_localT0Resolution_slot" + str(iSlot + 1),
507  "Laser peak position after the MC correction in slot " + str(iSlot + 1),
508  400,
509  meanCorrTime - 2.,
510  meanCorrTime + 2.) for iSlot in range(16)]
511  # laser peak resolution in the whole detector (10 ps binnig)
512  h_singleChannelResolution = TH1F("h_singleChannelResolution", "Resolution of the laser peak.", 40, 0., 0.4)
513 
514  # laser peak resolution as function of the channel (10 ps binnig)
515  h_singleChannelResolutionVSChannel = TH1F(
516  "singleChannelResolutionVSChannel",
517  "Resolution of the laser peak VS channel number",
518  512 * 16,
519  0,
520  512 * 16)
521 
522  # laser peak position as function of the channel (10 ps binnig)
523  h_localT0ResolutionVSChannel = TH1F(
524  "localT0ResolutionVSChannel",
525  "Resolution of the laser peak VS channel number",
526  512 * 16,
527  0,
528  512 * 16)
529 
530  # Distribution of the late peak
531  h_latePeakVSChannel = TH1F("latePeakVSChannel", "Late Peak position", 512 * 16, 0, 512 * 16)
532  # Distribution of the early peak
533  h_earlyPeakVSChannel = TH1F("earlyPeakVSChannel", "Early Peak position", 512 * 16, 0, 512 * 16)
534 
535  for entry in treeData:
536  globalChannel = entry.hwChannel + entry.slotID * 512
537 
538  h_localT0Resolution[entry.slotID - 1].Fill(entry.correctedTime)
539  h_singleChannelResolution.Fill(entry.sigma)
540 
541  h_singleChannelResolutionVSChannel.SetBinContent(globalChannel + 1, entry.sigma)
542  h_singleChannelResolutionVSChannel.SetBinError(globalChannel + 1, entry.errSigma)
543 
544  h_localT0ResolutionVSChannel.SetBinContent(globalChannel + 1, entry.correctedTime)
545  h_localT0ResolutionVSChannel.SetBinError(globalChannel + 1, entry.errPeakTime)
546 
547  h_latePeakVSChannel.SetBinContent(globalChannel + 1, entry.latePeak)
548  if entry.earlyPeak > 0:
549  h_earlyPeakVSChannel.SetBinContent(globalChannel + 1, entry.earlyPeak)
550 
551  # loops over the MC distributions for shift the according to the time offset of each slot
552  for iSlot in range(16):
553  laserTime = histoData.GetBinContent(iSlot * 512 + 1) # channle 0 of the slot
554  mcTime = h_latePeakVSChannel.GetBinContent(iSlot * 512 + 1) # channle 0 of the slot
555  timeShift = laserTime - mcTime
556  for iBin in range(512):
557  lateTime = h_latePeakVSChannel.GetBinContent(iBin + iSlot * 512 + 1) + timeShift
558  earlyTime = h_earlyPeakVSChannel.GetBinContent(iBin + iSlot * 512 + 1) + timeShift
559  h_latePeakVSChannel.SetBinContent(iBin + iSlot * 512 + 1, lateTime)
560  h_earlyPeakVSChannel.SetBinContent(iBin + iSlot * 512 + 1, earlyTime)
561 
562  # Prepares the output file and tree
563  tFileOut = TFile('plots.root', "RECREATE")
564 
565  histoData.Write()
566  h_latePeakVSChannel.Write()
567  h_earlyPeakVSChannel.Write()
568  h_singleChannelResolution.Write()
569  h_singleChannelResolutionVSChannel.Write()
570  h_localT0ResolutionVSChannel.Write()
571 
572  for i in range(16):
573  h_localT0Resolution[i].Write()
574 
575  tFileOut.Close()