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