Add IPM data wrapper

.gitignore

@@ -2,6 +2,7 @@
 # Data files
 *.bin
 *.tdf
+*.tgz
 
 
 

common/utils.py

 import numpy as np
+import scipy.signal
 from dataclasses import dataclass
+import tarfile
+import re
 
 
 def avg(*data, n=100):
@@ -36,9 +39,10 @@ class LassieSpillData(Trace):
     unit: str = 'a.u.'
     
     @classmethod
-    def from_file(cls, fname, calibrated=True, from_event=None):
+    def from_file(cls, fname, calibrated=True, from_event=None, verbose=1):
         """Read in time series data from a *.tdf file saved with lassiespill application
            
+        :param fname: path to filename
         :param calibrated: use calibrated or raw data
         :param from_event: return data from this event, time will also be relative to this event                    
         
@@ -70,13 +74,14 @@ class LassieSpillData(Trace):
                 
             elif bi.get_title() == 'TraceInfo':
                 for row in block.get_rows():
+                    if verbose: print(row)
                     if row.get_tag() == 'sampling frequency':
                         assert row.get_unit().lower() == 'hz', 'sampling frequency has unexpected unit {}'.format(row.get_unit())
                         fs = row.get_value()
                         
             elif bi.get_title() == 'Integrals':
                 for row in block.get_rows():
-                    print('>'+row.get_tag()+'<')
+                    if verbose: print(row)
                     if row.get_tag() == ('Calibrated' if calibrated else 'Raw') + ': beam in integral':
                         unit = row.get_unit()
                         
@@ -108,7 +113,60 @@ class LassieSpillData(Trace):
 # IPM data
 ###########
 
-... # TODO
+@dataclass
+class IPMData(Trace):
+    """Container for profile data from IPM
+    
+    t - time array in seconds
+    p - profile position array
+    p_unit - unit of displacement array values
+    x - horizontal profile amplitudes
+    y - vertical profile amplitudes
+    unit - unit of profile amplitudes
+    """
+    t: np.ndarray
+    p: np.ndarray
+    x: np.ndarray
+    y: np.ndarray
+    unit: str = 'a.u.'
+    p_unit: str = 'mm'
+    
+    @classmethod
+    def from_file(cls, fname, clean=False, verbose=1):
+        """Read in profile data from a *.tgz file saved with IPM application
+           
+        :param fname: path to filename
+        :param clean: use adc_clean.dat instead of adc.dat
+        :param from_event: return data from this event, time will also be relative to this event
+        
+        References: Giacomini, Tino <T.Giacomini@gsi.de>
+                One beam profile consist of 64 values.
+                The distance between two wires is 0.6 mm. The wires diameter is 1.5 mm. So the center to center distance of the wires is 2.1 mm.
+        
+        """
+        tar = tarfile.open(fname, "r:gz")
+        
+        # profile data
+        f = tar.extractfile('adc_clean.dat' if clean else 'adc.dat')
+        dat = np.loadtxt(f, skiprows=1)
+        p = 2.1 * (np.arange(64)-31.5)
+        x = dat[:,1].reshape((-1, 64))
+        y = dat[:,2].reshape((-1, 64))
+        
+        # settings
+        dt = 0
+        f = tar.extractfile('param.dat')
+        for line in f:
+            if verbose: print(line.decode('ascii').strip())
+            m = re.match(r'intt\s*=\s*(\d+)', 'intt =  0')
+            if m:
+                interval = int(m.group(1))
+                if interval == 0: dt = 5e-3 # 5 ms
+                #if interval == 1: dt = 5e-4 # 0.5 ms
+        assert dt, 'Could not determine time interval of profiles'
+        t = dt * np.arange(x.shape[0])         
+        
+        return IPMData(t, p, x, y)
 
 
 
@@ -126,16 +184,16 @@ class LiberaBBBData(Trace):
     x: np.ndarray
     y: np.ndarray
     s: np.ndarray
-    f: np.ndarray
     x_unit: str = 'mm'
     y_unit: str = 'mm'
     s_unit: str = 'a.u.'
-    f_unit: str = 'Hz'
+    
     
     @classmethod
     def from_file(cls, fname):
         """Read in bunch-by-bunch data from a *.bin file saved with libera-ireg
         
+        :param fname: path to filename
         References: Libera_Hadron_User_Manual_1.04.pdf
         """
         bunch = np.memmap(fname, dtype=np.dtype([
@@ -148,9 +206,7 @@ class LiberaBBBData(Trace):
                 ('t1', '<u2'),
                 ('status', '<u4')]))
         time = (bunch['TS'] + bunch['t1'])*4e-9 # in s
-        freq = 1/np.diff(time)
-        freq = np.concatenate((freq, [freq[-1]]))
-        return LiberaBBBData(time, bunch['X']*1e-6, bunch['Y']*1e-6, bunch['S'], freq)
+        return LiberaBBBData(time, bunch['X']*1e-6, bunch['Y']*1e-6, bunch['S'])
 
     
     

test/test.py

@@ -15,10 +15,9 @@ mpl.rcParams['axes.prop_cycle'] = mpl.cycler(color=CB_color_cycle)
     
     
 
-fig, ax = plt.subplots(5,1, sharex='col', figsize=(5,8))
+fig, ax = plt.subplots(6,1, sharex='col', figsize=(6,8), constrained_layout=True)
 
-data = LassieSpillData.from_file('test_2022-05-04_13 16 50_asl340.acc.gsi.de_GS01DT_ML_pniederm_test_42c_00014.tdf',
-                                 from_event=CMD_BEAM_INJECTION)
+data = LassieSpillData.from_file('test_2022-05-04_13 16 50_asl340.acc.gsi.de_GS01DT_ML_pniederm_test_42c_00014.tdf', from_event=CMD_BEAM_INJECTION)
 ax[0].plot(data.t, data.v)
 ax[0].set(title='BCT', ylabel=data.unit)
 
@@ -29,9 +28,18 @@ ax[1].plot(data.t, data.y, label=f'y / {data.y_unit}')
 ax[1].legend()
 ax[1].set(title='BPM', ylabel=data.x_unit)
 ax[2].plot(data.t, data.s)
-ax[2].set(title='BPM sum', ylabel=data.s_unit)
-ax[3].plot(*avg(data.t, data.f, n=1000))
-ax[3].set(title='BPM bunch frequency', ylabel=data.f_unit)
+ax[2].set(title='BPM sum', ylabel=data.s_unit)        
+ax[3].plot(*avg(data.t[:-1], 1/np.diff(data.t), n=1000))
+ax[3].set(title='BPM bunch frequency', ylabel='Hz')
+
+
+data = IPMData.from_file('ipmData_SIS1809_220427_111520.tgz', clean=False)
+cm = ax[4].pcolormesh(data.t, data.p, data.x.T, cmap='gist_ncar_r')
+ax[4].set(title='IPM X', ylabel=data.p_unit)
+fig.colorbar(cm, ax=ax[4])
+cm = ax[5].pcolormesh(data.t, data.p, data.y.T, cmap='gist_ncar_r')
+ax[5].set(title='IPM Y', ylabel=data.p_unit)
+fig.colorbar(cm, ax=ax[5])
 
 
 ax[-1].set(xlabel='time / s')