From ee077b19f7d1b7d0b40762fa9caa8a0eb06e15ec Mon Sep 17 00:00:00 2001 From: jochenstahn Date: Wed, 25 Sep 2024 10:37:43 +0200 Subject: [PATCH] cleaned and filtering for low proton current --- libeos/file_reader.py | 75 +++++++++++++++++++++---------------------- libeos/nb_helpers.py | 3 +- 2 files changed, 39 insertions(+), 39 deletions(-) diff --git a/libeos/file_reader.py b/libeos/file_reader.py index 355f5bd..11e1424 100644 --- a/libeos/file_reader.py +++ b/libeos/file_reader.py @@ -178,14 +178,16 @@ class AmorData: self.read_event_stream() totalNumber = np.shape(self.tof_e)[0] - self.sort_events_by_pulse() + self.sort_pulses() + + self.associate_pulse_with_current() self.define_monitor() - # sort the events into the related pulses - self.extract_walltime(norm) + self.monitor_threshold() + self.filter_strange_times() self.merge_frames() @@ -200,7 +202,7 @@ class AmorData: logging.info(f' number of events: total = {totalNumber:7d}, filtered = {np.shape(self.lamda_e)[0]:7d}') - def sort_events_by_pulse(self): + def sort_pulses(self): chopperPeriod = int(2*self.tau*1e9) pulseTime = np.sort(self.dataPacketTime_p) pulseTime = pulseTime[np.abs(pulseTime[:]-np.roll(pulseTime, 1)[:])>5] @@ -221,33 +223,22 @@ class AmorData: self.pulseTimeS = np.append(self.pulseTimeS, nxt) nxt += chopperPeriod self.pulseTimeS = np.append(self.pulseTimeS, tt) - # remove 'partially filled' pulses - self.pulseTimeS = self.pulseTimeS[1:-1] def associate_pulse_with_current(self): if self.monitorType == 'protonCharge': + lowCurrentThreshold = 0.05 # mA self.currentTime -= self.seriesStartTime - currentInterpolator = sp.interpolate.interp1d(self.currentTime, self.current, kind='previous', bounds_error=False) - charge = np.array(currentInterpolator(self.pulseTimeS) * 2*self.tau *1e-3, dtype=float) + currentInterpolator = sp.interpolate.interp1d(self.currentTime, self.current, kind='previous', bounds_error=False, fill_value=0) + self.charge = np.array(currentInterpolator(self.pulseTimeS) * 2*self.tau *1e-3, dtype=float) # filter low-current pulses - charge = np.where(charge > 2*self.tau * 1e-1, charge, 0) - #self.event - # TODO: activate the following filter AND remove the respective events : - # remove events (wallTime, tof and pixelID) from stream for pulses with too low monitor signal - # probably using np.isin() - + self.charge = np.where(self.charge > 2*self.tau *lowCurrentThreshold, self.charge, 0) + # remove 'partially filled' pulses + self.charge[0] = 0 + self.charge[-1] = 0 def define_monitor(self): if self.monitorType == 'protonCharge': - # associate each pulse with a proton current - self.currentTime -= self.seriesStartTime - currentInterpolator = sp.interpolate.interp1d(self.currentTime, self.current, kind='previous', bounds_error=False, fill_value=0) - charge = np.array(currentInterpolator(self.pulseTimeS) * 2*self.tau *1e-3, dtype=float) - # TODO: activate the following filter AND remove the respective events : - # remove events (wallTime, tof and pixelID) from stream for pulses with too low monitor signal - # probably using np.isin() - #charge = np.where(charge > 2*self.tau * 1e-1, charge, 0) - chargeSum = np.sum(charge) + chargeSum = np.sum(self.charge) logging.warning(f' proton charge = {chargeSum:9.3f} mC') self.monitor = chargeSum elif self.monitorType == 'countingTime': @@ -255,6 +246,29 @@ class AmorData: else: self.monitor = 1. + def extract_walltime(self, norm): + #self.dataPacketTime_p = np.array(self.dataPacketTime_p, dtype=float) / 1e9 + if nb_helpers: + self.wallTime_e = nb_helpers.extract_walltime(self.tof_e, self.dataPacket_p, self.dataPacketTime_p) + else: + self.wallTime_e = np.empty(np.shape(self.tof_e)[0], dtype=int) + for i in range(len(self.dataPacket_p)-1): + self.wallTime_e[self.dataPacket_p[i]:self.dataPacket_p[i+1]] = self.dataPacketTime_p[i] + self.wallTime_e[self.dataPacket_p[-1]:] = self.dataPacketTime_p[-1] + #if not self.startTime and not norm: + # self.startTime = self.wallTime_e[0] + self.wallTime_e -= self.seriesStartTime + logging.debug(f' wall time from {self.wallTime_e[0]/1e9} to {self.wallTime_e[-1]/1e9}') + + def monitor_threshold(self): + goodTimeS = self.pulseTimeS[self.charge!=0] + filter_e = np.where(np.isin(self.wallTime_e, goodTimeS), True, False) + self.tof_e = self.tof_e[filter_e] + self.pixelID_e = self.pixelID_e[filter_e] + self.wallTime_e = self.wallTime_e[filter_e] + logging.warning(f' rejected {np.shape(self.charge)[0]-np.shape(goodTimeS)[0]} pulses due to low beam current') + logging.warning(f' rejected {np.shape(filter_e)[0]-np.shape(self.tof_e)[0]} events due to low beam current') + def filter_qz_range(self, norm): if self.config.qzRange[1]<0.3 and not norm: self.mask_e = np.logical_and(self.mask_e, @@ -331,21 +345,6 @@ class AmorData: if np.shape(filter_e)[0]-np.shape(self.tof_e)[0]>0.5: logging.warning(f'# strange times: {np.shape(filter_e)[0]-np.shape(self.tof_e)[0]}') - def extract_walltime(self, norm): - self.dataPacketTime_p = np.array(self.dataPacketTime_p, dtype=float) / 1e9 - if nb_helpers: - self.wallTime_e = nb_helpers.extract_walltime(self.tof_e, self.dataPacket_p, self.dataPacketTime_p) - else: - totalNumber = np.shape(self.tof_e)[0] - self.wallTime_e = np.empty(totalNumber) - for i in range(len(self.dataPacket_p)-1): - self.wallTime_e[self.dataPacket_p[i]:self.dataPacket_p[i+1]] = self.dataPacketTime_p[i] - self.wallTime_e[self.dataPacket_p[-1]:] = self.dataPacketTime_p[-1] - if not self.startTime and not norm: - self.startTime = self.wallTime_e[0] - self.wallTime_e -= self.startTime - logging.debug(f' wall time from {self.wallTime_e[0]} to {self.wallTime_e[-1]}') - def read_event_stream(self): self.tof_e = np.array(self.hdf['/entry1/Amor/detector/data/event_time_offset'][:])/1.e9 self.pixelID_e = np.array(self.hdf['/entry1/Amor/detector/data/event_id'][:], dtype=int) diff --git a/libeos/nb_helpers.py b/libeos/nb_helpers.py index 52efc95..c72d3b5 100644 --- a/libeos/nb_helpers.py +++ b/libeos/nb_helpers.py @@ -16,7 +16,8 @@ def merge_frames(tof_e, tofCut, tau, total_offset): def extract_walltime(tof_e, dataPacket_p, dataPacketTime_p): # assigning every event the wall time of the event packet (absolute time of pulse ?start?) totalNumber = np.shape(tof_e)[0] - wallTime_e = np.empty(totalNumber, dtype=np.float64) + #wallTime_e = np.empty(totalNumber, dtype=np.float64) + wallTime_e = np.empty(totalNumber, dtype=int) for i in nb.prange(len(dataPacket_p)-1): for j in range(dataPacket_p[i], dataPacket_p[i+1]): wallTime_e[j] = dataPacketTime_p[i]