Functions
	qdc_calibration (board_id, tofpet_id, channel, tac, v_coarse, v_fine, tf)

	qdc_chi2 (board_id, tofpet_id, channel, tac, TDC=0)

	qdc_sat (board_id, tofpet_id, channel, tac, v_fine)

	time_calibration (board_id, tofpet_id, channel, tac, t_coarse, t_fine, TDC=0)

	comb_calibration (board_id, tofpet_id, channel, tac, v_coarse, v_fine, t_coarse, t_fine, GQDC=1.0, TDC=0)

	calibrationReport ()

	channel (tofpet_id, tofpet_channel, position)

	run (nEvent)

	getMapEvent2Time ()

	asynInfo (run)

	dumpMap ()

Variables
dict	h = {}

int	chi2Max = 2000.

float	saturationLimit = 0.95

int	mufi_hitThreshold = 2

	parser = ArgumentParser()

	dest

	help

	type

	int

	required

	False

	default

bool	withGeoFile = False

	options = parser.parse_args()

	runNr = str(options.runNumber).zfill(6)

str	path = options.path+'run_'+ runNr+'/'

str	part = ""

str	inFile = 'data_'+part+'.root'

str	outFile = "sndsw_raw_"+runNr+'-'+part+'.root'

bool	local = False

	fqdc_cal = open(path+'/qdc_cal.csv')

	Lqdc = fqdc_cal.readlines()

	ftdc_cal = open(path+'/tdc_cal.csv')

	Ltdc = ftdc_cal.readlines()

	server = os.environ['EOSSHIP']

	status

	L = Lqdc

dict	SiPMmap = {}

dict	TofpetMap = {}

str	infile = "/eos/experiment/sndlhc/testbeam/MuFilter/SiPM_mappings/"+system+"_SiPM_mapping.csv"

	reader = csv.DictReader(f, delimiter=',')

	tmp = L[l].replace('\n','').split(',')

	row = x.split(',')

dict	qdc_cal = {}

	board_id = int(tmp[0])

	fe_id = int(tmp[1])

	channel = int(tmp[2])

	tac = int(tmp[3])

dict	X = qdc_cal[board_id][fe_id][channel][tac]

	tdc = int(tmp[4])

dict	stations = {}

dict	boardMaps = {}

str	board = 'board_'+str(stations[station][mat])

dict	slots = {0:'A',1:'A',2:'B',3:'B',4:'C',5:'C',6:'D',7:'D'}

dict	MufiSystem = {}

int	s = 0

dict	offMap = {}

str	dataName = 'data'+part+'.root'

	f0 = ROOT.TFile.Open(X+path+dataName)

	nEvent

dict	boards = {}

	name = b.GetName()

	fSink = ROOT.FairRootFileSink(outFile)

	sTree = ROOT.TTree('rawConv','raw data converted')

	header = ROOT.FairEventHeader()

	eventSND = sTree.Branch("EventHeader",header,32000,-1)

	digiSciFi = ROOT.TClonesArray("sndScifiHit")

	digiSciFiBranch = sTree.Branch("Digi_ScifiHits",digiSciFi,32000,1)

	digiMuFilter = ROOT.TClonesArray("MuFilterHit")

	digiMuFilterHitBranch = sTree.Branch("Digi_MuFilterHits",digiMuFilter,32000,1)

	clusScifi = ROOT.TClonesArray("sndCluster")

	clusScifiBranch = sTree.Branch("Cluster_Scifi",clusScifi,32000,1)

	B = ROOT.TList()

dict	theMap = {}

str	test = 'map'+str(options.runNumber)+'.pkl'

	fin = open('map'+str(options.runNumber)+'.pkl','rb')

	asyn = asynInfo(options.runNumber)

Function Documentation

◆ asynInfo()

convertRawData_muTestbeam.asynInfo ( run )

Definition at line 541 of file convertRawData_muTestbeam.py.

def asynInfo(run):
    server = os.environ['EOSSHIP']
    with client.File() as f:
      location = "/eos/experiment/sndlhc/raw_data/commissioning/TB_H8_october/run_logs/log_"+str(run)+".txt"
      f.open(server+location)
      status, L = f.read()
      A = L.decode().split('\n')
      f.close()
    info = {}
    for x in A:
      X = x.replace('\n','')
      if X=='':continue
      if X.find('2021')==0:
          tag = X
          info[tag]={}
      else:
          if X.find('disc')>0: info[tag]['summary'] = eval(X)
          else: info[tag]['detail'] = eval(X)
# correction procedure
    A = {41:{},43:{},59:{},60:{}}
    C = {}
    prev = {41:None,43:None,59:None,60:None}
    errorCounter = 0
    T = list(info.keys())
    T.sort()
    for tag in T:
          test = info[tag]['summary']
          lastT = test['last_timestamp']
          if not 'desync_info' in test: continue
          C[tag] = {41:{},43:{},59:{},60:{}}
          for L in test['desync_info']:
             board = L[0]
             t = L[1][0]
             c = L[1][1]
             C[tag][board][c]=t
    # check that desync is constant over the complete run
          for c in C[tag][41]:
              tref = C[tag][41][c]
              break
          for board in C[tag]:
              if not c in C[tag][board]:
                  #print('different trigger counts, problematic',tag,C[tag])
                  errorCounter+=1
              else:
                  if not prev[board] == None:
                         if not C[tag][board][c]-tref == prev[board]:
                              print('desync changed, problematic',tag,prev,C[tag])
                  prev[board] = C[tag][board][c]-tref
    
    T = list(C.keys())
    for b in C[T[0]]:
        for c in C[T[0]][b]:
            A[b]= C[T[0]][b][c]
    return A
 

◆ calibrationReport()

convertRawData_muTestbeam.calibrationReport ( )

Definition at line 182 of file convertRawData_muTestbeam.py.

def calibrationReport():
    ut.bookHist(h,'chi2','chi2',1000,0.,10000)
    report = {}
    TDC = 0
    for b in qdc_cal:
         for t in  qdc_cal[b]:
              for c in qdc_cal[b][t]:
                   for tac in qdc_cal[b][t][c]:
                             par    = qdc_cal[b][t][c][tac]
                             if 'chi2Ndof' in par:                             chi2 = par['chi2Ndof']
                             else: chi2=-1
                             parT = qdc_cal[b][t][c][tac][TDC]
                             if 'chi2Ndof' in parT:   chi2T = parT['chi2Ndof']
                             else: chi2T=-1
                             key = tac +10*c + t*10*100 + b*10*100*100
                             if not key in report: report[key]=[chi2,chi2T]
    for key in report:
         rc=h['chi2'].Fill(report[key][0])
         rc=h['chi2'].Fill(report[key][1])
    h['chi2'].Draw()
    return report
 
# station mapping for SciFi

◆ channel()

convertRawData_muTestbeam.channel	(	tofpet_id,
		tofpet_channel,
		position
	)

Definition at line 261 of file convertRawData_muTestbeam.py.

def channel(tofpet_id,tofpet_channel,position):
    return (64 * tofpet_id + 63 - tofpet_channel + 512*position) # 512 channels per mat, 1536 channels per plane
                                                                                                                                     # one channel covers all 6 layers of fibres
 
# reading hits and converting to event information

◆ comb_calibration()

convertRawData_muTestbeam.comb_calibration	(	board_id,
		tofpet_id,
		channel,
		tac,
		v_coarse,
		v_fine,
		t_coarse,
		t_fine,
		GQDC = `1.0`,
		TDC = `0`
	)

Definition at line 158 of file convertRawData_muTestbeam.py.

def comb_calibration(board_id,tofpet_id,channel,tac,v_coarse,v_fine,t_coarse,t_fine,GQDC = 1.0, TDC=0): # max gain QDC = 3.6
    par  = qdc_cal[board_id][tofpet_id][channel][tac]
    parT = par[TDC]
    x    = t_fine
    ftdc = (-parT['b']-ROOT.TMath.Sqrt(parT['b']**2-4*parT['a']*(parT['c']-x)))/(2*parT['a']) #   Ettore 28/01/2022 +parT['d']
    timestamp = t_coarse + ftdc
    tf = timestamp - t_coarse
    x = v_coarse - tf
    fqdc = -par['c']*ROOT.TMath.Log(1+ROOT.TMath.Exp( par['a']*(x-par['e'])**2-par['b']*(x-par['e']) )) + par['d']
    value = (v_fine-fqdc)/GQDC
    return timestamp,value,max(par['chi2Ndof'],parT['chi2Ndof']),v_fine/par['d']
 

◆ dumpMap()

convertRawData_muTestbeam.dumpMap ( )

Definition at line 597 of file convertRawData_muTestbeam.py.

def dumpMap():
 fout = open('map'+str(options.runNumber)+'.pkl','wb')
 pickle.dump(theMap,fout)
 

◆ getMapEvent2Time()

convertRawData_muTestbeam.getMapEvent2Time ( )

Definition at line 534 of file convertRawData_muTestbeam.py.

def getMapEvent2Time():
     k = 0
     for event in f0.event:  
          theMap[event.timestamp]=k
          k+=1
          if k%10000000==0: print('key time map: now at ',k)
 

◆ qdc_calibration()

convertRawData_muTestbeam.qdc_calibration	(	board_id,
		tofpet_id,
		channel,
		tac,
		v_coarse,
		v_fine,
		tf
	)

Definition at line 134 of file convertRawData_muTestbeam.py.

def qdc_calibration(board_id,tofpet_id,channel,tac,v_coarse,v_fine,tf):
    GQDC = 1.0 # or 3.6
    par = qdc_cal[board_id][tofpet_id][channel][tac]
    x = v_coarse - tf
    fqdc = -par['c']*ROOT.TMath.Log(1+ROOT.TMath.Exp( par['a']*(x-par['e'])**2-par['b']*(x-par['e']) )) + par['d']
    value = (v_fine-fqdc)/GQDC
    return value
 

◆ qdc_chi2()

convertRawData_muTestbeam.qdc_chi2	(	board_id,
		tofpet_id,
		channel,
		tac,
		TDC = `0`
	)

Definition at line 142 of file convertRawData_muTestbeam.py.

def qdc_chi2(board_id,tofpet_id,channel,tac,TDC=0):
    par    = qdc_cal[board_id][tofpet_id][channel][tac]
    parT = qdc_cal[board_id][tofpet_id][channel][tac][TDC]
    return max(par['chi2Ndof'],parT['chi2Ndof'])
 

◆ qdc_sat()

convertRawData_muTestbeam.qdc_sat	(	board_id,
		tofpet_id,
		channel,
		tac,
		v_fine
	)

Definition at line 147 of file convertRawData_muTestbeam.py.

def qdc_sat(board_id,tofpet_id,channel,tac,v_fine):
    par = qdc_cal[board_id][tofpet_id][channel][tac]
    return v_fine/par['d']
 

◆ run()

convertRawData_muTestbeam.run ( nEvent )

Definition at line 305 of file convertRawData_muTestbeam.py.

def run(nEvent):
 event = f0.event
 for eventNumber in range(options.nStart,nEvent):
   ncreated = 0
   rc = event.GetEvent(eventNumber)
   T = event.timestamp
   if eventNumber%options.heartBeat==0:
      tt = time.ctime()
      print('run ',options.runNumber, ' event',eventNumber," ",tt)
   # find entry numbers for all boards
   otherEvents={}
   maxHits = 0
   found = False
   for board in boards:
         board_id = int(board.split('_')[1])
         bt = boards[board]
         rc  = bt.GetEvent(eventNumber)
         if bt.n_hits < 1: continue
         if bt.n_hits > maxHits: 
             maxHits = bt.n_hits
             maxBoard = board_id
         t1 = asyn[board_id]
         for boardx in boards:
             boardx_id = int(boardx.split('_')[1])
             t2 = asyn[boardx_id]
             dt = t2 - t1
             newT = T + dt
             if newT<0 or dt==0: continue
             if newT in theMap:
                if options.debug: print("found",eventNumber)
                key = theMap[newT]
                if key < 0 : continue    # event time had already been used
                found = True
                if not board_id in otherEvents: 
                   otherEvents[board_id]={}                   
                otherEvents[board_id][boardx_id]=[newT,theMap[newT]]
# problem if more than one board has links to another event.
   if found and options.debug:     print("+++",otherEvents)
 
   if len(otherEvents)>1:
           print("Houston we have a problem, more than one board match with other events",eventNumber,otherEvents)
           for b in otherEvents:
               rc = boards["board_"+str(b)].GetEvent(eventNumber)
               txt = "mult "+str( boards["board_"+str(b)].n_hits )
               for bx in otherEvents[b]:
                   rc = boards["board_"+str(bx)].GetEvent(otherEvents[b][bx][1])
                   txt+= " " +str( boards["board_"+str(bx)].n_hits )
               print(txt)
           continue
 
   if len(otherEvents)==0:
           if maxHits>mufi_hitThreshold: 
                otherEvents[maxBoard]={}
                otherEvents[maxBoard][maxBoard] = [T,eventNumber]
           else: continue
   else:
           thisBoard = list(otherEvents.keys())[0]
           otherEvents[thisBoard][thisBoard] = [T,eventNumber]
   listOfBoards = {}
   for x in otherEvents:
      for b in otherEvents[x]:
         tmp = otherEvents[x][b]
         listOfBoards[b]=tmp[1]
         theMap[tmp[0]] = -1    # use event only once
 
   if options.debug: print("---",listOfBoards)
 
   header.SetEventTime(event.timestamp)
   header.SetRunId(options.runNumber)
   if options.debug: print('event:',eventNumber,event.timestamp)
   indexSciFi=0
   digiSciFi.Delete()
   digiSciFiStore = {}
   indexMuFilter=0
   digiMuFilter.Delete()
   digiMuFilterStore = {}
# special part for synchronizing the MuFilter testbeam data
   for board_id in listOfBoards:
      board = "board_"+str(board_id)
      scifi = True
      if board in boardMaps['Scifi']:
          station,mat = boardMaps['Scifi'][board]
      elif board in boardMaps['MuFilter']:
          scifi = False
      else: 
           print(board+' not known. Serious error, stop')
           1/0
      bt = boards[board]
      rc  = bt.GetEvent(listOfBoards[board_id])
      for n in range(bt.n_hits):
             mask = False
             if options.debug:
                  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])
             tofpet_id = bt.tofpet_id[n]
             tofpet_channel = bt.tofpet_channel[n]
             tac = bt.tac[n]
             TDC,QDC,Chi2ndof,satur = 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])
             if Chi2ndof > chi2Max:
                       if QDC>1E20:    QDC = 997.   # checking for inf
                       if QDC != QDC:  QDC = 998.   # checking for nan
                       if QDC>0: QDC = -QDC
                       mask = True
             elif satur  > saturationLimit or QDC>1E20 or QDC != QDC:
                       if QDC>1E20:    QDC = 987.   # checking for inf
                       if options.debug: 
                              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)
                       if QDC != QDC:  QDC = 988.   # checking for nan
                       if options.debug: 
                               print('nan',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],TDC,bt.t_coarse[n],eventNumber,Chi2ndof)
                       A = int(min( QDC,1000))
                       B = min(satur,999)/1000.
                       QDC = A+B
                       mask = True
             elif Chi2ndof > chi2Max:
                       if QDC>0: QDC = -QDC
                       mask = True
             if options.debug:
                  print('calibrated: tdc = ',TDC,'  qdc = ',QDC)  # TDC clock cycle = 160 MHz or 6.25ns
 
             if not scifi:
# mufi encoding
                slot = bt.tofpet_id[n]
                tmp = boardMaps['MuFilter'][board][slots[slot]]
                system = tmp.split('_')[0]
                key = (slot%2)*1000 + bt.tofpet_channel[n]
                if options.debug: print(system,key,board,bt.tofpet_id[n],bt.tofpet_id[n]%2,bt.tofpet_channel[n])
                sipmChannel = 99
                if not key in TofpetMap[system]:
                        print('key does not exist',key)
                        print(system, key, TofpetMap[system])
                else:
                       sipmChannel = TofpetMap[system][key]-1
                nSiPMs = abs(offMap[tmp][1])
                nSides =   abs(offMap[tmp][2])
                direction            = int(offMap[tmp][1]/nSiPMs)
                detID                   = offMap[tmp][0] + direction*(sipmChannel//(nSiPMs))
                sipm_number = sipmChannel%(nSiPMs)
                if tmp.find('Right')>0:
                       sipm_number += nSiPMs
                if not detID in digiMuFilterStore:
                        digiMuFilterStore[detID] =  ROOT.MuFilterHit(detID,nSiPMs,nSides)
                test = digiMuFilterStore[detID].GetSignal(sipm_number)
                digiMuFilterStore[detID].SetDigi(QDC,TDC,sipm_number)
                if mask: digiMuFilterStore[detID].SetMasked(sipm_number)
 
                if options.debug:
                    ncreated +=1
                    print('create mu hit: ',detID,tmp,system,slot,offMap[tmp],sipmChannel,nSiPMs,nSides,test,slot)
                    print('                ',detID,sipm_number,QDC,TDC)
                if test>0 or detID%1000>200 or sipm_number>15:
                    print('what goes wrong?',detID,sipm_number,system,key,board,bt.tofpet_id[n],bt.tofpet_channel[n],test)
             else:
# scifi encoding
                chan = channel(bt.tofpet_id[n],bt.tofpet_channel[n],mat)
                orientation = 1
                if station[2]=="Y": orientation = 0
                sipmLocal = (chan - mat*512)
                sipmID = 1000000*int(station[1]) + 100000*orientation + 10000*mat + 1000*(sipmLocal//128) + chan%128
                if not sipmID in digiSciFiStore: digiSciFiStore[sipmID] =  ROOT.sndScifiHit(sipmID)
                digiSciFiStore[sipmID].SetDigi(QDC,TDC)
                if mask: digiSciFiStore[sipmID].setInvalid()
                if options.debug:
                    print('create scifi hit: tdc = ',board,sipmID,QDC,TDC)
                    print('tofpet:', bt.tofpet_id[n],bt.tofpet_channel[n],mat,chan)
                    print(station[1],station[2],mat,chan,int(chan/128)%4,chan%128)
 
   for sipmID in digiSciFiStore:
             if digiSciFi.GetSize() == indexSciFi: digiSciFi.Expand(indexSciFi+100)
             digiSciFi[indexSciFi]=digiSciFiStore[sipmID]
             indexSciFi+=1
   for detID in digiMuFilterStore:
             if digiMuFilter.GetSize() == indexMuFilter: digiMuFilter.Expand(indexMuFilter+100)
             digiMuFilter[indexMuFilter]=digiMuFilterStore[detID]
             indexMuFilter+=1
# make simple clustering for scifi, only works with geometry file. Don't think it is a good idea at the moment
   if withGeoFile:
    index = 0
    hitDict = {}
    for k in range(digiSciFi.GetEntries()):
        d = digiSciFi[k]
        if not d.isValid(): continue 
        hitDict[d.GetDetectorID()] = k
        hitList = list(hitDict.keys())
        if len(hitList)>0:
              hitList.sort()
              tmp = [ hitList[0] ]
              cprev = hitList[0]
              ncl = 0
              last = len(hitList)-1
              hitlist = ROOT.std.vector("sndScifiHit*")()
              for i in range(len(hitList)):
                   if i==0 and len(hitList)>1: continue
                   c=hitList[i]
                   if (c-cprev)==1: 
                        tmp.append(c)
                   if (c-cprev)!=1 or c==hitList[last]:
                        first = tmp[0]
                        N = len(tmp)
                        hitlist.clear()
                        for aHit in tmp: 
                                hitlist.push_back( digiSciFi[hitDict[aHit]])
                                print(aHit,hitDict[aHit],digiSciFi[hitDict[aHit]].GetDetectorID())
                        print(hitlist.size())
                        aCluster = ROOT.sndCluster(first,N,hitlist,scifiDet)
                        print("cluster created")
                        if  clusScifi.GetSize() == index: clusScifi.Expand(index+10)
                        clusScifi[index]=aCluster
                        index+=1
                        if c!=hitList[last]:
                             ncl+=1
                             tmp = [c]
                   cprev = c
# fill TTree
   sTree.Fill()
   if options.debug: 
        print('====> number of created hits',ncreated)
 if options.debug:
     print('number of events processed',sTree.GetEntries(),f0.event.GetEntries())
 sTree.Write()
 
# https://gitlab.cern.ch/snd-scifi/software/-/wikis/Raw-data-format 
# tac: 0-3, identifies the Time-To-Analogue converter used for this hit (each channel has four of them and they require separate calibration).
# t_coarse: Coarse timestamp of the hit, running on a 4 times the LHC clock
# t_fine: 0-1023, fine timestamp of the hit. It contains the raw value of the charge digitized by the TDC and requires calibration.
# v_coarse: 0-1023, QDC mode: it represents the number of clock cycles the charge integration lasted.
# v_fine = 0-1023, QDC mode: represents the charge measured. Requires calibration.
 

◆ time_calibration()

convertRawData_muTestbeam.time_calibration	(	board_id,
		tofpet_id,
		channel,
		tac,
		t_coarse,
		t_fine,
		TDC = `0`
	)

Definition at line 151 of file convertRawData_muTestbeam.py.

def time_calibration(board_id,tofpet_id,channel,tac,t_coarse,t_fine,TDC=0):  #??
    par = qdc_cal[board_id][tofpet_id][channel][tac][TDC]
    x = t_fine
    ftdc = (-par['b']-ROOT.TMath.Sqrt(par['b']**2-4*par['a']*(par['c']-x)))/(2*par['a']) #   Ettore 28/01/2022 +par['d']
    timestamp = t_coarse+ftdc
    return timestamp
 