Add counts to twofold coincidences for a given `tmax`.

    Parameters
    ----------
    times: np.ndarray of hit times (int32)
    tdcs: np.ndarray of channel_ids (uint8)
    mat: ref to a np.array((465, tmax * 2 + 1))
    tmax: int (time window)

    Returns
    -------
    mat: coincidence matrix (np.array((465, tmax * 2 + 1)))

def add_to_twofold_matrix(times, tdcs, mat, tmax=10):
    """Add counts to twofold coincidences for a given `tmax`.

    Parameters
    ----------
    times: np.ndarray of hit times (int32)
    tdcs: np.ndarray of channel_ids (uint8)
    mat: ref to a np.array((465, tmax * 2 + 1))
    tmax: int (time window)

    Returns
    -------
    mat: coincidence matrix (np.array((465, tmax * 2 + 1)))

    """
    h_idx = 0    # index of initial hit
    c_idx = 0    # index of coincident candidate hit
    n_hits = len(times)
    multiplicity = 0
    while h_idx <= n_hits:
        c_idx = h_idx + 1
        if (c_idx < n_hits) and (times[c_idx] - times[h_idx] <= tmax):
            multiplicity = 2
            c_idx += 1
            while (c_idx < n_hits) and (times[c_idx] - times[h_idx] <= tmax):
                c_idx += 1
                multiplicity += 1
            if multiplicity != 2:
                h_idx = c_idx
                continue
            c_idx -= 1
            h_tdc = tdcs[h_idx]
            c_tdc = tdcs[c_idx]
            h_time = times[h_idx]
            c_time = times[c_idx]
            if h_tdc != c_tdc:
                dt = int(c_time - h_time)
                if h_tdc > c_tdc:
                    mat[get_comb_index(c_tdc, h_tdc), -dt + tmax] += 1
                else:
                    mat[get_comb_index(h_tdc, c_tdc), dt + tmax] += 1
        h_idx = c_idx