code: def add_to_twofold_matrix(times, tdcs, mat, tmax=10):
      
      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
docs: 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)))