replaced scipy interpolation by dedicated routine (by Artur)

libeos/file_reader.py

@@ -7,7 +7,6 @@ from typing import List
 
 import h5py
 import numpy as np
-import scipy as sp
 from orsopy import fileio
 from orsopy.fileio.model_language import SampleModel
 
@@ -21,6 +20,17 @@ try:
 except Exception:
     nb_helpers = None
 
+def get_current_per_pulse(t_pulses, t_currents, currents):
+    # add currents for early pulses and current time value after last pulse (j+1)
+    t_currents = np.hstack([[0], t_currents, [t_pulses[-1]+1]])
+    currents = np.hstack([[0], currents])
+    pulse_currents = np.zeros(t_pulses.shape[0], dtype=float)
+    j = 0
+    for i, ti in enumerate(t_pulses):
+        if ti>=t_currents[j+1]: 
+            j+=1
+        pulse_currents[i] = currents[j]
+    return pulse_currents
 
 class AmorData:
     """read meta-data and event streams from .hdf file(s), apply filters and conversions"""
@@ -236,8 +246,7 @@ class AmorData:
         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, fill_value=0)
-            self.monitorPerPulse = np.array(currentInterpolator(self.pulseTimeS) * 2*self.tau *1e-3, dtype=float)
+            self.monitorPerPulse = get_current_per_pulse(self.pulseTimeS, self.currentTime, self.current) * 2*self.tau * 1e-3
             # filter low-current pulses
             self.monitorPerPulse = np.where(self.monitorPerPulse > 2*self.tau *lowCurrentThreshold, self.monitorPerPulse, 0)
             # remove 'partially filled' pulses
@@ -248,16 +257,6 @@ class AmorData:
         else: 
             self.monitorPerPulse = 1./np.shape(pulseTimeS)[1]
 
-    #def define_monitor(self):
-    #    if self.monitorType == 'protonCharge':
-    #        chargeSum = np.sum(self.monitorPerPulse)
-    #        logging.warning(f'      proton charge = {chargeSum:9.3f} mC')
-    #        self.monitor = chargeSum
-    ##    elif self.monitorType == 'countingTime':
-    #        self.monitor = self.stopTime - self.seriesStartTime
-    #    else:
-    #        self.monitor = 1.
-
     def extract_walltime(self, norm):
         if nb_helpers:
             self.wallTime_e = nb_helpers.extract_walltime(self.tof_e, self.dataPacket_p, self.dataPacketTime_p)
@@ -266,8 +265,6 @@ class AmorData:
             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 -= np.int64(self.seriesStartTime)
         logging.debug(f'      wall time from {self.wallTime_e[0]/1e9:6.1f} s to {self.wallTime_e[-1]/1e9:6.1f} s')
 

libeos/reduction.py

@@ -184,15 +184,16 @@ class AmorReduction:
             stop = wallTime_e[-1]
         # make overwriting log lines possible by removing newline at the end
         #logging.StreamHandler.terminator = "\r"
+        logging.warning(f'    time slizing')
         for ti, time in enumerate(np.arange(start, stop, interval)):
-            logging.warning(f'    time slize {ti:4d}')
+            logging.info(f'      slize {ti:4d}  t = {time:5.0f}')
 
             filter_e = np.where((time<wallTime_e) & (wallTime_e<time+interval), True, False)
             lamda_e = self.file_reader.lamda_e[filter_e]
             detZ_e = self.file_reader.detZ_e[filter_e]
             filter_m = np.where((time<pulseTimeS) & (pulseTimeS<time+interval), True, False)
             self.monitor = np.sum(self.file_reader.monitorPerPulse[filter_m])
-            logging.warning(f'      monitor = {self.monitor:7.2f}')
+            logging.info(f'        monitor = {self.monitor:7.2f}')
 
             qz_lz, qx_lz, ref_lz, err_lz, res_lz, lamda_lz, theta_lz, int_lz, mask_lz = self.project_on_lz(
                     self.file_reader, self.norm_lz, self.normAngle, lamda_e, detZ_e)
@@ -226,7 +227,7 @@ class AmorReduction:
             self.datasetsRqz.append(orso_data)
         # reset normal logging behavior
         #logging.StreamHandler.terminator = "\n"
-        logging.warning(f'    time slizing, done')
+        logging.info(f'      done')
 
     def save_Rqz(self):
         fname = os.path.join(self.reader_config.dataPath, f'{self.output_config.outputName}.Rqz.ort')
@@ -342,7 +343,6 @@ class AmorReduction:
             lamda_e          = fromHDF.lamda_e
             detZ_e           = fromHDF.detZ_e
             self.normMonitor = np.sum(fromHDF.monitorPerPulse)
-            print(self.normMonitor)
             self.norm_lz, bins_l, bins_z = np.histogram2d(lamda_e, detZ_e, bins = (self.grid.lamda(), self.grid.z()))
             self.norm_lz = np.where(self.norm_lz>2, self.norm_lz, np.nan)
             # correct for the SM reflectivity