Create a simulation of the background stellar population.
        Because some stars have been clipped to generate the CMD,
        this function tends to slightly underestimate (~1%) the 
        background as compared to the true catalog.

        The simulation of background object colors relies on the
        data-derived CMD. As such, it is a binned random generator
        and thus has some fundamental limitations.
        - The expected number of counts per bin is drawn ra

        There are a few limitations of this procedure:
        - Colors are drawn from the CMD of the background annulus
        - The number of stars per CMD bin is randomized according to the CMD
        - The colors/mags are then uniformly distributed within the bin
        - This leads to trouble with large bins when the cloud-in-cells 
        algorithm is applied to the simulated data
        - The positions are chosen randomly over the spherical cap of the ROI
        - Objects that are outside of the 

        WARNING: The cloud-in-cells method of generating
        the CMD leads to some difficulties since it disperses
        objects from high-density zones to low density zones.

        - Magnitudes are not randomized according to their errors

def background(self,mc_source_id=2,seed=None):
        """
        Create a simulation of the background stellar population.
        Because some stars have been clipped to generate the CMD,
        this function tends to slightly underestimate (~1%) the 
        background as compared to the true catalog.

        The simulation of background object colors relies on the
        data-derived CMD. As such, it is a binned random generator
        and thus has some fundamental limitations.
        - The expected number of counts per bin is drawn ra

        There are a few limitations of this procedure:
        - Colors are drawn from the CMD of the background annulus
        - The number of stars per CMD bin is randomized according to the CMD
        - The colors/mags are then uniformly distributed within the bin
        - This leads to trouble with large bins when the cloud-in-cells 
        algorithm is applied to the simulated data
        - The positions are chosen randomly over the spherical cap of the ROI
        - Objects that are outside of the 

        WARNING: The cloud-in-cells method of generating
        the CMD leads to some difficulties since it disperses
        objects from high-density zones to low density zones.

        - Magnitudes are not randomized according to their errors
        """
        if seed is not None: np.random.seed(seed)
        self._setup_cmd()

        # Randomize the number of stars per bin according to Poisson distribution
        nstar_per_bin = np.random.poisson(lam=self.bkg_lambda)
        nstar = nstar_per_bin.sum()

        logger.info("Simulating %i background stars..."%nstar)

        if not self.config['simulate'].get('uniform'):
            logger.info("Generating colors from background CMD.")

            # Distribute the stars within each CMD bin
            delta_color = self.bkg_centers_color[1]-self.bkg_centers_color[0]
            delta_mag   = self.bkg_centers_mag[1]-self.bkg_centers_mag[0]

            # Distribute points within each color-mag bins
            xx,yy = np.meshgrid(self.bkg_centers_color,self.bkg_centers_mag)
            color = np.repeat(xx.flatten(),repeats=nstar_per_bin.flatten())
            color += np.random.uniform(-delta_color/2.,delta_color/2.,size=nstar)
            mag_1 = np.repeat(yy.flatten(),repeats=nstar_per_bin.flatten())
            mag_1 += np.random.uniform(-delta_mag/2.,delta_mag/2.,size=nstar)
        else:
            # Uniform color-magnitude distribution
            logger.info("Generating uniform CMD.")
            mag_1 = np.random.uniform(self.config['mag']['min'],self.config['mag']['max'],size=nstar)
            color = np.random.uniform(self.config['color']['min'],self.config['color']['max'],size=nstar)

        mag_2 = mag_1 - color

        # Random points drawn from healpix subpixels
        logger.info("Generating uniform positions...")
        idx = np.random.randint(0,len(self.subpix)-1,size=nstar)
        lon,lat = pix2ang(self.nside_subpixel,self.subpix[idx])

        nside_pixel = self.nside_pixel
        pix = ang2pix(nside_pixel, lon, lat)

        # There is probably a better way to do this step without creating the full HEALPix map
        mask = -1. * np.ones(hp.nside2npix(nside_pixel))
        mask[self.roi.pixels] = self.mask.mask_1.mask_roi_sparse
        mag_lim_1 = mask[pix]
        mask = -1. * np.ones(hp.nside2npix(nside_pixel))
        mask[self.roi.pixels] = self.mask.mask_2.mask_roi_sparse
        mag_lim_2 = mask[pix]

        mag_err_1 = self.photo_err_1(mag_lim_1 - mag_1)
        mag_err_2 = self.photo_err_2(mag_lim_2 - mag_2)
        mc_source_id = mc_source_id * np.ones(len(mag_1))

        # ADW: Should magnitudes be randomized by the erros?
        #mag_1 += (np.random.normal(size=len(mag_1)) * mag_err_1)
        #mag_2 += (np.random.normal(size=len(mag_2)) * mag_err_2)

        select = (mag_lim_1>mag_1)&(mag_lim_2>mag_2)

        ### # Make sure objects lie within the original cmd (should be done later...)
        ### select &= (ugali.utils.binning.take2D(self.mask.solid_angle_cmd, color, mag_1,
        ###                                       self.roi.bins_color, self.roi.bins_mag) > 0)

        logger.info("Clipping %i simulated background stars..."%(~select).sum())
        
        hdu = ugali.observation.catalog.makeHDU(self.config,mag_1[select],mag_err_1[select],
                                                mag_2[select],mag_err_2[select],
                                                lon[select],lat[select],mc_source_id[select])
        catalog = ugali.observation.catalog.Catalog(self.config, data=hdu.data)
        return catalog