swabian_nfold_coincidences.py

def find_nfold_coincidences(timestamps, channels, nfold = 4, max_coincidence_time_ps = 20e3):
    """ Given a list of timestamps, and a list of channels those timestamps
    appeared on (both of length N), finds nfold coincidences.  Coincidences are
    limited to cases where:
        (1) There are `nfold` timestamps within a timeframe `max_coincidence_time_ps`
        (2) Those timestamps all appear on different channels """

    td = np.diff(timestamps)
    long_delays_idx = np.where(td > max_coincidence_time_ps)[0]
    coincidence_t = np.hsplit(timestamps, long_delays_idx+1)
    coincidence_c = np.hsplit(channels, long_delays_idx+1)
    coincidences = []
    for n,t in enumerate(coincidence_t):
        c = coincidence_c[n]
        # Check if (1) There are nfold close-together timestamps
        #          (2) Those timestamps all appear on different channels
        if (len(t) == nfold) and (len(np.unique(c)) == nfold):
            coincidences.append(t[np.argsort(c)])
    coincidences = np.array(coincidences)
    return coincidences

def fourfold_coincidences_to_xy(
    coincidences,
    bbox = (-np.inf, np.inf, -np.inf, np.inf) # left right bottom top
    ):
    c = coincidences
    xc = c[:,0] - c[:,1]
    yc = c[:,2] - c[:,3]
    sel_x = (bbox[0] < xc) & (xc < bbox[1])
    sel_y = (bbox[2] < yc) & (yc < bbox[3])
    sel = sel_x & sel_y

    return x[sel], y[sel]

timestamps = data['timestamps']
channels = np.abs(np.array(data['channels']))
coincidences = find_nfold_coincidences(timestamps, channels, nfold = 4, max_coincidence_time_ps = 20e-9*1e12)
xc,yc = fourfold_coincidences_to_xy(coincidences, bbox = (-4000, 9000, -4000, 9000))
plt.figure()
plt.plot(xc, yc,'.', alpha = 0.3)