13 from ROOT
import Belle2
14 from ROOT
import TH1F, TH2F, TFile, TTree, TF1
15 from array
import array
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'
27 stagingTags = [
'staging_data_reprocessing']
30 output_dir =
'commonT0'
37 print(
"usage: basf2", argvs[0],
"experiment run")
39 experiment = int(argvs[1])
42 expNo =
'e' +
'{:0=4d}'.format(experiment)
43 runNo =
'r' +
'{:0=5d}'.format(run)
47 for typ
in [
'4S',
'Continuum',
'Scan']:
48 folder = data_dir +
'/' + expNo +
'/' + typ +
'/' + runNo +
'/' + skim_dir
49 files += glob.glob(folder +
'/cdst.*.root')
51 B2ERROR(
'No cdst files found')
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)
61 fileName = output_folder +
'/commonT0-' + expNo +
'-' + runNo +
'.root'
62 print(
'Output file:', fileName)
66 ''' exclude (mask-out) BS 13d '''
69 ''' event processing '''
72 if digit.getModuleID() == 13
and digit.getBoardstackNumber() == 3:
73 digit.setHitQuality(Belle2.TOPDigit.c_Junk)
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.
87 * The histogram of the bunch finder time offset
89 * The fit parameters ina tree format (usefuly when merging several files)
98 Creates the histogram used for the commoT0 calculation, and
99 takes the necessary objects from the DataStore
103 self.
file = TFile.Open(fileName,
'recreate')
106 if not geo.isValid():
107 B2FATAL(
'TOP geometry not available in database')
113 self.
h1 = TH1F(
"offset",
"current offset; offset [ns]", 600, -9.0, 9.0)
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]")
124 self.
h1a = TH1F(
"offset_a",
"current offset; offset [ns]", 200, xmi, xma)
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]")
135 self.
h1b = TH1F(
"offset_b",
"current offset; offset [ns]", 200, xmi, xma)
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]")
146 self.
expNo = evtMetaData.getExperiment()
149 self.
run =
'{:0=5d}'.format(evtMetaData.getRun())
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))
162 B2INFO(
'No common T0 calibration done yet')
164 B2ERROR(
'Common T0 not available in database')
168 Fills the histogram of the time offset
172 if not recBunch.isValid():
174 if recBunch.isReconstructed()
and recBunch.getNumTracks() == 2:
175 offset = recBunch.getCurrentOffset()
176 evtNum = evtMetaData.getEvent()
178 self.
h2.Fill(evtNum, offset)
182 self.
h2a.Fill(evtNum, a)
186 self.
h2b.Fill(evtNum, b)
190 Selects a histogram with the peak closest to histogram center
193 halfbins = self.
h1a.GetNbinsX() / 2
194 if abs(self.
h1a.GetMaximumBin() - halfbins) < abs(self.
h1b.GetMaximumBin() - halfbins):
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
208 maximum = h_to_fit.GetBinCenter(h_to_fit.GetMaximumBin())
209 hmax = h_to_fit.GetMaximum()
211 xmin = h_to_fit.GetXaxis().GetXmin()
212 xmax = h_to_fit.GetXaxis().GetXmax()
217 '([0]/TMath::Sqrt(2*TMath::Pi() * [2]*[2])) * exp(-0.5*((x-[1])/[2])**2) \
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)
224 func.SetParameter(3, 0.)
225 func.SetParameter(4, hmax * 0.1)
226 status = h_to_fit.Fit(func,
'L R S')
229 tree = TTree(
'tree',
'')
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])
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')
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
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)
276 B2RESULT(
'Calibration for exp' + str(self.
expNo) +
'-run' + self.
run)
278 B2RESULT(
'Old common T0 [ns]: ' + str(round(self.
t0, 3)) +
' +- ' +
279 str(round(self.
t0Err, 3)))
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)
288 use_central_database(globalTag)
289 for tag
in stagingTags:
290 use_central_database(tag)
292 if os.path.isfile(db):
293 use_local_database(db, invertLogging=
True)
295 B2ERROR(db +
": local database not found")
302 main.add_module(
'RootInput', inputFileNames=files)
305 main.add_module(
'TOPGeometryParInitializer')
308 main.add_module(
'TOPTimeRecalibrator', subtractBunchTime=
False)
311 main.add_module(
'TOPChannelMasker')
317 main.add_module(
'TOPBunchFinder', usePIDLikelihoods=
True, subtractRunningOffset=
False)
323 main.add_module(
'TOPCommonT0Calibrator', sample=sampleType, outputFileName=fileName)
325 B2ERROR(
'unknown method ' + method)
329 progress = register_module(
'Progress')
330 main.add_module(progress)