Add TDC data plotting

plotting.py

@@ -485,3 +485,96 @@ def plot_beam_size(ax, ipm_data, xy, time_range=None, smoothing=None, **kwargs):
     ls, = ax.plot(ipm_data.t[mask], v, **kwargs)
     ax.fill_between(ipm_data.t[mask], v-ve, v+ve, alpha=0.3, color=ls.get_color(), zorder=-10)    
     ax.set(ylabel=f'Beam size $\\sigma_{xy}$ / mm', xlabel='Time / s')
+
+    
+
+# TDC data
+###########
+
+def plot_spill_intensity(ax, spill_data, bin_time=10e-6, *, cumulative=False, rate=False):
+    """Plot spill intensity (binned particle counts) over extraction time
+    
+    :param ax: Axis to plot onto
+    :param spill_data: Instance of TDCSpill
+    :param bin_time: width of bins for particle counting in sec
+    """
+    nbins = int(np.ceil(spill_data.length/bin_time))
+    weights = np.ones_like(spill_data.hittimes)/bin_time if rate else None
+    ax.hist(spill_data.hittimes, range=(spill_data.time_offset, spill_data.time_offset+bin_time*nbins), bins=nbins,
+            histtype='step', cumulative=cumulative, weights=weights)
+    ax.set(xlabel='Extraction time / s', xlim=(spill_data.time_offset, spill_data.time_offset+bin_time*nbins),
+           ylabel='Extracted particles' + ('' if cumulative else ' / s' if rate else f' / ${(bin_time*SI.s).to_compact():~gL}$ interval'))
+
+
+def plot_spill_duty_factor(ax, spill_data, count_bin_time=10e-6, evaluate_bin_time=10e-3, *, show_poisson_limit=False, **kwargs):
+    """Plot spill intensity (binned particle counts) over extraction time
+    
+    :param ax: Axis to plot onto
+    :param spill_data: Instance of TDCSpill
+    :param count_bin_time: width of bins for particle counting in sec
+    :param evaluate_bin_time: width of bins for duty factor evaluation (value is rounded to the next multiple of count_bin_time)
+                              
+    :param show_poisson_limit: if true, show poison limit of duty factor
+    
+    """
+    cbins = int(np.ceil(spill_data.length/count_bin_time))
+    hist = np.histogram(spill_data.hittimes, bins=cbins, 
+                        range=(spill_data.time_offset, spill_data.time_offset+count_bin_time*cbins))
+    ebins = int(round(evaluate_bin_time/count_bin_time))
+    counts = hist[0][:int(len(hist[0])/ebins)*ebins].reshape((-1, ebins))
+    edges = hist[1][:int(len(hist[0])/ebins+1)*ebins:ebins]
+    
+    # spill duty factor
+    F = np.mean(counts, axis=1)**2 / np.mean(counts**2, axis=1)
+    ax.stairs(F, edges, **kwargs)
+    
+    if show_poisson_limit:
+        Fp = np.mean(counts, axis=1) / ( np.mean(counts, axis=1) + 1 )
+        ax.stairs(Fp, edges, color='gray', label='Poisson limit', zorder=-1)
+        
+    ax.set(xlabel='Extraction time / s', ylim=(0,1.1),
+           ylabel=f'Spill duty factor F\n(${(evaluate_bin_time*SI.s).to_compact():~gL}$ / ${(count_bin_time*SI.s).to_compact():~gL}$ intervals)'
+    )
+
+
+def plot_spill_consecutive_particle_delay(ax, spill_data, limit=5e-6, resolution=10e-9):
+    """Plot consecutive particle separation (particle-interval histograms)
+    
+    :param ax: Axis to plot onto
+    :param spill_data: Instance of TDCSpill
+    :param limit: maximum delay to plot (range of histogram)
+    :param resolution: time resolution (bin width of histogram)
+    """
+    nbins = int(np.ceil(limit/resolution))
+    delay = np.diff(spill_data.hittimes)
+    ax.hist(delay, range=(0, resolution*nbins), bins=nbins, histtype='step') #, weights=np.ones_like(delay)/resolution)
+    ax.set(xlabel='Delay between consecutive particles / s', xlim=(resolution, resolution*nbins),
+           ylabel=f'particle count / ${(resolution*SI.s).to_compact():~gL}$ delay interval')
+
+    #TODO: plot poison statistics
+
+
+def plot_spill_frequency_spectrum(ax, spill_data, fmax=100e3):
+    """Plot frequency spectrum of particle arrival times
+    
+    :param ax: Axis to plot onto
+    :param spill_data: Instance of TDCSpill
+    :param fmax: maximum frequency for analysis in Hz
+    :param bins: number of bins for fft
+    """
+    
+    # re-sample timestamps into equally sampled time series
+    dt = 1/2/fmax
+    bins = int(spill_data.length/dt)
+    timeseries, _ = np.histogram(spill_data.hittimes, bins=bins, range=(spill_data.time_offset, spill_data.time_offset + bins*dt))
+    
+    # frequency analysis
+    freq = np.fft.rfftfreq(len(timeseries), d=dt)[1:]
+    mag = np.abs(np.fft.rfft(timeseries))[1:]
+    
+    # plot
+    ax.plot(freq, mag)
+    ax.set(xlabel=f'Frequency / Hz', xlim=(0, fmax),
+           ylabel='Magnitude / a.u.', ylim=(0, 1.1*max(mag)))
+
+