Belle II Software  release-05-01-25
cdst_calibrateCommonT0.py
1 #!/usr/bin/env python3
2 # -*- coding: utf-8 -*-
3 
4 # ---------------------------------------------------------------------------------------
5 # Calibrate common T0 with Bhabha (or dimuon) events using new constants from DB
6 # (M. Staric, 2019-07-10)
7 #
8 # usage: basf2 cdst_calibrateCommonT0.py experiment run
9 # job: bsub -q s "basf2 cdst_calibrateCommonT0.py experiment run"
10 # ---------------------------------------------------------------------------------------
11 
12 from basf2 import *
13 from ROOT import Belle2
14 from ROOT import TH1F, TH2F, TFile, TTree, TF1
15 from array import array
16 import math
17 import sys
18 import glob
19 import os
20 
21 # ----- those need to be adjusted before running --------------------------------------
22 #
23 sampleType = 'bhabha' # sample type: 'bhabha' or 'dimuon'
24 data_dir = '/group/belle2/dataprod/Data/release-03-02-02/DB00000635/proc00000009_nofilter'
25 skim_dir = 'skim/hlt_bhabha/cdst/sub00'
26 globalTag = 'data_reprocessing_prompt' # base global tag (fall-back)
27 stagingTags = ['staging_data_reprocessing'] # list of global tags with new calibration
28 localDB = [] # list of local databases with new calibration
29 method = 'BF' # BF: bunch offset fit, LL: likelihood method
30 output_dir = 'commonT0' # main output folder
31 #
32 # -------------------------------------------------------------------------------------
33 
34 # Argument parsing
35 argvs = sys.argv
36 if len(argvs) < 3:
37  print("usage: basf2", argvs[0], "experiment run")
38  sys.exit()
39 experiment = int(argvs[1])
40 run = int(argvs[2])
41 
42 expNo = 'e' + '{:0=4d}'.format(experiment)
43 runNo = 'r' + '{:0=5d}'.format(run)
44 
45 # Make list of files
46 files = []
47 for typ in ['4S', 'Continuum', 'Scan']:
48  folder = data_dir + '/' + expNo + '/' + typ + '/' + runNo + '/' + skim_dir
49  files += glob.glob(folder + '/cdst.*.root')
50 if len(files) == 0:
51  B2ERROR('No cdst files found')
52  sys.exit()
53 
54 # Output folder
55 output_folder = output_dir + '/' + expNo + '/' + sampleType + '/' + method
56 if not os.path.isdir(output_folder):
57  os.makedirs(output_folder)
58  print('New folder created: ' + output_folder)
59 
60 # Output file name
61 fileName = output_folder + '/commonT0-' + expNo + '-' + runNo + '.root'
62 print('Output file:', fileName)
63 
64 
65 class Mask_BS13d(Module):
66  ''' exclude (mask-out) BS 13d '''
67 
68  def event(self):
69  ''' event processing '''
70 
71  for digit in Belle2.PyStoreArray('TOPDigits'):
72  if digit.getModuleID() == 13 and digit.getBoardstackNumber() == 3:
73  digit.setHitQuality(Belle2.TOPDigit.c_Junk)
74 
75 
77  """
78  ** Description **
79  Module to calculate the commont time offset of the TOP detector.
80  It selects two-track events, saves the bunch finder time offset in an histrogram and
81  fits it to extract the global time offset with respect to the RF clock.
82  It should be run in an individual job for each run, but in future it may be further
83  authomatize to run on multiple runs in the same job.
84 
85  **Output**
86  Root file containing:
87  * The histogram of the bunch finder time offset
88  * The fit function
89  * The fit parameters ina tree format (usefuly when merging several files)
90 
91  **Contributors**
92  * Marko Staric
93  * Umberto Tamponi
94  """
95 
96  def initialize(self):
97  """
98  Creates the histogram used for the commoT0 calculation, and
99  takes the necessary objects from the DataStore
100  """
101 
102 
103  self.file = TFile.Open(fileName, 'recreate')
104 
105  geo = Belle2.PyDBObj('TOPGeometry')
106  if not geo.isValid():
107  B2FATAL('TOP geometry not available in database')
108 
109 
110  self.bunchTimeSep = geo.getNominalTDC().getSyncTimeBase() / 24
111 
112 
113  self.h1 = TH1F("offset", "current offset; offset [ns]", 600, -9.0, 9.0)
114 
115 
116  self.h2 = TH2F("offset_vs_event", "current offset versus event number",
117  100, 0.0, 1000000.0, 200, -3.0, 3.0)
118  self.h2.SetXTitle("event number")
119  self.h2.SetYTitle("offset [ns]")
120 
121  xmi = -self.bunchTimeSep / 2
122  xma = self.bunchTimeSep / 2
123 
124  self.h1a = TH1F("offset_a", "current offset; offset [ns]", 200, xmi, xma)
125 
126 
127  self.h2a = TH2F("offset_vs_event_a", "current offset versus event number",
128  100, 0.0, 1000000.0, 200, xmi, xma)
129  self.h2a.SetXTitle("event number")
130  self.h2a.SetYTitle("offset [ns]")
131 
132  xmi = 0.0
133  xma = self.bunchTimeSep
134 
135  self.h1b = TH1F("offset_b", "current offset; offset [ns]", 200, xmi, xma)
136 
137 
138  self.h2b = TH2F("offset_vs_event_b", "current offset versus event number",
139  100, 0.0, 1000000.0, 200, xmi, xma)
140  self.h2b.SetXTitle("event number")
141  self.h2b.SetYTitle("offset [ns]")
142 
143  evtMetaData = Belle2.PyStoreObj('EventMetaData')
144 
145 
146  self.expNo = evtMetaData.getExperiment()
147 
148 
149  self.run = '{:0=5d}'.format(evtMetaData.getRun())
150 
151 
152  self.t0 = 0
153 
154  self.t0Err = 0
155  commonT0 = Belle2.PyDBObj('TOPCalCommonT0')
156  if commonT0.isValid():
157  if commonT0.isCalibrated():
158  self.t0 = commonT0.getT0()
159  self.t0Err = commonT0.getT0Error()
160  B2INFO('Common T0 used in data processing: ' + str(self.t0))
161  else:
162  B2INFO('No common T0 calibration done yet')
163  else:
164  B2ERROR('Common T0 not available in database')
165 
166  def event(self):
167  """
168  Fills the histogram of the time offset
169  """
170  recBunch = Belle2.PyStoreObj('TOPRecBunch')
171  evtMetaData = Belle2.PyStoreObj('EventMetaData')
172  if not recBunch.isValid():
173  return
174  if recBunch.isReconstructed() and recBunch.getNumTracks() == 2:
175  offset = recBunch.getCurrentOffset()
176  evtNum = evtMetaData.getEvent()
177  self.h1.Fill(offset)
178  self.h2.Fill(evtNum, offset)
179  # wrap-around into [-1/2, 1/2] of bunch cycle
180  a = offset - round(offset / self.bunchTimeSep, 0) * self.bunchTimeSep
181  self.h1a.Fill(a)
182  self.h2a.Fill(evtNum, a)
183  # wrap-around into [0, 1] of bunch cycle
184  b = offset - round(offset / self.bunchTimeSep - 0.5, 0) * self.bunchTimeSep
185  self.h1b.Fill(b)
186  self.h2b.Fill(evtNum, b)
187 
188  def getHistogramToFit(self):
189  """
190  Selects a histogram with the peak closest to histogram center
191  """
192 
193  halfbins = self.h1a.GetNbinsX() / 2
194  if abs(self.h1a.GetMaximumBin() - halfbins) < abs(self.h1b.GetMaximumBin() - halfbins):
195  return self.h1a
196  else:
197  return self.h1b
198 
199  def terminate(self):
200  """
201  Performs the fit of the bunch finder offset distribution,
202  using a line + gaussian function.
203  The results of the fit are saved in a tree
204  """
205 
206  # Get histogram to fit
207  h_to_fit = self.getHistogramToFit()
208  maximum = h_to_fit.GetBinCenter(h_to_fit.GetMaximumBin())
209  hmax = h_to_fit.GetMaximum()
210  sigma0 = 0.16
211  xmin = h_to_fit.GetXaxis().GetXmin()
212  xmax = h_to_fit.GetXaxis().GetXmax()
213 
214  # Fit function
215  func = TF1(
216  'func',
217  '([0]/TMath::Sqrt(2*TMath::Pi() * [2]*[2])) * exp(-0.5*((x-[1])/[2])**2) \
218  + [3]*x + [4]',
219  xmin, xmax)
220  func.SetParameter(0, hmax * 0.9 * math.sqrt(2*math.pi) * sigma0)
221  func.SetParameter(1, maximum)
222  func.SetParameter(2, sigma0)
223  func.SetParLimits(2, 0.05, 0.3) # to avoid the sign ambiguity
224  func.SetParameter(3, 0.)
225  func.SetParameter(4, hmax * 0.1)
226  status = h_to_fit.Fit(func, 'L R S')
227 
228  # Tree creation, branches declaration....
229  tree = TTree('tree', '')
230 
231  expNum = array('i', [0])
232  runNum = array('i', [0])
233  fitted_offset = array('f', [0.])
234  offset = array('f', [0.])
235  offsetErr = array('f', [0.])
236  sigma = array('f', [0.])
237  integral = array('f', [0.])
238  nEvt = array('f', [0.])
239  chi2 = array('f', [0.])
240  fitStatus = array('i', [0])
241 
242  tree.Branch('expNum', expNum, 'expNum/I')
243  tree.Branch('runNum', runNum, 'runNum/I')
244  tree.Branch('fitted_offset', fitted_offset, 'fitted_offset/F')
245  tree.Branch('offset', offset, 'offset/F')
246  tree.Branch('offsetErr', offsetErr, 'offsetErr/F')
247  tree.Branch('sigma', sigma, 'sigma/F')
248  tree.Branch('chi2', chi2, 'chi2/F')
249  tree.Branch('integral', integral, 'integral/F')
250  tree.Branch('nEvt', nEvt, 'nEvt/F')
251  tree.Branch('fitStatus', fitStatus, 'fitStatus/I')
252 
253  # Dumps the fit results into the tree
254  expNum[0] = self.expNo
255  runNum[0] = int(self.run)
256  fitted_offset[0] = func.GetParameter(1)
257  new_t0 = fitted_offset[0] + self.t0
258  offset[0] = new_t0 - round(new_t0 / self.bunchTimeSep, 0) * self.bunchTimeSep
259  offsetErr[0] = func.GetParError(1)
260  sigma[0] = func.GetParameter(2)
261  integral[0] = func.GetParameter(0) / h_to_fit.GetBinWidth(1)
262  chi2[0] = func.GetChisquare() / float(func.GetNDF())
263  nEvt[0] = h_to_fit.Integral()
264  fitStatus[0] = int(status)
265 
266  tree.Fill()
267  tree.Write()
268  self.h1.Write()
269  self.h1a.Write()
270  self.h1b.Write()
271  self.h2.Write()
272  self.h2a.Write()
273  self.h2b.Write()
274  self.file.Close()
275 
276  B2RESULT('Calibration for exp' + str(self.expNo) + '-run' + self.run)
277  if self.t0Err > 0:
278  B2RESULT('Old common T0 [ns]: ' + str(round(self.t0, 3)) + ' +- ' +
279  str(round(self.t0Err, 3)))
280  else:
281  B2RESULT('Old common T0 [ns]: -- not calibrated -- ')
282  B2RESULT('New common T0 [ns]: ' + str(round(offset[0], 3)) + ' +- ' +
283  str(round(offsetErr[0], 3)))
284  B2RESULT("Output written to " + fileName)
285 
286 
287 # Database
288 use_central_database(globalTag)
289 for tag in stagingTags:
290  use_central_database(tag)
291 for db in localDB:
292  if os.path.isfile(db):
293  use_local_database(db, invertLogging=True)
294  else:
295  B2ERROR(db + ": local database not found")
296  sys.exit()
297 
298 # Create path
299 main = create_path()
300 
301 # Input (cdst files)
302 main.add_module('RootInput', inputFileNames=files)
303 
304 # Initialize TOP geometry parameters
305 main.add_module('TOPGeometryParInitializer')
306 
307 # Time Recalibrator
308 main.add_module('TOPTimeRecalibrator', subtractBunchTime=False)
309 
310 # Channel masking
311 main.add_module('TOPChannelMasker')
312 
313 # Exclude BS13d
314 main.add_module(Mask_BS13d())
315 
316 # Bunch finder
317 main.add_module('TOPBunchFinder', usePIDLikelihoods=True, subtractRunningOffset=False)
318 
319 # Common T0 calibration
320 if method == 'BF':
321  main.add_module(calibrateGlobalT0Offline())
322 elif method == 'LL':
323  main.add_module('TOPCommonT0Calibrator', sample=sampleType, outputFileName=fileName)
324 else:
325  B2ERROR('unknown method ' + method)
326  sys.exit()
327 
328 # Print progress
329 progress = register_module('Progress')
330 main.add_module(progress)
331 
332 # Process events
333 process(main)
334 
335 # Print statistics
336 print(statistics)
cdst_calibrateCommonT0.calibrateGlobalT0Offline.h1b
h1b
histogram of current offset, wrap-around into [0, 1] of bunch sep.
Definition: cdst_calibrateCommonT0.py:135
cdst_calibrateCommonT0.calibrateGlobalT0Offline.h1a
h1a
histogram of current offset, wrap-around into [-1/2, 1/2] of bunch sep.
Definition: cdst_calibrateCommonT0.py:124
cdst_calibrateCommonT0.calibrateGlobalT0Offline.h2
h2
histogram of current offset vs event number
Definition: cdst_calibrateCommonT0.py:116
cdst_calibrateCommonT0.calibrateGlobalT0Offline.initialize
def initialize(self)
Definition: cdst_calibrateCommonT0.py:96
Belle2::PyStoreObj
a (simplified) python wrapper for StoreObjPtr.
Definition: PyStoreObj.h:69
cdst_calibrateCommonT0.Mask_BS13d.event
def event(self)
Definition: cdst_calibrateCommonT0.py:68
cdst_calibrateCommonT0.calibrateGlobalT0Offline.h1
h1
histogram of current offset
Definition: cdst_calibrateCommonT0.py:113
cdst_calibrateCommonT0.calibrateGlobalT0Offline.getHistogramToFit
def getHistogramToFit(self)
Definition: cdst_calibrateCommonT0.py:188
cdst_calibrateCommonT0.calibrateGlobalT0Offline.bunchTimeSep
bunchTimeSep
bunch separation time
Definition: cdst_calibrateCommonT0.py:110
cdst_calibrateCommonT0.calibrateGlobalT0Offline.expNo
expNo
experiment number
Definition: cdst_calibrateCommonT0.py:146
cdst_calibrateCommonT0.calibrateGlobalT0Offline.h2b
h2b
histogram of current offset vs event number, wrap-around into [0, 1]
Definition: cdst_calibrateCommonT0.py:138
Belle2::PyDBObj
Class to access a DBObjPtr from Python.
Definition: PyDBObj.h:50
cdst_calibrateCommonT0.Mask_BS13d
Definition: cdst_calibrateCommonT0.py:65
cdst_calibrateCommonT0.calibrateGlobalT0Offline.t0Err
t0Err
common T0 uncertainty
Definition: cdst_calibrateCommonT0.py:154
cdst_calibrateCommonT0.calibrateGlobalT0Offline.t0
t0
common T0 used for the calibration of the input file
Definition: cdst_calibrateCommonT0.py:152
cdst_calibrateCommonT0.calibrateGlobalT0Offline.event
def event(self)
Definition: cdst_calibrateCommonT0.py:166
cdst_calibrateCommonT0.calibrateGlobalT0Offline
Definition: cdst_calibrateCommonT0.py:76
Belle2::PyStoreArray
a (simplified) python wrapper for StoreArray.
Definition: PyStoreArray.h:58
cdst_calibrateCommonT0.calibrateGlobalT0Offline.terminate
def terminate(self)
Definition: cdst_calibrateCommonT0.py:199
cdst_calibrateCommonT0.calibrateGlobalT0Offline.run
run
run number formatted as string with leading zeros
Definition: cdst_calibrateCommonT0.py:149
cdst_calibrateCommonT0.calibrateGlobalT0Offline.h2a
h2a
histogram of current offset vs event number, wrap-around into [-1/2, 1/2]
Definition: cdst_calibrateCommonT0.py:127
cdst_calibrateCommonT0.calibrateGlobalT0Offline.file
file
Output root file.
Definition: cdst_calibrateCommonT0.py:103