--- /dev/null
+import datetime
+import glob
+from hldLib import *
+import numpy as np
+
+def tdcTimeData(w):
+ if not (w & 0x80000000):
+ return None
+
+ return {
+ "coarseTime": (w >> 0) & 0x7ff,
+ "edge": (w >> 11) & 1,
+ "fineTime": (w >> 12) & 0x3ff,
+ "channelNo": (w >> 22) & 0x7f
+ }
+
+thresh = -1
+counter = 0
+TDCCHANS = 5
+
+files = glob.glob("/u/mpenschuck/*.hld")
+files.sort()
+files.pop()
+print "# Read file ", files[-1]
+
+hits = np.zeros( (TDCCHANS,2,1024) ).astype(int)
+hitTimes = np.zeros( TDCCHANS ).astype("double")
+
+invChan = 0
+
+with open(files[-1], "rb") as hld:
+ for evt in eventIterator(hld):
+ for sEvt in evt[1]:
+ #print dumpEvt(evt)
+ for ssEvt in sEvt[1]:
+ if ssEvt[0]['subsubEvtId'] == 0xf3c0:
+ cts, remain = extractCTS(ssEvt[1])
+ sync, tdc = extractSync(remain)
+ #tdc, remain = extractTDC(remain)
+ assert(len(sync) in {1, 5, 11})
+
+ if 0:
+ print "CTS: ", dumpArray(cts)
+ print "Sync: ", dumpArray(sync)
+ print "TDC: ", dumpArray(tdc)
+
+ for w in tdc:
+ if w & 0x80000000:
+ tdcData = tdcTimeData(w)
+ if not tdcData: continue
+ #print tdcData["channelNo"]
+ if tdcData["channelNo"] < TDCCHANS:
+ hits[tdcData["channelNo"], tdcData["edge"], tdcData["fineTime"]] += 1
+ else:
+ invChan += 1
+
+ elif (w >> 29) == 1: # tdc header
+ #assert( (w >> 16) & 0xff == evt[0]['trgCode'] )
+ #print "hdr"
+ pass
+
+
+ counter += 1
+ if thresh > 0 and thresh <= counter:
+ break
+
+print "#Processed: %d events" % counter
+print "#Found %d valid and %d invalid hits: %s" % (np.sum(hits), invChan, ", ".join(["(chan: %d, fall: %d, rise: %d)" % (i, np.sum(hits[i,0,:]), np.sum(hits[i,1,:])) for i in range(TDCCHANS)]))
+
+calibInNs = np.cumsum(hits[:,1,:], axis=1).astype('double') / (np.sum(hits[:,1,:], axis=1)[:,np.newaxis] + 1e-100)
+calibInNs[:-1] += 0.5 * (calibInNs[1:]-calibInNs[:-1])
+calibInNs *= 5.0 / np.max(calibInNs, axis=1)[:,np.newaxis]
+np.save("calib.dat", calibInNs)
+
+np.savetxt("hist.txt", hits[:,1,:].T, fmt='%d')
+np.savetxt("calib.txt", calibInNs.T)
\ No newline at end of file
--- /dev/null
+import datetime
+import glob
+from hldLib import *
+import numpy as np
+
+def slopeWithoutAvg(X):
+ slope = X[1:]-X[:-1]
+ return slope - np.average(slope)
+
+def slope(X):
+ return X[1:] - X[:-1]
+
+def tdcTimeData(w):
+ if not (w & 0x80000000):
+ return None
+
+ return {
+ "coarseTime": (w >> 0) & 0x7ff,
+ "edge": (w >> 11) & 1,
+ "fineTime": (w >> 12) & 0x3ff,
+ "channelNo": (w >> 22) & 0x7f
+ }
+
+thresh = 150
+counter = 0
+TDCCHANS = 5
+
+files = glob.glob("/u/mpenschuck/*.hld")
+files.sort()
+
+epoch = 0
+epochTimeInNs = 0.0
+
+calibDataTDC = np.load("calib.dat.npy")
+
+firstEpoch = -1
+events = []
+
+with open(files[-1], "rb") as hld:
+ lastEvent = None
+ for evt in eventIterator(hld):
+ if not len(evt[1]): continue
+
+ event = {
+ 'trgCode': evt[1][0][0]['trgCode'],
+ 'sync': None,
+ 'hits': [[] for i in range(TDCCHANS)]
+ }
+
+ for sEvt in evt[1]:
+ for ssEvt in sEvt[1]:
+ if ssEvt[0]['subsubEvtId'] == 0xf3c0:
+ cts, remain = extractCTS(ssEvt[1])
+ sync, tdc = extractSync(remain)
+ assert(len(sync) in {1, 5, 11})
+
+ event['sync'] = sync
+
+ if lastEvent:
+ for w in tdc:
+ if w & 0x80000000:
+ tdcData = tdcTimeData(w)
+ if tdcData["edge"] == 1 and tdcData["channelNo"] < TDCCHANS and tdcData["fineTime"] != 0x3ff:
+ fineTimeInNs = calibDataTDC[tdcData["channelNo"], tdcData["fineTime"]]
+ assert( 0 <= fineTimeInNs <= 5.0 )
+ coarseTimeInNs = tdcData["coarseTime"] * 5.0
+ lastEvent['hits'][ tdcData["channelNo"] ].append( coarseTimeInNs - fineTimeInNs + epochTimeInNs )
+ #if tdcData["channelNo"] in {0,3}: print tdcData, fineTimeInNs, coarseTimeInNs - fineTimeInNs
+
+ elif (w >> 29) == 3: # epoch counter
+ epoch = w & 0x0fffffff
+ if firstEpoch < 0:
+ firstEpoch = epoch
+ epochTimeInNs = (epoch - firstEpoch) * 10240.0
+
+ elif (w >> 29) == 1: # tdc header
+ if (w >> 16) & 0xff != lastEvent['trgCode']:
+ break
+
+ elif (w >> 29) == 2: # debug
+ pass
+ else:
+ print "Unknown TDC word type: 0x%08x" % w
+
+ if lastEvent and lastEvent['sync'] != None and len(lastEvent['hits'][0]) == 1 and len(lastEvent['hits'][3]) == 1:
+ events.append(lastEvent)
+
+ lastEvent = event
+
+ counter += 1
+ if thresh > 0 and thresh <= counter:
+ break
+
+cbmEvtTimestamps = np.array( [1.0*evt['sync'][2] for evt in events] )
+cbmEvtTimestamps[1:] = np.cumsum( cbmEvtTimestamps[1:] - cbmEvtTimestamps[:-1] ) * 8.0
+cbmEvtTimestamps[0] = 0.0
+
+trbEvtTimestamps = np.array( [1.0*evt['sync'][1] for evt in events] )
+trbEvtTimestamps[1:] = np.cumsum( trbEvtTimestamps[1:] - trbEvtTimestamps[:-1] ) * 10.0
+trbEvtTimestamps[0] = 0.0
+
+factor = np.average( slope(trbEvtTimestamps) / slope(cbmEvtTimestamps) )
+
+tdcTrbRefTimestamp = np.array( [evt['hits'][0][0] for evt in events] ) * factor
+tdcCbmRefTimestamp = np.array( [evt['hits'][3][0] for evt in events] ) * factor
+timestampDelta = tdcCbmRefTimestamp - tdcTrbRefTimestamp
+
+offset=0
+if offset > 0:
+ cbmEvtTimestamps = cbmEvtTimestamps[offset:]
+ timestampDelta = timestampDelta[:-offset]
+elif offset < 0:
+ timestampDelta = timestampDelta[-offset:]
+ cbmEvtTimestamps = cbmEvtTimestamps[:offset]
+
+
+timestamps = cbmEvtTimestamps - timestampDelta + timestampDelta[0]
+
+interval = timestamps[1:] - timestamps[:-1]
+
+print "Std of interval", np.std(interval)
+
+
+
+plotData = np.vstack([ slopeWithoutAvg(cbmEvtTimestamps), slopeWithoutAvg(timestampDelta), slope(timestamps) ])
+np.savetxt("jitter.txt", plotData.T)
+
+pulserTS = np.array( [ p for evt in events for p in evt['hits'][1]] )
+pulserInterval = pulserTS[1:] - pulserTS[:-1]
+
+#np.savetxt("timing_trigger.txt", np.vstack([slope(tdcTrbRefTimestamp), slope(tdcCbmRefTimestamp)]).T)
+
+np.savetxt("timing_trigger.txt", np.vstack([slope(tdcTrbRefTimestamp), slope(tdcCbmRefTimestamp)]).T)
+np.savetxt("timing_trigger_delta.txt", np.vstack([timestampDelta]).T)
--- /dev/null
+set terminal pdf
+set output "hist.pdf"
+
+set xlabel "Fine-Time"
+set ylabel "Hits"
+set key left
+set xrange[0:511]
+
+plot \
+ "hist.txt" u 1 w lines title "Timing Trg (Trb)", \
+ "hist.txt" u 4 w lines title "Timing Trg (CBM)"
+# "hist.txt" u 2 w lines title "Pulser", \
+# "hist.txt" u 3 w lines title "DLM", \
+# , \
+# "hist.txt" u 5 w lines title "External"
+
+set output "calib.pdf"
+
+set ylabel "Time in ns"
+
+set yrange [-0.1 : 5.1]
+
+plot \
+ "calib.txt" u 1 w lines title "Timing Trg (Trb)", \
+ "calib.txt" u 2 w lines title "Pulser", \
+ "calib.txt" u 3 w lines title "DLM", \
+ "calib.txt" u 4 w lines title "Timing Trg (CBM)", \
+ "calib.txt" u 5 w lines title "External"
+
+
+set output "jitter.pdf"
+
+set ylabel "Offset in ns"
+set xlabel "Events"
+set xrange [*:*]
+set yrange [*:*]
+
+plot "jitter.txt" u 3 w linespoints title "CBM TS - (TDC(TimingTrg@CBM) - TDC(Tinimg@Trb))"
+
+set output "timing_trigger.pdf"
+set xlabel "Events"
+set ylabel "Time in ns"
+plot \
+ "timing_trigger.txt" u 1 w linespoints title "Timing Trigger (Trb)", \
+ "timing_trigger.txt" u 2 w linespoints title "Timing Trigger (CBM)"
+
+set output "timing_trigger_delta.pdf"
+plot \
+ "timing_trigger_delta.txt" u 1 w linespoints title "Timing Trigger (CBM-Trb)"
+
+
--- /dev/null
+import struct
+import datetime
+
+def eventIterator(f):
+ f.seek(0)
+ pos = 0
+
+ while(1):
+ # we might need to skip some padding
+ assert(f.tell() <= pos)
+ f.seek(pos)
+
+ evtHdrBuf = f.read(8*4)
+ if (len(evtHdrBuf) != 8*4):
+ break #eof
+
+ hdr = struct.unpack("<" + ("I" * 8), evtHdrBuf)
+ if (hdr[1] > 0xf0000): hdr = struct.unpack(">" + ("I" * 8), evtHdrBuf)
+
+ evtSize = hdr[0]
+
+ pos += evtSize
+ if (evtSize % 8):
+ pos += 8 - (evtSize % 8)
+
+ # decode hdr
+ hdrDec = {
+ 'size': evtSize,
+ 'triggerType': hdr[2] & 0xffff,
+ 'evtNumFile': hdr[3],
+ 'timestamp': datetime.datetime(
+ (hdr[4] >> 16) & 0xff, (hdr[4] >> 8) & 0xff, (hdr[4] >> 0) & 0xff,
+ (hdr[5] >> 16) & 0xff, (hdr[5] >> 8) & 0xff, (hdr[5] >> 0) & 0xff ),
+ 'runNum': hdr[6]
+ }
+
+ # load payload
+ subEvents = []
+ innerPos = 8*4
+ while(innerPos < evtSize):
+ subEvtHdrBuf = f.read(4 * 4)
+ endian = "<"
+ subEvtHdr = struct.unpack("<" + ("I" * 4), subEvtHdrBuf)
+ if (subEvtHdr[1] > 0xf0000):
+ subEvtHdr = struct.unpack(">" + ("I" * 4), subEvtHdrBuf)
+ endian = ">"
+
+ subEvtSize = subEvtHdr[0]
+ subEvtHdrDec = {
+ 'size': subEvtSize,
+ 'subEvtId': subEvtHdr[2] & 0xffff,
+ 'trgNum': (subEvtHdr[3] >> 8) & 0xffffff,
+ 'trgCode': subEvtHdr[3] & 0xff
+ }
+
+ subsubEvents = []
+ ssPos = 4*4
+ while(ssPos < subEvtSize):
+ sseHdrBuf = f.read(4)
+ sseHdr = struct.unpack(endian + "I", sseHdrBuf)
+ sseSize = ((sseHdr[0] >> 16) & 0xffff) * 4
+ sseHdrDec = {
+ 'size': sseSize,
+ 'subsubEvtId': sseHdr[0] & 0xffff
+ }
+
+ sseBuf = f.read(sseSize)
+ sseCont = struct.unpack(endian + ("I" * (sseSize/4)), sseBuf)
+ subsubEvents.append( (sseHdrDec, sseCont) )
+
+ ssPos += sseSize + 4
+
+ subEvents.append( (subEvtHdrDec, subsubEvents) )
+
+ innerPos += subEvtSize
+
+ yield (hdrDec, subEvents)
+
+def dumpEvt(evt):
+ res = str(evt[0]) + "\n"
+ for subEvt in evt[1]:
+ res += " " + dumpSubEvt(subEvt).replace("\n", "\n ") + "\n"
+ return res
+
+def dumpSubEvt(sEvt):
+ h = sEvt[0]
+ res = "subEvtId: 0x%04x, trgNum: % 9d, trgCode: % 4d, size: % 4d\n" % (h['subEvtId'], h['trgNum'], h['trgCode'], h['size'])
+ for ssEvt in sEvt[1]:
+ res += " " + dumpSubSubEvt(ssEvt).replace("\n", "\n ") + "\n"
+ return res
+
+def dumpSubSubEvt(ssEvt):
+ res = "ID: 0x%04x, Size: %d\n" % (ssEvt[0]['subsubEvtId'], ssEvt[0]['size'])
+ res += dumpArray(ssEvt[1])
+ return res
+
+def dumpArray(arr):
+ res = ""
+ for i in range(0, len(arr)):
+ if i != 0 and i % 8 == 0: res += "\n"
+ res += " [% 3d] 0x%08x" % (i+1, arr[i])
+
+ return res
+
+
+def extractCTS(ssEvtData):
+ hdr = ssEvtData[0]
+ length = 1 + \
+ ((hdr >> 16) & 0xf) * 2 + \
+ ((hdr >> 20) & 0x1f) * 2 + \
+ ((hdr >> 25) & 0x1) * 2 + \
+ ((hdr >> 26) & 0x1) * 3 + \
+ ((hdr >> 27) & 0x1) * 1
+
+ if (((hdr >> 28) & 0x3) == 0x1): length += 1
+ if (((hdr >> 28) & 0x3) == 0x2): length += 4
+
+ return (ssEvtData[:length], ssEvtData[length:])
+
+def extractSync(ssEvtData):
+ hdr = ssEvtData[0]
+ assert(hdr >> 28 == 0x1)
+ packType = (hdr >> 26) & 0x3
+ length = 1
+ if (packType == 1): length += 4
+ if (packType == 3): length += 10
+
+ return (ssEvtData[:length], ssEvtData[length:])
+
+
+def extractTDC(ssEvtData):
+ length = ssEvtData[0]
+ assert(length >= len(ssEvtData))
+ return (ssEvtData[:length+1], ssEvtData[length:])
+
jspc29 / trb3.git / commitdiff