docs: def cutIsochronePath(g, r, g_err, r_err, isochrone, radius=0.1, return_all=False):
      """
      Cut to identify objects within isochrone cookie-cutter.
  
      ADW: This should be moved into the isochrone class.
      """
      import scipy.interpolate
      from ugali.isochrone import CompositeIsochrone
  
      if isinstance(isochrone, CompositeIsochrone):
          isochrone = isochrone.isochrones[0]
  
      if len(g) == 0:
          return np.array([],dtype=bool)
  
      try:
          if np.all(isochrone.stage == 'Main'):
              # Dotter case
              index_transition = len(isochrone.stage)
          else:
              # Other cases
              index_transition = np.nonzero(isochrone.stage > 3)[0][0] + 1    
      except AttributeError:
          index_transition = 1
  
      mag_1_rgb = isochrone.mag_1[0: index_transition] + isochrone.distance_modulus
      mag_2_rgb = isochrone.mag_2[0: index_transition] + isochrone.distance_modulus
      mag_1_rgb = mag_1_rgb[::-1]
      mag_2_rgb = mag_2_rgb[::-1]
      
      # Cut one way...
      f_isochrone = scipy.interpolate.interp1d(mag_2_rgb, mag_1_rgb - mag_2_rgb, bounds_error=False, fill_value = 999.)
      color_diff = np.fabs((g - r) - f_isochrone(r))
      cut_2 = (color_diff < np.sqrt(0.1**2 + r_err**2 + g_err**2))
  
       # ...and now the other
      f_isochrone = scipy.interpolate.interp1d(mag_1_rgb, mag_1_rgb - mag_2_rgb, bounds_error=False, fill_value = 999.)
      color_diff = np.fabs((g - r) - f_isochrone(g))
      cut_1 = (color_diff < np.sqrt(0.1**2 + r_err**2 + g_err**2))
  
      cut = np.logical_or(cut_1, cut_2)
      
      # Include horizontal branch if it exists
      if not np.any(isochrone.stage == isochrone.hb_stage):
          index_transition = np.nonzero(isochrone.stage==isochrone.hb_stage)[0][0]+1
          mag_1_hb = isochrone.mag_1[index_transition:] + isochrone.distance_modulus
          mag_2_hb = isochrone.mag_2[index_transition:] + isochrone.distance_modulus
          path_hb = makePath(mag_1_hb, mag_2_hb)
          cut_hb = path_hb.contains_points(list(zip(g, r)), radius=0.1)
          logger.debug('Applying HB selection')
          logger.debug(np.sum(cut))
          cut = np.logical_or(cut, cut_hb)
          logger.debug(np.sum(cut))
  
      mag_bins = np.arange(16., 24.1, 0.1)
      mag_centers = 0.5 * (mag_bins[1:] + mag_bins[0:-1])
      magerr = np.tile(0., len(mag_centers))
      for ii in range(0, len(mag_bins) - 1):
          cut_mag_bin = (g > mag_bins[ii]) & (g < mag_bins[ii + 1])
          magerr[ii] = np.median(np.sqrt(0.1**2 + r_err[cut_mag_bin]**2 + g_err[cut_mag_bin]**2))
  
      if return_all:
          return cut, mag_centers[f_isochrone(mag_centers) < 100], (f_isochrone(mag_centers) + magerr)[f_isochrone(mag_centers) < 100], (f_isochrone(mag_centers) - magerr)[f_isochrone(mag_centers) < 100]
      else:
          return cut