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