examples/InteractiveStreamingExample/streamDetSource.py

0001 #!/usr/bin/python3.6
0002
0003 # Copyright 2020, Jefferson Science Associates, LLC.
0004 # Subject to the terms in the LICENSE file found in the top-level directory.
0005
0006 import argparse
0007 import numpy as np
0008 import matplotlib.pyplot as plt
0009 from scipy.stats import gumbel_r
0010 import pandas as pd
0011
0012 parser = argparse.ArgumentParser()
0013 parser.add_argument('-sr', '--sampleRate', metavar='SAMPLERATE', type=int, nargs=1, required=True,
0014                     help='sampling rate in MHz of SAMPA chips -> [5], [10], or [20]')
0015 parser.add_argument('-nc', '--numChans', metavar='NUMCHANS', type=int, nargs=1, required=True,
0016                     help='number of channels to simulate')
0017 parser.add_argument('-ne', '--numEvents', metavar='NUMEVENTS', type=int, nargs=1, required=True,
0018                     help='number of events to simulate')
0019 parser.add_argument('-spec', '--spectra', metavar='SPECTRA', type=str, nargs=1, required=True,
0020                     help='ADC spectra to simulate -> [gumbel] or [sampa]')
0021 parser.add_argument('-m', '--mode', metavar='MODE', type=str, nargs=1,
0022                     help='SAMPA DAQ mode to simulate -> [das] or [dsp]')
0023 args = parser.parse_args()
0024
0025 sampleRate = args.sampleRate[0]
0026 numChans = args.numChans[0]
0027 numEvents = args.numEvents[0]
0028 specMode = args.spectra[0]
0029 if args.mode : adcMode = args.mode[0]
0030
0031 def sampaHitFunc(sample, peak, startTime, decayTime, baseLine):
0032     adcSample = np.piecewise(sample, [sample < startTime, sample >= startTime], [lambda sample: baseLine,
0033                              lambda sample: (peak * np.power(((sample - startTime) / decayTime), 4) *
0034                                              np.exp((-4) * ((sample - startTime) / decayTime)) + baseLine)])
0035     return adcSample
0036
0037 class DataStream:
0038
0039     def __init__(self, sampleRate, numEvents):
0040         # sampling rate (s^-1), number of events to simulate
0041         self.sampleRate, self.numEvents = sampleRate, numEvents
0042         # number of samples in read-out window, up-to 1024 samples for SAMPA chips
0043         self.numSamples = 1024
0044         # read-out window length (s)
0045         self.windowWidth = (1.0 / self.sampleRate) * self.numSamples
0046         # length of run (s)
0047         self.runTime = self.numEvents * self.windowWidth
0048         # initial time sample
0049         self.initTime = 0.0
0050
0051     # method to simulate time samples
0052     def simTimeSamples(self, windowStart, windowEnd):
0053         self.timeSamples = np.linspace(windowStart, windowEnd, int(self.numSamples))
0054         self.initTime = self.timeSamples[-1]
0055         self.randTimeIndex = np.random.randint(int(len(self.timeSamples) * 0.25),
0056                                                high=int(len(self.timeSamples) * 0.75))
0057
0058     # method to simulate arbitrary adc samples via a gumbel distribution
0059     # def simAdcSignals(self, numSamples, randTimeIndex) :
0060     def gumbelAdcSignals(self):
0061         self.baseLine = np.random.rand(int(self.numSamples)) * 1.0e-2
0062         self.hitJudge = np.random.random()
0063         if (self.hitJudge < 0.5):  # simulate no hit
0064             self.adcSamples = self.baseLine
0065         elif (self.hitJudge >= 0.5):  # simulate hit
0066             self.gumbelMean = np.random.random()
0067             self.gumbelBeta = np.random.random()
0068             self.gumbelPpf = np.linspace(gumbel_r.ppf(0.001, loc=self.gumbelMean, scale=self.gumbelBeta),
0069                                          gumbel_r.ppf(0.999, loc=self.gumbelMean, scale=self.gumbelBeta), 100)
0070             self.gumbelPdf = gumbel_r.pdf(self.gumbelPpf, loc=self.gumbelMean, scale=self.gumbelBeta)
0071             self.adcSamples = np.insert(self.baseLine, self.randTimeIndex, self.gumbelPdf)[:int(self.numSamples)]
0072
0073     # method to simulate sampa adc samples via piecewise fourth order semi-gaussian
0074     def sampaAdcSignals(self):
0075         self.baseLine = np.random.randint(60, 81, size=self.numSamples)
0076         self.hitJudge = np.random.randint(0, 2, size=1)[0]
0077         if (self.hitJudge == 0):  # simulate no hit
0078             self.adcSamples = self.baseLine
0079         elif (self.hitJudge == 1):  # simulate hit
0080             self.peak = np.random.randint(5000, 50001, size=1)[0]
0081             self.baseLine = np.random.randint(60, 81, size=self.numSamples)
0082             self.startTime = np.random.randint(4, 7, size=1)[0]
0083             self.decayTime = np.random.randint(2, 5, size=1)[0]
0084             self.numHitSamples = np.random.randint(10, 16, size=1)[0]
0085             self.samples = np.arange(0, self.numHitSamples + 11)
0086             self.hitSamples = sampaHitFunc(self.samples, self.peak, self.startTime, self.decayTime, np.average(self.baseLine))
0087             self.adcSamples = np.insert(self.baseLine, self.randTimeIndex, self.hitSamples)[:int(self.numSamples)]
0088
0089     def __iter__(self):
0090         return self
0091
0092     def __next__(self):
0093         if (self.initTime == 0.0):
0094             # simulate time and adc samples
0095             self.simTimeSamples(self.initTime, self.windowWidth)
0096             if specMode == 'gumbel': self.gumbelAdcSignals()
0097             if specMode == 'sampa': self.sampaAdcSignals()
0098         elif (self.initTime > 0.0 and self.initTime < self.runTime):
0099             # simulate the time and adc signals
0100             self.simTimeSamples(self.initTime, self.initTime + self.windowWidth)
0101             if specMode == 'gumbel': self.gumbelAdcSignals()
0102             if specMode == 'sampa': self.sampaAdcSignals()
0103         elif (self.initTime >= self.runTime):
0104             raise StopIteration
0105
0106 # # data file for streaming analysis
0107 # datFile = open('run-%2.0d-mhz-%d-chan-%d-ev.dat' % (sampleRate, numChans, numEvents), 'w+')
0108 # for chan in range(1, numChans+1) :
0109 #     datFile.write('# channel = %d\n' % chan)
0110 #     dataObj = DataStream(sampleRate*1.0e+6, numEvents)
0111 #     for event in dataObj :
0112 #         np.savetxt(datFile, (dataObj.timeSamples, dataObj.adcSamples), fmt='%.9f')
0113 #    # plt.plot(dataObj.timeSamples, dataObj.adcSamples)
0114 #    # plt.xlabel('TDC Sample (arb. units)')
0115 #    # plt.ylabel('ADC Sample (arb. units)')
0116 #    # plt.title('Channel %d' % chan)
0117 #    # plt.show()
0118 # datFile.close()
0119
0120 # gumbel data file for event based analysis
0121 if specMode == 'gumbel' :
0122     datFile = open('run-%d-mhz-%d-chan-%d-ev.dat' % (sampleRate, numChans, numEvents), 'w+')
0123     for ievent in range(1, numEvents + 1):
0124         datFile.write('@ event = %d\n' % ievent)
0125         for chan in range(1, numChans + 1):
0126             datFile.write('# channel = %d\n' % chan)
0127             dataObj = DataStream(sampleRate * 1.0e+6, 1)
0128             for event in dataObj:
0129                 if len(dataObj.timeSamples) != len(dataObj.adcSamples): print("!!!Something terrible is amiss!!!")
0130                 np.savetxt(datFile, (np.add(dataObj.timeSamples, dataObj.windowWidth * (ievent - 1)), dataObj.adcSamples), fmt='%.9f')
0131     # plt.plot(dataObj.timeSamples, dataObj.adcSamples)
0132     # plt.xlabel('TDC Sample (arb. units)')
0133     # plt.ylabel('ADC Sample (arb. units)')
0134     # plt.title('Channel %d' % chan)
0135     # plt.show()
0136     datFile.close()
0137
0138 if specMode == 'sampa' :
0139     datFile = open('run-%d-mhz-%d-chan-%d-ev.dat' % (sampleRate, numChans, numEvents), 'w+')
0140     dfl = []; esl = []  # data frame list, event series list
0141     for ievent in range(1, numEvents + 1):
0142         csl = []  # channel series list
0143         for chan in range(1, numChans + 1):
0144             dataObj = DataStream(sampleRate * 1.0e+6, 1)
0145             for event in dataObj:
0146                 if len(dataObj.timeSamples) != len(dataObj.adcSamples): print("!!!Something terrible is amiss!!!")
0147                 csl.append(pd.Series(dataObj.adcSamples, name = 'chan_%d' % chan))
0148         esl.append(csl)
0149         dfl.append(pd.concat(esl[ievent-1], axis = 1))
0150     df = pd.concat(dfl, ignore_index=True)
0151     np.savetxt(datFile, df.values, fmt='%04d')
0152     # df.plot(y = 'chan_1', use_index = True, marker = 'o', c = 'tab:blue', ls = '')
0153     # plt.xlabel('TDC Sample (arb. units)')
0154     # plt.ylabel('ADC Sample (arb. units)')
0155     # plt.title('Channel %d' % chan)
0156     # plt.show()
0157     datFile.close()