diff --git a/xgm.py b/xgm.py index 5f7c3c092e51fa3e31c6df7e9613d865d188450d..c9b717c864e92c485511302c823df048fc729f42 100644 --- a/xgm.py +++ b/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