xgm.py

@@ -535,7 +535,7 @@ def calibrateTIM(data, rollingWindow=200, mcp=1, plot=False, use_apd=True, intst
         
         ax = plt.subplot(234)
         xgm_fast = selectSASEinXGM(data)
-        ax.scatter(filteredTIM, xgm_fast, s=5, alpha=0.1)
+        ax.scatter(filteredTIM, xgm_fast, s=5, alpha=0.1, rasterize=True)
         fit, cov = np.polyfit(filteredTIM.values.flatten(),xgm_fast.values.flatten(),1, cov=True)
         y=np.poly1d(fit)
         x=np.linspace(filteredTIM.min(), filteredTIM.max(), 10)
@@ -644,3 +644,81 @@ def timFactorFromTable(voltage, photonEnergy, mcp=1):
     poly = np.poly1d(tim_calibration_table[photonEnergy][mcp-1])
     f = -np.exp(poly(voltage))
     return f
+
+
+def checkTimApdWindow(data, mcp=1, use_apd=True, intstart=None, intstop=None):
+    ''' Plot the first and last pulses in MCP trace together with 
+        the window of integration to check if the pulse integration
+        is properly calculated. If the number of pulses changed during
+        the run, it selects a train where the number of pulses was 
+        maximum.
+        
+        Inputs:
+            data: xarray Dataset
+            mcp: MCP channel (1, 2, 3 or 4)
+            use_apd: if True, gets the APD parameters from the digitizer
+                device. If False, uses intstart and intstop as boundaries
+                and uses the bunch pattern to determine the separation
+                between two pulses.
+            intstart: trace index of integration start of the first pulse
+            intstop: trace index of integration stop of the first pulse
+            
+        Output:
+            Plot    
+    '''
+    npulses_max = data['npulses_sase3'].max().values
+    tid = data['npulses_sase3'].where(data['npulses_sase3'] == npulses_max,
+                                      drop=True)[0].trainId.values
+    if 'MCP{}raw'.format(mcp) not in data:
+        tid, data_from_train = data.attrs['run'].train_from_id(tid)
+        trace = data_from_train['SCS_UTC1_ADQ/ADC/1:network']['digitizers.channel_1_D.raw.samples']
+        print('no raw data for MCP{}. Loading trace from MCP1'.format(mcp))
+        label_trace='MCP1 Voltage [V]'
+    else:
+        idx = np.argwhere(data['MCP{}raw'.format(mcp)].trainId.values == tid)[0]
+        trace = data['MCP{}raw'.format(mcp)][idx].T
+        label_trace='MCP{} Voltage [V]'.format(mcp)
+    if use_apd:
+        pulseStart = data.attrs['run'].get_array(
+            'SCS_UTC1_ADQ/ADC/1', 'board1.apd.channel_0.pulseStart.value')[0].values
+        pulseStop = data.attrs['run'].get_array(
+            'SCS_UTC1_ADQ/ADC/1', 'board1.apd.channel_0.pulseStop.value')[0].values
+        initialDelay = data.attrs['run'].get_array(
+                        'SCS_UTC1_ADQ/ADC/1', 'board1.apd.channel_0.initialDelay.value')[0].values
+        upperLimit = data.attrs['run'].get_array(
+                        'SCS_UTC1_ADQ/ADC/1', 'board1.apd.channel_0.upperLimit.value')[0].values
+        nsamples = upperLimit - initialDelay
+    else:
+        pulseStart = intstart
+        pulseStop = intstop
+        if npulses_max > 1:
+            sa3 = data['sase3'].where(data['sase3']>1)
+            step = sa3.where(data['npulses_sase3']>1, drop=True)[0,:2].values
+            step = int(step[1] - step[0])
+            nsamples = 440 * step
+        else:
+            nsamples = 0
+
+    fig, ax = plt.subplots(figsize=(5,3))
+    ax.plot(trace[:pulseStop+25], color='C1', label='first pulse')
+    ax.axvspan(pulseStart, pulseStop, color='k', alpha=0.1, label='APD region')
+    ax.axvline(pulseStart, color='gray', ls='--')
+    ax.axvline(pulseStop, color='gray', ls='--')
+    ax.set_xlim(pulseStart-25, pulseStop+25)
+    ax.locator_params(axis='x', nbins=4)
+    ax.set_ylabel(label_trace)
+    ax.set_xlabel('First pulse sample #')
+    if npulses_max > 1:
+        pulseStart = pulseStart + nsamples*(npulses_max-1)
+        pulseStop = pulseStop + nsamples*(npulses_max-1)
+        ax2 = ax.twiny()
+        ax2.plot(range(pulseStart-25,pulseStop+25), trace[pulseStart-25:pulseStop+25],
+                color='C4', label='last pulse')
+        ax2.locator_params(axis='x', nbins=4)
+        ax2.set_xlabel('Last pulse sample #')
+        lines, labels = ax.get_legend_handles_labels()
+        lines2, labels2 = ax2.get_legend_handles_labels()
+        ax2.legend(lines + lines2, labels + labels2, loc=0)
+    else:
+        ax.legend(loc='lower left')
+    plt.tight_layout()
\ No newline at end of file