DSSC1module.py

@@ -106,7 +106,7 @@ class DSSC1module:
         except:
             self.delay_mm = 0*self.nrj
         self.t0 = t0
-        self.delay_ps = tb.positionToDelay(self.delay_mm, origin=self.t0)
+        self.delay_ps = tb.positionToDelay(self.delay_mm, origin=self.t0, invert=True)
                     
     def collect_dssc_module_file(self):
         """ Collect the raw DSSC module h5 files.

Load.py

@@ -63,6 +63,10 @@ mnemonics = {
     "UND": {'source':'SA3_XTD10_UND/DOOCS/PHOTON_ENERGY',
                 'key':'actualPosition.value',
                 'dim':None},
+    #DPS imagers
+    "DPS2CAM2": {'source':'SCS_BLU_DPS-2/CAM/IMAGER2CAMERA:daqOutput',
+                'key':'data.image.pixels',
+                'dim':['dps2cam2_y', 'dps2cam2_x']},
     # XTD10 XGM
     ## keithley
     "XTD10_photonFlux": {'source':'SA3_XTD10_XGM/XGM/DOOCS',
@@ -209,6 +213,10 @@ mnemonics = {
     "PP800_TeleLens": {'source':'SCS_ILH_LAS/MOTOR/LT7',
                  'key':'actualPosition.value',
                  'dim':None},
+    "ILH_8CAM1": {'source':'SCS_ILH_LAS/CAM/8CAM1:daqOutput',
+                'key':'data.image.pixels',
+                'dim':['8cam1_y', '8cam1_x']},
+
     
     # FFT
     "scannerX": {'source':'SCS_CDIFFT_SAM/LMOTOR/SCANNERX',
@@ -229,6 +237,10 @@ mnemonics = {
     "magnet_old": {'source':'SCS_CDIFFT_MAG/SUPPLY/CURRENT',
                'key':'actualCurrent.value',
                'dim':None},
+    
+    "Vertical_FDM": {'source':'SCS_CDIFFT_LDM/CAM/CAMERA1A:daqOutput',
+                'key':'data.image.pixels',
+                'dim':['vfdm_y', 'vfdm_x']},
 
     # FastCCD, if in raw folder, raw images
     #          if in proc folder, dark substracted and relative gain corrected
@@ -471,6 +483,6 @@ def concatenateRuns(runs):
             return
     
     result = xr.concat(orderedRuns, dim='trainId')
-    result.attrs['run'] = [run.attrs['run'] for run in orderedRuns]
-    result.attrs['runFolder'] = [run.attrs['runFolder'] for run in orderedRuns]
+    for k in orderedRuns[0].attrs.keys():
+        result.attrs[k] = [run.attrs[k] for run in orderedRuns]
     return result

XAS.py

@@ -12,6 +12,7 @@
 import numpy as np
 import matplotlib.gridspec as gridspec
 import matplotlib.pyplot as plt
+import re
 
 def absorption(T, Io):
     """ Compute the absorption A = -ln(T/Io)
@@ -190,9 +191,9 @@ def xas(nrun, bins=None, Iokey='SCS_SA3', Itkey='MCP3apd', nrjkey='nrj', Iooffse
                 color='C1', alpha=0.2)
         ax1_twin.set_ylabel('Io')
         try:
-            proposalNB=int(nrun.attrs['runFolder'].split('/')[-4][1:])
-            runNB=int(nrun.attrs['runFolder'].split('/')[-2][1:])
-            ax1.set_title('run {:d} p{:}'.format(runNB, proposalNB))
+            proposalNB=int(re.findall(r'p(\d{6})', nrun.attrs['runFolder'])[0])
+            runNB=int(re.findall(r'r(\d{4})', nrun.attrs['runFolder'])[0])
+            ax1.set_title(f'run {runNB} p{proposalNB}')
         except:
             f.suptitle(nrun.attrs['plot_title'])
         

__init__.py

@@ -6,3 +6,4 @@ from ToolBox.Laser_utils import *
 from ToolBox.DSSC import DSSC
 from ToolBox.azimuthal_integrator import *
 from ToolBox.DSSC1module import *
+from ToolBox.FastCCD import *

azimuthal_integrator.py

 import numpy as np
 
 class azimuthal_integrator(object):
-    def __init__(self, imageshape, center, polar_range, dr=2):
+    def __init__(self, imageshape, center, polar_range, dr=2, aspect=204/236):
         '''
         Create a reusable integrator for repeated azimuthal integration of similar
         images. Calculates array indices for a given parameter set that allows
@@ -21,6 +21,9 @@ class azimuthal_integrator(object):
         dr : int, default 2
             radial width of the integration slices. Takes non-square DSSC pixels into account.
         
+        aspect: float, default 204/236 for DSSC
+            aspect ratio of the pixel pitch
+        
         Returns
         =======
         ai : azimuthal_integrator instance
@@ -39,7 +42,7 @@ class azimuthal_integrator(object):
         ycoord -= cy
 
         # distance from center, hexagonal pixel shape taken into account
-        dist_array = np.hypot(xcoord * 204 / 236, ycoord)
+        dist_array = np.hypot(xcoord * aspect, ycoord)
 
         # array of polar angles
         if np.abs(polar_range[1]-polar_range[0]) > 180:
@@ -53,7 +56,7 @@ class azimuthal_integrator(object):
         polar_array = np.mod(polar_array, np.pi)
         self.polar_mask = (polar_array > tmin) * (polar_array < tmax)
 
-        self.maxdist = min(sx  - cx, sy  - cy)
+        self.maxdist = max(sx  - cx, sy  - cy)
 
         ix, iy = np.indices(dimensions=(sx, sy))
         self.index_array = np.ravel_multi_index((ix, iy), (sx, sy))

knife_edge.py

@@ -8,10 +8,12 @@ import matplotlib.pyplot as plt
 import numpy as np
 from scipy.special import erfc
 from scipy.optimize import curve_fit
+import bisect
 
-def knife_edge(nrun, axisKey='scannerX', signalKey='FastADC4peaks', p0=None, full=False, plot=False):
+def knife_edge(nrun, axisKey='scannerX', signalKey='FastADC4peaks', 
+               axisRange=[None,None], p0=None, full=False, plot=False):
     ''' Calculates the beam radius at 1/e^2 from a knife-edge scan by fitting with
-        erfc function: f(a, u) = a*erfc(u) or f(a, u) = a*erfc(-u) where 
+        erfc function: f(a,b,u) = a*erfc(u)+b or f(a,b,u) = a*erfc(-u)+b where 
         u = sqrt(2)*(x-x0)/w0 with w0 the beam radius at 1/e^2 and x0 the beam center.
         Inputs:
             nrun: xarray Dataset containing the detector signal and the motor 
@@ -19,22 +21,24 @@ def knife_edge(nrun, axisKey='scannerX', signalKey='FastADC4peaks', p0=None, ful
             axisKey: string, key of the axis against which the knife-edge is 
                   performed.
             signalKey: string, key of the detector signal.
-            p0: list, initial parameters used for the fit: x0, w0, a. If None, a beam
+            axisRange: list of length 2, minimum and maximum values between which to apply
+                       the fit.
+            p0: list, initial parameters used for the fit: x0, w0, a, b. If None, a beam
                 radius of 100 um is assumed.
             full: bool: If False, returns the beam radius and standard error. If True,
                 returns the popt, pcov list of parameters and covariance matrix from 
-                curve_fit.
+                curve_fit as well as the fitting function.
             plot: bool: If True, plots the data and the result of the fit.
         Outputs:
             If full is False, ndarray with beam radius at 1/e^2 in mm and standard 
             error from the fit in mm. If full is True, returns popt and pcov from 
             curve_fit function.
     '''
-    def integPowerUp(x, x0, w0, a):
-        return a*erfc(-np.sqrt(2)*(x-x0)/w0)
+    def integPowerUp(x, x0, w0, a, b):
+        return a*erfc(-np.sqrt(2)*(x-x0)/w0) + b
 
-    def integPowerDown(x, x0, w0, a):
-        return a*erfc(np.sqrt(2)*(x-x0)/w0)
+    def integPowerDown(x, x0, w0, a, b):
+        return a*erfc(np.sqrt(2)*(x-x0)/w0) + b
 
     #get the number of pulses per train from the signal source:
     dim = nrun[signalKey].dims[1]
@@ -46,22 +50,38 @@ def knife_edge(nrun, axisKey='scannerX', signalKey='FastADC4peaks', p0=None, ful
     sortIdx = np.argsort(positions)
     positions = positions[sortIdx]
     intensities = nrun[signalKey].values.flatten()[sortIdx]
-
+    
+    if axisRange[0] is None or axisRange[0] < positions[0]:
+        idxMin = 0
+    else:
+        if axisRange[0] >= positions[-1]:
+            raise ValueError('The minimum value of axisRange is too large')
+        idxMin = bisect.bisect(positions, axisRange[0])
+    if axisRange[1] is None or axisRange[1] > positions[-1]:
+        idxMax = None
+    else:
+        if axisRange[1] <= positions[0]:
+            raise ValueError('The maximum value of axisRange is too small')
+        idxMax = bisect.bisect(positions, axisRange[1]) + 1
+    positions = positions[idxMin:idxMax]
+    intensities = intensities[idxMin:idxMax]
+       
     # estimate a linear slope fitting the data to determine which function to fit
     slope = np.cov(positions, intensities)[0][1]/np.var(positions) 
     if slope < 0:
         func = integPowerDown
-        funcStr = 'a*erfc(np.sqrt(2)*(x-x0)/w0)'
+        funcStr = 'a*erfc(np.sqrt(2)*(x-x0)/w0) + b'
     else:
         func = integPowerUp
-        funcStr = 'a*erfc(-np.sqrt(2)*(x-x0)/w0)'
+        funcStr = 'a*erfc(-np.sqrt(2)*(x-x0)/w0) + b'
     if p0 is None:
-        p0 = [np.mean(positions), 0.1, np.max(intensities)/2]
+        p0 = [np.mean(positions), 0.1, np.max(intensities)/2, 0]
     popt, pcov = curve_fit(func, positions, intensities, p0=p0)
     print('fitting function:', funcStr)
     print('w0 = (%.1f +/- %.1f) um'%(popt[1]*1e3, pcov[1,1]**0.5*1e3))
-    print('x0 = (%.3f +/- %.3f) mm'%(popt[0], pcov[0,0]**0.5*1e3))
-    print('a = %e +/- %e '%(popt[2], pcov[2,2]**0.5*1e3))
+    print('x0 = (%.3f +/- %.3f) mm'%(popt[0], pcov[0,0]**0.5))
+    print('a = %e +/- %e '%(popt[2], pcov[2,2]**0.5))
+    print('b = %e +/- %e '%(popt[3], pcov[3,3]**0.5))
     
     if plot:
         xfit = np.linspace(positions.min(), positions.max(), 1000)
@@ -78,6 +98,6 @@ def knife_edge(nrun, axisKey='scannerX', signalKey='FastADC4peaks', p0=None, ful
         plt.title(nrun.attrs['runFolder'])
         plt.tight_layout()
     if full:
-        return popt, pcov
+        return popt, pcov, func
     else:
         return np.array([popt[1], pcov[1,1]**0.5])

xgm.py

@@ -8,6 +8,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import xarray as xr
+from scipy.signal import find_peaks
 
 # Machine
 def pulsePatternInfo(data, plot=False):
@@ -573,7 +574,7 @@ def getTIMapd(data, mcp=1, use_apd=True, intstart=None, intstop=None,
             print(f'Warning: apd parameter was set to record 1 pulse out of {period} @ 4.5 MHz ' +
                   f'but XFEL delivered 1 pulse out of {period_from_bunch_pattern}.')
         maxPulses = data['npulses_sase3'].max().values
-        if period*pulseIdDim < period_from_bunch_pattern*maxPulses:
+        if period*pulseIdDim < period_from_bunch_pattern*(maxPulses-1):
             print(f'Warning: Number of pulses and/or rep. rate in apd parameters were set ' +
                   f'too low ({pulseIdDim})to record the {maxPulses} SASE 3 pulses')
         peaks = data[f'MCP{mcp}apd']
@@ -967,6 +968,51 @@ def fastAdcPeaks(data, channel, intstart, intstop, bkgstart, bkgstop, period=Non
         results[:,i] = integ
     return results
 
+def autoFindFastAdcPeaks(data, channel=5, threshold=35000, display=False, plot=False):
+    ''' Automatically finds positive peaks in channel of Fast ADC trace, assuming
+        a minimum absolute height of 'threshold' counts and a minimum width of 4 
+        samples. The find_peaks function and determination of the peak region and 
+        baseline subtraction is optimized for typical photodiode signals of the 
+        SCS instrument (ILH, FFT reflectometer, FFT diag stage).
+        Inputs:
+            data: xarray Dataset containing Fast ADC traces
+            key: data key of the array of traces
+            threshold: minimum height of the peaks
+            display: bool, displays info on the pulses found
+            plot: plots regions of integration of the first pulse in the trace
+        Output:
+            peaks: DataArray of the integrated peaks 
+    '''
+    
+    key = f'FastADC{channel}raw'
+    if key not in data:
+        raise ValueError(f'{key} not found in data set')
+    trace = data[key].where(data['npulses_sase3']>0, drop=True).isel(trainId=0).values
+    centers, peaks = find_peaks(trace, height=threshold, width=(4, None))
+    c = centers[0]
+    w = np.average(peaks['widths']).astype(int)
+    period = np.median(np.diff(centers)).astype(int)
+    npulses = centers.shape[0]
+    intstart = int(c - w/4) + 1
+    intstop = int(c + w/4) + 1
+    bkgstop = int(peaks['left_ips'][0])-5
+    bkgstart = bkgstop - 10
+    if display:
+        print(f'Found {npulses} pulses, avg. width={w}, period={period} samples, ' +
+              f'rep. rate={1e6/(9.230769*period):.3f} kHz')
+    fAdcPeaks = fastAdcPeaks(data, channel=channel, intstart=intstart, intstop=intstop,
+                         bkgstart=bkgstart, bkgstop=bkgstop, period=period, npulses=npulses)
+    if plot:
+        plt.figure()
+        plt.plot(trace, 'o-', ms=3)
+        for i in range(npulses):
+            plt.axvline(intstart+i*period, ls='--', color='g')
+            plt.axvline(intstop+i*period, ls='--', color='r')
+            plt.axvline(bkgstart+i*period, ls='--', color='lightgrey')
+            plt.axvline(bkgstop+i*period, ls='--', color='grey')
+        plt.title(f'Fast ADC {channel} trace')
+        plt.xlim(bkgstart-10, intstop + 50)
+    return fAdcPeaks
 
 def mergeFastAdcPeaks(data, channel, intstart, intstop, bkgstart, bkgstop, 
                       period=None, npulses=None, dim='lasPulseId'):