Cut to identify objects within isochrone cookie-cutter.

    ADW: This should be moved into the isochrone class.

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