docs: 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)