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