SND@LHC Software
Loading...
Searching...
No Matches
ConvRawData.py
Go to the documentation of this file.
1
2#!/usr/bin/env python
3import ROOT,os,sys
4import boardMappingParser
5import csv
6import time
7import calendar
8from rootpyPickler import Unpickler
9import json
10
11# raw data from Ettore: https://cernbox.cern.ch/index.php/s/Ten7ilKuD3qdnM2
12
13class ConvRawDataPY(ROOT.FairTask):
14 " raw data conversion by event "
15 def Init(self,options):
16 self.debug = options.debug
17 self.monitoring = options.online
18 self.auto = options.auto
19 local = False
20 if (options.path.find('eos')<0 and options.server.find('root://')<0 and not options.online) or os.path.isdir(options.path): local = True
21 if local: server = ''
22 else:
23 server = options.server
24 from XRootD import client
25 # https://xrootd.slac.stanford.edu/doc/python/xrootd-python-0.1.0/examples.html
26 if options.runNumber>0:
27 runNr = str(options.runNumber).zfill(6)
28 path = options.path+'run_'+ runNr+'/'
29 if options.partition < 0:
30 inFile = 'data.root'
31 self.outFile = "sndsw_raw_"+runNr+'.root'
32 else:
33 part = str(options.partition).zfill(4)
34 inFile = 'data_'+part+'.root'
35 self.outFile = "sndsw_raw_"+runNr+'-'+part+'.root'
36 if self.monitoring: # this is an online stream
37 server = options.server
38 runNr = str( abs(options.runNumber) ).zfill(6)
39 path = options.path+'run_'+ runNr+'/'
40 inFile = 'data_'+part+'.root'
41 self.outFile = ROOT.TMemFile('monitorRawData', 'recreate')
42
43# get filling scheme per run
44 self.fsdict = False
45 try:
46 if options.path.find('2022')!=-1: fpath = "/eos/experiment/sndlhc/convertedData/physics/2022/"
47 elif options.path.find('2023')!=-1: fpath = "/eos/experiment/sndlhc/convertedData/physics/2023/"
48 elif options.path.find('2024')!=-1: fpath = "/eos/experiment/sndlhc/convertedData/physics/2024/"
49 else: fpath = "/eos/experiment/sndlhc/convertedData/commissioning/TI18/"
50 fg = ROOT.TFile.Open(options.server+fpath+"/FSdict.root")
51 pkl = Unpickler(fg)
52 FSdict = pkl.load('FSdict')
53 fg.Close()
54 if options.runNumber in FSdict:
55 if 'B1' in FSdict[options.runNumber]: self.fsdict = FSdict[options.runNumber]
56 except:
57 print('continue without knowing filling scheme',options.server+options.path)
58
59 # put the run's FS in format to be passed to FairTasks as input
60 self.FSmap = ROOT.TMap()
61 if self.fsdict:
62 # For LHC ion runs the phase shift of beam 2 could not be
63 # determined and is assigned 0 or 1718(=1782-64) in the FS dict
64 # For pp runs phase2 is 3564-129.
65 if self.fsdict['phaseShift2'] == 0 or self.fsdict['phaseShift2'] == 1718:
66 Nbunches = 1782
67 else: Nbunches = 3564 # proton runs
68 for bunchNumber in range (0, Nbunches):
69 nb1 = (Nbunches + bunchNumber - self.fsdict['phaseShift1'])%Nbunches
70 nb2 = (Nbunches + bunchNumber - self.fsdict['phaseShift1']- self.fsdict['phaseShift2'])%Nbunches
71 b1 = nb1 in self.fsdict['B1']
72 b2 = nb2 in self.fsdict['B2']
73 IP1 = False
74 IP2 = False
75 if b1:
76 IP1 = self.fsdict['B1'][nb1]['IP1']
77 if b2:
78 IP2 = self.fsdict['B2'][nb2]['IP2']
79 self.FSmap.Add(ROOT.TObjString(str(bunchNumber)), ROOT.TObjString(str(int(IP2))+str(int(IP1))+str(int(b2))+str(int(b1))))
80 else: self.FSmap.Add(ROOT.TObjString("0"), ROOT.TObjString("-1"))
81
82 self.run = ROOT.FairRunAna()
83 self.ioman = ROOT.FairRootManager.Instance()
84 ioman = self.ioman
85# open input file
86 if options.auto:
87 mpath = options.path+options.monitorTag+'run_'+ runNr+'/'
88 self.fiN=ROOT.TFile.Open(server+mpath+inFile)
89 else:
90 self.fiN=ROOT.TFile.Open(server+path+inFile)
91 self.nEvents = 1
92 self.newFormat = True
93 if self.fiN.Get('event'): self.newFormat = False # old format
94 if not self.monitoring:
95 if self.newFormat:
96 if options.nEvents<0: self.nEvents = self.fiN.data.GetEntries()
97 else: self.nEvents = min(options.nEvents,self.fiN.data.GetEntries())
98 else:
99 if options.nEvents<0: self.nEvents = self.fiN.event.GetEntries()
100 else: self.nEvents = min(options.nEvents,self.fiN.event.GetEntries())
101 print('converting ',self.nEvents,' events ',' of run',options.runNumber)
102 # Pass input parameters to the task - runN, paths, etc.
103 ioman.RegisterInputObject('runN', ROOT.TObjString(str(options.runNumber)))
104 ioman.RegisterInputObject('pathCalib', ROOT.TObjString(server+path))
105 ioman.RegisterInputObject('pathJSON', ROOT.TObjString(server+path))
106 ioman.RegisterInputObject('nEvents', ROOT.TObjString(str(self.nEvents)))
107 ioman.RegisterInputObject('nStart', ROOT.TObjString(str(options.nStart)))
108 ioman.RegisterInputObject('debug', ROOT.TObjString(str(options.debug)))
109 ioman.RegisterInputObject('stop', ROOT.TObjString(str(options.stop)))
110 ioman.RegisterInputObject('heartBeat', ROOT.TObjString(str(options.heartBeat)))
111 ioman.RegisterInputObject('makeCalibration', ROOT.TObjString(str(int(options.makeCalibration))))
112 ioman.RegisterInputObject('chi2Max', ROOT.TObjString(str(options.chi2Max)))
113 ioman.RegisterInputObject('saturationLimit', ROOT.TObjString(str(options.saturationLimit)))
114 ioman.RegisterInputObject('local', ROOT.TObjString(str(int(local))))
115 ioman.RegisterInputObject('newFormat', ROOT.TObjString(str(int(self.newFormat))))
116 ioman.RegisterInputObject('FSmap', self.FSmap)
117 self.options = options
118
119 # Initialize logger: set severity and verbosity
120 logger = ROOT.FairLogger.GetLogger()
121 logger.SetColoredLog(True)
122 logger.SetLogVerbosityLevel('low')
123 logger.SetLogScreenLevel('error')
124 logger.SetLogToScreen(True)
125 if options.debug:
126 logger.SetLogToFile(True)
127 logger.SetLogFileName('run_'+str(options.runNumber)+'-'+part+'_CPP.log')
128 logger.SetLogFileLevel('info')
129 logger.SetLogScreenLevel('info')
130
131 # Pass raw data file as input object
132 ioman.RegisterInputObject("rawData", self.fiN)
133
134 # Set output
135 if self.monitoring: self.outfile = ROOT.FairRootFileSink(self.outFile)
136 elif options.FairTask_convRaw: self.outfile = ROOT.FairRootFileSink(self.outFile.replace('.root','_CPP.root'))
137 else: self.outfile = ROOT.FairRootFileSink(self.outFile)
138 self.run.SetSink(self.outfile)
139
140 # Don't use FairRoot's default event header settings
141 self.run.SetEventHeaderPersistence(False)
142 self.xrdb = ROOT.FairRuntimeDb.instance()
143 self.xrdb.getContainer("FairBaseParSet").setStatic()
144 self.xrdb.getContainer("FairGeoParSet").setStatic()
145# Fair convRawData task
146 if options.FairTask_convRaw:
147 self.run.AddTask(ROOT.ConvRawData())
148 self.fSink = self.outfile
149 #X = fSink.GetRootFile()
150 #X.SetCompressionLevel(ROOT.CompressionSettings(ROOT.kLZMA, 5)) # it seems it has no effect
151 self.run.Init()
152
153#-------end of init for cpp ------------------------------------
154 else:
155 if options.makeCalibration:
156 if local:
157 fqdc_cal = open(path+'/qdc_cal.csv')
158 Lqdc = fqdc_cal.readlines()
159 ftdc_cal = open(path+'/tdc_cal.csv')
160 Ltdc = ftdc_cal.readlines()
161 else:
162 with client.File() as f:
163 f.open(server+path+"/qdc_cal.csv")
164 status, L = f.read()
165 Lqdc = L.decode().split('\n')
166 f.close()
167 f.open(server+path+"/tdc_cal.csv")
168 status, L = f.read()
169 Ltdc = L.decode().split('\n')
170 f.close()
171 # calibration data
172 self.qdc_cal = {}
173 qdc_cal = self.qdc_cal
174 L = Lqdc
175 for l in range(1,len(L)):
176 tmp = L[l].replace('\n','').split(',')
177 if len(tmp)<10:continue
178 board_id = int(tmp[0])
179 if not board_id in qdc_cal: qdc_cal[board_id]={}
180 fe_id = int(tmp[1])
181 if not fe_id in qdc_cal[board_id]: qdc_cal[board_id][fe_id]={}
182 channel = int(tmp[2])
183 if not channel in qdc_cal[board_id][fe_id]: qdc_cal[board_id][fe_id][channel]={}
184 tac = int(tmp[3])
185 if not tac in qdc_cal[board_id][fe_id][channel]: qdc_cal[board_id][fe_id][channel][tac]={}
186 X = qdc_cal[board_id][fe_id][channel][tac]
187 X['a']=float(tmp[4])
188 X['b']=float(tmp[5])
189 X['c']=float(tmp[6])
190 X['d']=float(tmp[8])
191 X['e']=float(tmp[10])
192 if float(tmp[9]) < 2: X['chi2Ndof'] = 999999.
193 else: X['chi2Ndof']=float(tmp[7])/float(tmp[9])
194
195 L=Ltdc
196 for l in range(1,len(L)):
197 tmp = L[l].replace('\n','').split(',')
198 if len(tmp)<9:continue
199 board_id = int(tmp[0])
200 if not board_id in qdc_cal: qdc_cal[board_id]={}
201 fe_id = int(tmp[1])
202 if not fe_id in qdc_cal[board_id]: qdc_cal[board_id][fe_id]={}
203 channel = int(tmp[2])
204 if not channel in qdc_cal[board_id][fe_id]: qdc_cal[board_id][fe_id][channel]={}
205 tac = int(tmp[3])
206 if not tac in qdc_cal[board_id][fe_id][channel]: qdc_cal[board_id][fe_id][channel][tac]={}
207 tdc = int(tmp[4])
208 if not tdc in qdc_cal[board_id][fe_id][channel][tac]: qdc_cal[board_id][fe_id][channel][tac][tdc]={}
209 X = qdc_cal[board_id][fe_id][channel][tac][tdc]
210 X['a']=float(tmp[5])
211 X['b']=float(tmp[6])
212 X['c']=float(tmp[7])
213 X['d']=float(tmp[9])
214 if float(tmp[10]) < 2: X['chi2Ndof'] = 999999.
215 else: X['chi2Ndof']=float(tmp[8])/float(tmp[10])
216
217 # read MuFilter mappings of SiPM channel to tofpet_id and tofpet_channel
218 SiPMmap={}
219 self.TofpetMap = {}
220 key = {'DS':2,'US':1,'Veto':0}
221 sndswPath = os.environ['SNDSW_ROOT']
222 for system in key:
223 infile = sndswPath+"/geometry/"+system+"_SiPM_mapping.csv"
224 SiPMmap[system] = {}
225 if os.path.isdir(sndswPath+"/geometry"):
226 with open(infile, "r") as f:
227 reader = csv.DictReader(f, delimiter=',')
228 for i, row in enumerate(reader):
229 SiPMmap[system][int(row['SiPM'])] = [int(row['CN']), int(row['pin']), int(row['TOFPET']), int(row['CHANNEL'])]
230 self.TofpetMap[key[system]] = {}
231 for channel in SiPMmap[system]:
232 row = SiPMmap[system][channel]
233 self.TofpetMap[key[system]][row[2]*1000+row[3]] = channel
234 if local:
235 with open(path+'board_mapping.json') as f:
236 jsonStr = f.read()
237 else:
238 with client.File() as f:
239 f.open(server+path+"/board_mapping.json")
240 status, jsonStr = f.read()
241 # pass the read string to getBoardMapping()
242 self.boardMaps = boardMappingParser.getBoardMapping(jsonStr)
243
244 self.slots = {0:'A',1:'A',2:'B',3:'B',4:'C',5:'C',6:'D',7:'D'}
245
246 self.MufiSystem = {}
247 for b in self.boardMaps['MuFilter']:
248 board_id = int(b.split('_')[1])
249 self.MufiSystem[board_id]={}
250 for x in self.boardMaps['MuFilter'][b]:
251 for slot in self.slots:
252 s = 0
253 tmp = self.boardMaps['MuFilter'][b][x].split('_')[0]
254 if tmp=='US': s = 1
255 elif tmp=='DS': s = 2
256 if self.slots[slot]==x: self.MufiSystem[board_id][slot]=s
257
258 self.offMap={}
259 # first bar, number of sipm channels / bar and direction
260 for s in range(1,3):
261 for o in ['Left','Right']:
262 self.offMap['Veto_'+str(s)+o] =[10000 + (s-1)*1000+ 0,8,2] # first channel, nSiPMs, nSides, from bottom to top
263 self.offMap['Veto_3Vert'] = [10000 + 2*1000+ 6,-8,1]
264 for s in range(1,6):
265 for o in ['Left','Right']:
266 self.offMap['US_'+str(s)+o] =[20000 + (s-1)*1000+ 9,-8,2] # from top to bottom
267 for s in range(1,5):
268 for o in ['Vert']:
269 self.offMap['DS_'+str(s)+o] =[30000 + (s-1)*1000+ 119, -1,1] # direction not known
270 if s>3: continue
271 for o in ['Left','Right']:
272 self.offMap['DS_'+str(s)+o] =[30000 + (s-1)*1000+ 59,-1,2] # direction not known
273
274 self.boards = {}
275 for b in self.fiN.GetListOfKeys():
276 name = b.GetName()
277 if name.find('board')!=0: continue
278 self.boards[name]=self.fiN.Get(name)
279
280 self.fSink = self.outfile
281 self.sTree = ROOT.TTree('rawConv','raw data converted')
282 ROOT.gDirectory.pwd()
283 self.header = ROOT.SNDLHCEventHeader()
284 eventSND = self.sTree.Branch("EventHeader.",self.header,32000,-1)
285
286 self.digiSciFi = ROOT.TClonesArray("sndScifiHit")
287 self.digiSciFiBranch = self.sTree.Branch("Digi_ScifiHits",self.digiSciFi,32000,1)
288 self.digiMuFilter = ROOT.TClonesArray("MuFilterHit")
289 self.digiMuFilterHitBranch = self.sTree.Branch("Digi_MuFilterHits",self.digiMuFilter,32000,1)
290 self.geoVersion = 0
291
292 B = ROOT.TList()
293 B.SetName('BranchList')
294 B.Add(ROOT.TObjString('sndScifiHit'))
295 B.Add(ROOT.TObjString('MuFilterHit'))
296 B.Add(ROOT.TObjString('SNDLHCEventHeader'))
297 self.fSink.WriteObject(B,"BranchList", ROOT.TObject.kSingleKey)
298 self.fSink.SetRunId(options.runNumber)
299 self.fSink.SetOutTree(self.sTree)
300
301 if self.newFormat: self.run_startUTC = self.getStartTime()
302
303#-------end of init for py ------------------------------------
304
305 def channel(self,tofpet_id,tofpet_channel,position):
306 return (64 * tofpet_id + 63 - tofpet_channel + 512*position) # 512 channels per mat, 1536 channels per plane
307 # one channel covers all 6 layers of fibres
308
309 # following procedure in https://gitlab.cern.ch/scifi-telescope/scifi-rec/-/blob/ettore-tofpet/lib/processors/applycalibration.cpp
310 def qdc_chi2(self,board_id,tofpet_id,channel,tac,TDC=0):
311 qdc_cal = self.qdc_cal
312 par = qdc_cal[board_id][tofpet_id][channel][tac]
313 parT = qdc_cal[board_id][tofpet_id][channel][tac][TDC]
314 return max(par['chi2Ndof'],parT['chi2Ndof'])
315
316 def qdc_sat(self,board_id,tofpet_id,channel,tac,v_fine):
317 qdc_cal = self.qdc_cal
318 par = qdc_cal[board_id][tofpet_id][channel][tac]
319 return v_fine/par['d']
320
321 def comb_calibration(self,board_id,tofpet_id,channel,tac,v_coarse,v_fine,t_coarse,t_fine,GQDC = 1.0, TDC=0): # max gain QDC = 3.6
322 qdc_cal = self.qdc_cal
323 par = qdc_cal[board_id][tofpet_id][channel][tac]
324 parT = par[TDC]
325 x = t_fine
326 ftdc = (-parT['b']-ROOT.TMath.Sqrt(parT['b']**2-4*parT['a']*(parT['c']-x)))/(2*parT['a']) # Ettore 28/01/2022 +parT['d']
327 timestamp = t_coarse + ftdc
328 tf = timestamp - t_coarse
329 x = v_coarse - tf
330 fqdc = -par['c']*ROOT.TMath.Log(1+ROOT.TMath.Exp( par['a']*(x-par['e'])**2-par['b']*(x-par['e']) )) + par['d']
331 value = (v_fine-fqdc)/GQDC
332 return timestamp,value,max(par['chi2Ndof'],parT['chi2Ndof']),v_fine/par['d']
333
334 def calibrationReport(self):
335 qdc_cal = self.qdc_cal
336 report = {}
337 TDC = 0
338 for b in qdc_cal:
339 for t in qdc_cal[b]:
340 for c in qdc_cal[b][t]:
341 for tac in qdc_cal[b][t][c]:
342 par = qdc_cal[b][t][c][tac]
343 if 'chi2Ndof' in par: chi2 = par['chi2Ndof']
344 else: chi2=-1
345 parT = qdc_cal[b][t][c][tac][TDC]
346 if 'chi2Ndof' in parT: chi2T = parT['chi2Ndof']
347 else: chi2T=-1
348 key = tac +10*c + t*10*100 + b*10*100*100
349 if not key in report: report[key]=[chi2,chi2T]
350 return report
351
352 def debugMapping(self,board,tofpet_id,tofpet_channel):
353 key = (tofpet_id%2)*1000 + tofpet_channel
354 tmp = boardMaps['MuFilter'][board][self.slots[tofpet_id]]
355 sipmChannel = TofpetMap[system][key]-1
356 nSiPMs = abs(offMap[tmp][1])
357 nSides = abs(offMap[tmp][2])
358 direction = int(offMap[tmp][1]/nSiPMs)
359 detID = offMap[tmp][0] + direction*(sipmChannel//(nSiPMs))
360 sipm_number = sipmChannel%(nSiPMs)
361 print(sipmChannel,nSiPMs,nSides,detID,sipm_number)
362
363 def printChannelInfo(self):
364 runNr = str( abs(self.options.runNumber) ).zfill(6)
365 path = self.options.path+'run_'+ runNr+'/'
366 from XRootD import client
367 with client.File() as f:
368 f.open(self.options.server+path+"/channels.csv")
369 status, L = f.read()
370 self.Lchannel = L.decode().split('\n')
371 f.close()
372 for l in self.Lchannel: print(l)
373
374 def printCurrentRun(self):
375 filename = "currently_processed_file.txt"
376 path = self.options.path
377 from XRootD import client
378 with client.File() as f:
379 f.open(self.options.server+path+filename)
380 status, L = f.read()
381 self.Lcrun = L.decode().split('\n')
382 f.close()
383 for l in self.Lcrun: print(l)
384
385 def getStartTime(self):
386 runNr = str( abs(self.options.runNumber) ).zfill(6)
387 path = self.options.path+'run_'+ runNr+'/'
388 jname = "run_timestamps.json"
389 from XRootD import client
390 with client.File() as f:
391 f.open(self.options.server+path+jname)
392 status, jsonStr = f.read()
393 f.close()
394 date = json.loads(jsonStr)
395 time_str = date['start_time'].replace('Z','')
396 time_obj = time.strptime(time_str, '%Y-%m-%dT%H:%M:%S')
397 startTimeOfRun = calendar.timegm(time_obj)
398 return startTimeOfRun
399
400 # reading hits and converting to event information
401
402 # https://gitlab.cern.ch/snd-scifi/software/-/wikis/Raw-data-format
403 # tac: 0-3, identifies the Time-To-Analogue converter used for this hit (each channel has four of them and they require separate calibration).
404 # t_coarse: Coarse timestamp of the hit, running on a 4 times the LHC clock
405 # t_fine: 0-1023, fine timestamp of the hit. It contains the raw value of the charge digitized by the TDC and requires calibration.
406 # v_coarse: 0-1023, QDC mode: it represents the number of clock cycles the charge integration lasted.
407 # v_fine = 0-1023, QDC mode: represents the charge measured. Requires calibration.
408
409 def Run(self):
410 if self.options.FairTask_convRaw:
411 self.run.Run(self.options.nStart, self.nEvents)
412 else:
413 for eventNumber in range(self.options.nStart,self.nEvents):
414 self.executeEvent(eventNumber)
415 # fill TTree
416 self.sTree.Fill()
417 def executeEvent(self,eventNumber):
418 if self.newFormat: self.executeEvent1(eventNumber)
419 else: self.executeEvent0(eventNumber)
420 def executeEvent1(self,eventNumber):
421 if self.options.FairTask_convRaw:
422 # update source and run conversion starting from a custom eventN
423 if self.auto:
424 self.run.GetTask("ConvRawData").UpdateInput(eventNumber)
425 self.run.Run(self.options.nStart, self.nEvents)
426 Fout = self.outfile.GetRootFile()
427 self.sTree = Fout.Get('cbmsim')
428 return
429 if eventNumber%self.options.heartBeat==0 or self.debug:
430 print('run ',self.options.runNumber, ' event',eventNumber," ",time.ctime())
431
432 event = self.fiN.data
433 event.GetEvent(eventNumber)
434 self.header.SetEventTime(event.evt_timestamp)
435 self.header.SetUTCtimestamp(int(event.evt_timestamp*6.23768*1e-9 + self.run_startUTC))
436 self.header.SetEventNumber(event.evt_number) # for new event header
437 self.header.SetFlags(event.evt_flags)
438 self.header.SetRunId( self.options.runNumber )
439 if self.FSmap.GetEntries()>1:
440 if self.header.GetAccMode()== 12: # ion runs
441 self.header.SetBunchType(int(str(self.FSmap.GetValue(str(int((event.evt_timestamp%(4*3564))/8))))))
442 else: # proton runs
443 self.header.SetBunchType(int(str(self.FSmap.GetValue(str(int((event.evt_timestamp%(4*3564))/4))))))
444 else:
445 self.header.SetBunchType(int(str(self.FSmap.GetValue("0"))))
446
447 indexSciFi=0
448 self.digiSciFi.Delete()
449 digiSciFiStore = {}
450 indexMuFilter = 0
451 self.digiMuFilter.Delete()
452 digiMuFilterStore = {}
453
454 for n in range(event.n_hits):
455 board_id = event.boardId[n]
456 board = 'board_'+str(board_id)
457 scifi = True
458 if board in self.boardMaps['Scifi']:
459 station,mat = self.boardMaps['Scifi'][board]
460 elif board in self.boardMaps['MuFilter']:
461 scifi = False
462 else:
463 print(board+' not known. Serious error, stop')
464 1/0
465 mask = False
466 tofpet_id = event.tofpetId[n]
467 tofpet_channel = event.tofpetChannel[n]
468 if self.options.makeCalibration:
469 if self.options.debug:
470 tac = event.tac[n]
471 print(scifi,board_id,tofpet_id,tofpet_channel,tac,event.tCoarse[n],event.tFine[n],event.vCoarse[n],event.vFine[n])
472 TDC,QDC,Chi2ndof,satur = self.comb_calibration(board_id,tofpet_id,tofpet_channel,tac,event.vCoarse[n],event.vFine[n],event.tCoarse[n],event.tFine[n])
473 else:
474 TDC = event.timestamp[n]
475 QDC = event.value[n]
476 QDCchi2N = event.valueCalChi2[n]/event.valueCalDof[n]
477 TDCchi2N = event.timestampCalChi2[n]/event.timestampCalDof[n]
478 Chi2ndof = max( QDCchi2N,TDCchi2N )
479 satur = event.value_saturation[n]
480 if Chi2ndof > self.options.chi2Max:
481 if QDC>1E20: QDC = 997. # checking for inf
482 if QDC != QDC: QDC = 998. # checking for nan
483 if QDC>0: QDC = -QDC
484 mask = True
485 elif satur > self.options.saturationLimit or QDC>1E20 or QDC != QDC:
486 if QDC>1E20: QDC = 987. # checking for inf
487 if self.options.debug:
488 print('inf',board_id,tofpet_id,tofpet_channel,event.tac[n],event.vCoarse[n],event.vFine[n],TDC-event.tCoarse[n],eventNumber,Chi2ndof)
489 if QDC != QDC: QDC = 988. # checking for nan
490 if self.options.debug:
491 print('nan',board_id,tofpet_id,tofpet_channel,event.tac[n],event.vCoarse[n],event.vFine[n],TDC-event.tCoarse[n],eventNumber,Chi2ndof)
492 A = int(min( QDC,1000))
493 B = min(satur,999)/1000.
494 QDC = A+B
495 mask = True
496 elif Chi2ndof > self.options.chi2Max:
497 if QDC>0: QDC = -QDC
498 mask = True
499 if self.options.debug:
500 print('calibrated: tdc = ',TDC,' qdc = ',QDC) # TDC clock cycle = 160 MHz or 6.25ns
501
502 if not scifi:
503 # mufi encoding
504 system = self.MufiSystem[board_id][tofpet_id]
505 key = (tofpet_id%2)*1000 + tofpet_channel
506 tmp = self.boardMaps['MuFilter'][board][self.slots[tofpet_id]]
507 if self.options.debug or not tmp.find('not')<0: print('debug',tmp,system,key,board,tofpet_id,tofpet_id%2,tofpet_channel)
508 sipmChannel = 99
509 if not key in self.TofpetMap[system]:
510 print('key does not exist',key)
511 print(system, key,board,tofpet_id, self.TofpetMap[system])
512 else:
513 sipmChannel = self.TofpetMap[system][key]-1
514 nSiPMs = abs(self.offMap[tmp][1])
515 nSides = abs(self.offMap[tmp][2])
516 direction = int(self.offMap[tmp][1]/nSiPMs)
517 detID = self.offMap[tmp][0] + direction*(sipmChannel//(nSiPMs))
518 sipm_number = sipmChannel%(nSiPMs)
519 if tmp.find('Right')>0:
520 sipm_number += nSiPMs
521 if not detID in digiMuFilterStore:
522 digiMuFilterStore[detID] = ROOT.MuFilterHit(detID,nSiPMs,nSides)
523 test = digiMuFilterStore[detID].GetSignal(sipm_number)
524 digiMuFilterStore[detID].SetDigi(QDC,TDC,sipm_number)
525 digiMuFilterStore[detID].SetDaqID(sipm_number,n,board_id, tofpet_id, tofpet_channel)
526 if mask: digiMuFilterStore[detID].SetMasked(sipm_number)
527
528 if self.options.debug:
529 print('create mu hit: ',detID,tmp,system,tofpet_id,self.offMap[tmp],sipmChannel,nSiPMs,nSides,test)
530 print(' ',detID,sipm_number,QDC,TDC)
531 if test>0 or detID%1000>200 or sipm_number>15:
532 print('what goes wrong?',detID,sipm_number,system,key,board,tofpet_id,tofpet_channel,test)
533
534 else:
535 # scifi encoding
536 chan = self.channelchannel(tofpet_id,tofpet_channel,mat)
537 orientation = 1
538 if station[2]=="Y": orientation = 0
539 sipmLocal = (chan - mat*512)
540 sipmID = 1000000*int(station[1]) + 100000*orientation + 10000*mat + 1000*(sipmLocal//128) + chan%128
541 if not sipmID in digiSciFiStore:
542 digiSciFiStore[sipmID] = ROOT.sndScifiHit(sipmID)
543 digiSciFiStore[sipmID].SetDigi(QDC,TDC)
544 digiSciFiStore[sipmID].SetDaqID(0,n, board_id, tofpet_id, tofpet_channel)
545 if mask: digiSciFiStore[sipmID].setInvalid()
546 if self.options.debug:
547 print('create scifi hit: tdc = ',board,sipmID,QDC,TDC)
548 print('tofpet:', tofpet_id,tofpet_channel,mat,chan)
549 print(station[1],station[2],mat,chan,int(chan/128)%4,chan%128)
550
551# copy hits to detector branches
552 for sipmID in digiSciFiStore:
553 if self.digiSciFi.GetSize() == indexSciFi: self.digiSciFi.Expand(indexSciFi+100)
554 self.digiSciFi[indexSciFi]=digiSciFiStore[sipmID]
555 indexSciFi+=1
556 for detID in digiMuFilterStore:
557 if self.digiMuFilter.GetSize() == indexMuFilter: self.digiMuFilter.Expand(indexMuFilter+100)
558 self.digiMuFilter[indexMuFilter]=digiMuFilterStore[detID]
559 indexMuFilter+=1
560 def executeEvent0(self,eventNumber):
561 if self.options.FairTask_convRaw:
562 # update source and run conversion starting from a custom eventN
563 if self.auto:
564 self.run.GetTask("ConvRawData").UpdateInput(eventNumber)
565 self.run.Run(self.options.nStart, self.nEvents)
566 Fout = self.outfile.GetRootFile()
567 self.sTree = Fout.Get('cbmsim')
568 return
569 if eventNumber%self.options.heartBeat==0 or self.debug:
570 print('run ',self.options.runNumber, ' event',eventNumber," ",time.ctime())
571 event = self.fiN.event
572 rc = event.GetEvent(eventNumber)
573 self.header.SetEventTime(event.timestamp)
574 self.header.SetRunId( self.options.runNumber )
575 self.header.SetInputFileId(self.geoVersion)
576
577 indexSciFi=0
578 self.digiSciFi.Delete()
579 digiSciFiStore = {}
580 indexMuFilter = 0
581 self.digiMuFilter.Delete()
582 digiMuFilterStore = {}
583
584 for board in self.boards:
585 board_id = int(board.split('_')[1])
586 scifi = True
587 if board in self.boardMaps['Scifi']:
588 station,mat = self.boardMaps['Scifi'][board]
589 elif board in self.boardMaps['MuFilter']:
590 scifi = False
591 else:
592 print(board+' not known. Serious error, stop')
593 1/0
594 bt = self.boards[board]
595 rc = bt.GetEvent(eventNumber)
596 for n in range(bt.n_hits):
597 mask = False
598 tofpet_id = bt.tofpet_id[n]
599 tofpet_channel = bt.tofpet_channel[n]
600 if self.options.makeCalibration:
601 if self.options.debug:
602 print(scifi,board_id,bt.tofpet_id[n],bt.tofpet_channel[n],bt.tac[n],bt.t_coarse[n],bt.t_fine[n],bt.v_coarse[n],bt.v_fine[n])
603 tac = bt.tac[n]
604 TDC,QDC,Chi2ndof,satur = self.comb_calibration(board_id,tofpet_id,tofpet_channel,tac,bt.v_coarse[n],bt.v_fine[n],bt.t_coarse[n],bt.t_fine[n])
605 else:
606 TDC = bt.timestamp[n]
607 QDC = bt.value[n]
608 Chi2ndof = 1
609 satur = 0.
610 if Chi2ndof > self.options.chi2Max:
611 if QDC>1E20: QDC = 997. # checking for inf
612 if QDC != QDC: QDC = 998. # checking for nan
613 if QDC>0: QDC = -QDC
614 mask = True
615 elif satur > self.options.saturationLimit or QDC>1E20 or QDC != QDC:
616 if QDC>1E20: QDC = 987. # checking for inf
617 if self.options.debug:
618 print('inf',board_id,bt.tofpet_id[n],bt.tofpet_channel[n],bt.tac[n],bt.v_coarse[n],bt.v_fine[n],TDC-bt.t_coarse[n],eventNumber,Chi2ndof)
619 if QDC != QDC: QDC = 988. # checking for nan
620 if self.options.debug:
621 print('nan',board_id,bt.tofpet_id[n],bt.tofpet_channel[n],bt.tac[n],bt.v_coarse[n],bt.v_fine[n],\
622 TDC-bt.t_coarse[n],TDC,bt.t_coarse[n],eventNumber,Chi2ndof)
623 A = int(min( QDC,1000))
624 B = min(satur,999)/1000.
625 QDC = A+B
626 mask = True
627 elif Chi2ndof > self.options.chi2Max:
628 if QDC>0: QDC = -QDC
629 mask = True
630 if self.options.debug:
631 print('calibrated: tdc = ',TDC,' qdc = ',QDC) # TDC clock cycle = 160 MHz or 6.25ns
632
633 if not scifi:
634 # mufi encoding
635 system = self.MufiSystem[board_id][tofpet_id]
636 key = (tofpet_id%2)*1000 + tofpet_channel
637 tmp = self.boardMaps['MuFilter'][board][self.slots[tofpet_id]]
638 if self.options.debug or not tmp.find('not')<0: print('debug',tmp,system,key,board,tofpet_id,tofpet_id%2,tofpet_channel)
639 sipmChannel = 99
640 if not key in self.TofpetMap[system]:
641 print('key does not exist',key)
642 print(system, key,board,tofpet_id, self.TofpetMap[system])
643 else:
644 sipmChannel = self.TofpetMap[system][key]-1
645 nSiPMs = abs(self.offMap[tmp][1])
646 nSides = abs(self.offMap[tmp][2])
647 direction = int(self.offMap[tmp][1]/nSiPMs)
648 detID = self.offMap[tmp][0] + direction*(sipmChannel//(nSiPMs))
649 sipm_number = sipmChannel%(nSiPMs)
650 if tmp.find('Right')>0:
651 sipm_number += nSiPMs
652 if not detID in digiMuFilterStore:
653 digiMuFilterStore[detID] = ROOT.MuFilterHit(detID,nSiPMs,nSides)
654 test = digiMuFilterStore[detID].GetSignal(sipm_number)
655 digiMuFilterStore[detID].SetDigi(QDC,TDC,sipm_number)
656 digiMuFilterStore[detID].SetDaqID(sipm_number, board_id, tofpet_id, tofpet_channel)
657 if mask: digiMuFilterStore[detID].SetMasked(sipm_number)
658
659 if self.options.debug:
660 print('create mu hit: ',detID,tmp,system,tofpet_id,self.offMap[tmp],sipmChannel,nSiPMs,nSides,test)
661 print(' ',detID,sipm_number,QDC,TDC)
662 if test>0 or detID%1000>200 or sipm_number>15:
663 print('what goes wrong?',detID,sipm_number,system,key,board,tofpet_id,tofpet_channel,test)
664
665 else:
666 # scifi encoding
667 chan = self.channelchannel(tofpet_id,tofpet_channel,mat)
668 orientation = 1
669 if station[2]=="Y": orientation = 0
670 sipmLocal = (chan - mat*512)
671 sipmID = 1000000*int(station[1]) + 100000*orientation + 10000*mat + 1000*(sipmLocal//128) + chan%128
672 if not sipmID in digiSciFiStore:
673 digiSciFiStore[sipmID] = ROOT.sndScifiHit(sipmID)
674 digiSciFiStore[sipmID].SetDigi(QDC,TDC)
675 digiSciFiStore[sipmID].SetDaqID(0, board_id, tofpet_id, tofpet_channel)
676 if mask: digiSciFiStore[sipmID].setInvalid()
677 if self.options.debug:
678 print('create scifi hit: tdc = ',board,sipmID,QDC,TDC)
679 print('tofpet:', tofpet_id,tofpet_channel,mat,chan)
680 print(station[1],station[2],mat,chan,int(chan/128)%4,chan%128)
681
682# copy hits to detector branches
683 for sipmID in digiSciFiStore:
684 if self.digiSciFi.GetSize() == indexSciFi: self.digiSciFi.Expand(indexSciFi+100)
685 self.digiSciFi[indexSciFi]=digiSciFiStore[sipmID]
686 indexSciFi+=1
687 for detID in digiMuFilterStore:
688 if self.digiMuFilter.GetSize() == indexMuFilter: self.digiMuFilter.Expand(indexMuFilter+100)
689 self.digiMuFilter[indexMuFilter]=digiMuFilterStore[detID]
690 indexMuFilter+=1
691
692
693 def Finalize(self):
694 if self.options.FairTask_convRaw:
695 # overwrite cbmsim
696 F = self.outfile.GetRootFile()
697 T = F.Get("cbmsim")
698 T.SetName('rawConv')
699 T.Write()
700 self.fiN.Close()
701 self.outfile.Close()
702 else:
703 if self.options.debug:
704 print('number of events processed',self.sTree.GetEntries(),self.fiN.event.GetEntries())
705 self.sTree.Write()
706 self.fiN.Close()
707 self.fSink.Close()
708 print("File closed ",self.fiN.GetName())
ConvRawData.ConvRawDataPY.executeEvent
executeEvent(self, eventNumber)
Definition ConvRawData.py:417
ConvRawData.ConvRawDataPY.qdc_sat
qdc_sat(self, board_id, tofpet_id, channel, tac, v_fine)
Definition ConvRawData.py:316
ConvRawData.ConvRawDataPY.executeEvent0
executeEvent0(self, eventNumber)
Definition ConvRawData.py:560
ConvRawData.ConvRawDataPY.qdc_chi2
qdc_chi2(self, board_id, tofpet_id, channel, tac, TDC=0)
Definition ConvRawData.py:310
ConvRawData.ConvRawDataPY.executeEvent1
executeEvent1(self, eventNumber)
Definition ConvRawData.py:420
ConvRawData.ConvRawDataPY.comb_calibration
comb_calibration(self, board_id, tofpet_id, channel, tac, v_coarse, v_fine, t_coarse, t_fine, GQDC=1.0, TDC=0)
Definition ConvRawData.py:321
ConvRawData.ConvRawDataPY.debugMapping
debugMapping(self, board, tofpet_id, tofpet_channel)
Definition ConvRawData.py:352
ConvRawData.ConvRawDataPY.channel
channel(self, tofpet_id, tofpet_channel, position)
Definition ConvRawData.py:305
boardMappingParser.getBoardMapping
getBoardMapping(jsonString)
Definition boardMappingParser.py:5