def writeColorLUT2(config,
                  outfile=None, isochrone=None, distance_modulus_array=None,
                  delta_mag=None, mag_err_array=None,
                  mass_steps=10000, plot=False):
    """
    Precompute a 4-dimensional signal color probability look-up table to speed up the likelihood evaluation.
    Inputs are a Config object (or file name), an Isochrone object,
    an array of distance moduli at which to evaluate the signal color probability,
    and an array of magnitude uncertainties which set the bin edges of those dimensions (zero implicity included).
    Finally there is an outfile name.
    """
    if plot: import ugali.utils.plotting

    if type(config) == str:
        config = ugali.utils.config.Config(config)
    if outfile is None:
        outfile = config.params['color_lut']['filename']
    if isochrone is None:
        isochrones = []
        for ii, name in enumerate(config.params['isochrone']['infiles']):
            isochrones.append(ugali.isochrone.Isochrone(config, name))
        isochrone = ugali.isochrone.CompositeIsochrone(isochrones, config.params['isochrone']['weights'])
    if distance_modulus_array is None:
        distance_modulus_array = config.params['color_lut']['distance_modulus_array']
    if delta_mag is None:
        delta_mag = config.params['color_lut']['delta_mag']
    if mag_err_array is None:
        mag_err_array = config.params['color_lut']['mag_err_array']

    mag_buffer = 0.5 # Safety buffer in magnitudes around the color-magnitude space defined by the ROI
    epsilon = 1.e-10
    if config.params['catalog']['band_1_detection']:
        bins_mag_1 = numpy.arange(config.params['mag']['min'] - mag_buffer,
                                  config.params['mag']['max'] + mag_buffer + epsilon,
                                  delta_mag)
        bins_mag_2 = numpy.arange(config.params['mag']['min'] - config.params['color']['max'] - mag_buffer,
                                  config.params['mag']['max'] - config.params['color']['min'] + mag_buffer + epsilon,
                                  delta_mag)
    else:
        bins_mag_1 = numpy.arange(config.params['mag']['min'] + config.params['color']['min'] - mag_buffer,
                                  config.params['mag']['max'] + config.params['color']['max'] + mag_buffer + epsilon,
                                  delta_mag)
        bins_mag_2 = numpy.arange(config.params['mag']['min'] - mag_buffer,
                                  config.params['mag']['max'] + mag_buffer + epsilon,
                                  delta_mag)

    # Output binning configuration
    #print config.params['catalog']['band_1_detection']
    #print config.params['mag']['min'], config.params['mag']['max']
    #print config.params['color']['min'], config.params['color']['max']

    #print bins_mag_1[0], bins_mag_1[-1], len(bins_mag_1)
    #print bins_mag_2[0], bins_mag_2[-1], len(bins_mag_2)

    isochrone_mass_init, isochrone_mass_pdf, isochrone_mass_act, isochrone_mag_1, isochrone_mag_2 = isochrone.sample(mass_steps=mass_steps)

    hdul = pyfits.HDUList()

    for index_distance_modulus, distance_modulus in enumerate(distance_modulus_array):

        logger.debug('(%i/%i)'%(index_distance_modulus, len(distance_modulus_array)))

        columns_array = []
        
        time_start = time.time()

        histo_isochrone_pdf = numpy.histogram2d(distance_modulus + isochrone_mag_1,
                                                distance_modulus + isochrone_mag_2,
                                                bins=[bins_mag_1, bins_mag_2],
                                                weights=isochrone_mass_pdf)[0]
        
        if plot:
            # Checked that axis are plotted correctly
            ugali.utils.plotting.twoDimensionalHistogram('Isochrone', 'mag_1', 'mag_2',
                                                         numpy.log10(histo_isochrone_pdf + epsilon).transpose(),
                                                         bins_mag_1, bins_mag_2,
                                                         lim_x=None, lim_y=None,
                                                         vmin=None, vmax=None)

            
        
        for index_mag_err_1, mag_err_1 in enumerate(mag_err_array):
            for index_mag_err_2, mag_err_2 in enumerate(mag_err_array):
                logger.debug('  Distance modulus = %.2f mag_err_1 = %.2f mag_err_2 = %.2f'%(distance_modulus, mag_err_1, mag_err_2))

                mag_1_sigma_step = delta_mag / mag_err_1
                n = int(numpy.ceil(4. / mag_1_sigma_step))
                mag_1_sigma = numpy.arange(-1. * (n + 0.5) * mag_1_sigma_step,
                                           ((n + 0.5) * mag_1_sigma_step) + epsilon,
                                           mag_1_sigma_step)
                mag_1_pdf_array = scipy.stats.norm.cdf(mag_1_sigma[1:]) - scipy.stats.norm.cdf(mag_1_sigma[0:-1])

                mag_2_sigma_step = delta_mag / mag_err_2
                n = int(numpy.ceil(4. / mag_2_sigma_step))
                mag_2_sigma = numpy.arange(-1. * (n + 0.5) * mag_2_sigma_step,
                                           ((n + 0.5) * mag_2_sigma_step) + epsilon,
                                           mag_2_sigma_step)
                mag_2_pdf_array = scipy.stats.norm.cdf(mag_2_sigma[1:]) - scipy.stats.norm.cdf(mag_2_sigma[0:-1])

                mag_1_pdf, mag_2_pdf = numpy.meshgrid(mag_2_pdf_array, mag_1_pdf_array)
                
                pdf = mag_1_pdf * mag_2_pdf

                histo_isochrone_pdf_convolve = scipy.signal.convolve2d(histo_isochrone_pdf, pdf, mode='same')

                if plot:
                    # Checked that axis are plotted correctly
                    ugali.utils.plotting.twoDimensionalHistogram('Convolved Isochrone', 'mag_1', 'mag_2',
                                                                 numpy.log10(histo_isochrone_pdf_convolve + epsilon).transpose(),
                                                                 bins_mag_1, bins_mag_2,
                                                                 lim_x=None, lim_y=None,
                                                                 vmin=None, vmax=None)

                columns_array.append(pyfits.Column(name = '%i%i'%(index_mag_err_1, index_mag_err_2),
                                                   format = '%iE'%(histo_isochrone_pdf_convolve.shape[1]),
                                                   array = histo_isochrone_pdf_convolve))

        hdu = pyfits.new_table(columns_array)
        hdu.header.update('DIST_MOD', distance_modulus)
        hdu.name = '%.2f'%(distance_modulus)
        hdul.append(hdu)

        time_end = time.time()
        logger.debug('%.2f s'%(time_end - time_start))

    # Store distance modulus info
    columns_array = [pyfits.Column(name = 'DISTANCE_MODULUS',
                                   format = 'E',
                                   array = distance_modulus_array)]
    hdu = pyfits.new_table(columns_array)
    hdu.name = 'DISTANCE_MODULUS'
    hdul.append(hdu)

    # Store magnitude error info
    columns_array = [pyfits.Column(name = 'BINS_MAG_ERR',
                                   format = 'E',
                                   array = numpy.insert(mag_err_array, 0, 0.))]
    hdu = pyfits.new_table(columns_array)
    hdu.name = 'BINS_MAG_ERR'
    hdul.append(hdu)

    # Store magnitude 1 info
    columns_array = [pyfits.Column(name = 'BINS_MAG_1',
                                   format = 'E',
                                   array = bins_mag_1)]
    hdu = pyfits.new_table(columns_array)
    hdu.name = 'BINS_MAG_1'
    hdul.append(hdu)

    # Store magnitude 2 info
    columns_array = [pyfits.Column(name = 'BINS_MAG_2',
                                   format = 'E',
                                   array = bins_mag_2)]
    hdu = pyfits.new_table(columns_array)
    hdu.name = 'BINS_MAG_2'
    hdul.append(hdu)

    logger.info('Writing look-up table to %s'%(outfile))
    hdul.writeto(outfile, clobber = True)