Merge pull request #552 from desihub/light_quickquasars

Light quickquasars

py/desisim/scripts/quickquasars.py

@@ -127,6 +127,8 @@ def parse(options=None):
 
     parser.add_argument('--save-continuum',action = "store_true", help="Save true continum to file")
 
+    parser.add_argument('--save-continuum-dwave',type=float, default=2, help="Delta wavelength to save true continum")
+
     parser.add_argument('--desi-footprint', action = "store_true" ,help="select QSOs in DESI footprint")
 
     parser.add_argument('--eboss',action = 'store_true', help='Setup footprint, number density, redshift distribution,\
@@ -143,6 +145,8 @@ def parse(options=None):
 
     parser.add_argument('--nmax', type=int, default=None, help="Max number of QSO per input file, for debugging")
 
+    parser.add_argument('--save-resolution',action='store_true', help="Save full resolution in spectra file. By default only one matrix is saved in the truth file.")
+
     if options is None:
         args = parser.parse_args()
     else:
@@ -551,7 +555,7 @@ def simulate_one_healpix(ifilename,args,model,obsconditions,decam_and_wise_filte
     tmp_qso_wave = trans_wave
 
     if args.save_continuum :
-        true_wave=np.linspace(args.wmin,args.wmax,int((args.wmax-args.wmin)/args.dwave)+1)
+        true_wave=np.linspace(args.wmin,args.wmax,int((args.wmax-args.wmin)/args.save_continuum_dwave)+1)
         true_flux=np.zeros((tmp_qso_flux.shape[0],true_wave.size))
         for q in range(tmp_qso_flux.shape[0]) :
             true_flux[q]=resample_flux(true_wave,tmp_qso_wave,tmp_qso_flux[q])
@@ -562,7 +566,7 @@ def simulate_one_healpix(ifilename,args,model,obsconditions,decam_and_wise_filte
         hdu_trueCont.name = "TRUE_CONT"
         hdu_trueCont.header['wmin'] = args.wmin
         hdu_trueCont.header['wmax'] = args.wmax
-        hdu_trueCont.header['dwave'] = args.dwave
+        hdu_trueCont.header['dwave'] = args.save_continuum_dwave
 
         del(continum_meta,true_wave,true_flux)
         log.info("True continum to be saved in {}".format(truth_filename))
@@ -717,10 +721,11 @@ def simulate_one_healpix(ifilename,args,model,obsconditions,decam_and_wise_filte
         specsim_config_file = 'desi'
 
     ### use Poisson = False to get reproducible results.
-    sim_spectra(qso_wave,qso_flux, args.program, obsconditions=obsconditions,spectra_filename=ofilename,
-        sourcetype="qso", skyerr=args.skyerr,ra=metadata["RA"],dec=metadata["DEC"],targetid=targetid,
-        meta=specmeta,seed=seed,fibermap_columns=fibermap_columns,use_poisson=False,
-        specsim_config_file=specsim_config_file, dwave_out=dwave_out)
+    ### use args.save_resolution = False to not save the matrix resolution per quasar in spectra files.
+    resolution=sim_spectra(qso_wave,qso_flux, args.program, obsconditions=obsconditions,spectra_filename=ofilename,
+                           sourcetype="qso", skyerr=args.skyerr,ra=metadata["RA"],dec=metadata["DEC"],targetid=targetid,
+                           meta=specmeta,seed=seed,fibermap_columns=fibermap_columns,use_poisson=False,
+                           specsim_config_file=specsim_config_file, dwave_out=dwave_out, save_resolution=args.save_resolution)
 
     ### Keep input redshift
     Z_spec = metadata['Z'].copy()
@@ -765,13 +770,22 @@ def simulate_one_healpix(ifilename,args,model,obsconditions,decam_and_wise_filte
     hduqso=pyfits.convenience.table_to_hdu(qsometa)
     hduqso.header['EXTNAME'] = 'TRUTH_QSO'
     hdulist=pyfits.HDUList([pyfits.PrimaryHDU(header=hdr),hdu,hduqso])
+
+
     if args.dla :
         hdulist.append(hdu_dla)
     if  args.balprob :
         hdulist.append(hdu_bal)
     if args.save_continuum :
         hdulist.append(hdu_trueCont)
-
+
+# Save one resolution matrix per camera to the truth file instead of one per quasar to the spectra files.
+    if not args.save_resolution:
+        for band in resolution.keys():
+            hdu = pyfits.ImageHDU(name="{}_RESOLUTION".format(band.upper()))
+            hdu.data = resolution[band].astype("f4")
+            hdulist.append(hdu)
+
     hdulist.writeto(truth_filename, overwrite=True)
     hdulist.close()
 

py/desisim/scripts/quickspectra.py

@@ -24,7 +24,7 @@
 
 def sim_spectra(wave, flux, program, spectra_filename, obsconditions=None,
                 sourcetype=None, targetid=None, redshift=None, expid=0, seed=0, skyerr=0.0, ra=None,
-                dec=None, meta=None, fibermap_columns=None, fullsim=False, use_poisson=True, specsim_config_file="desi", dwave_out=None):
+                dec=None, meta=None, fibermap_columns=None, fullsim=False, use_poisson=True, specsim_config_file="desi", dwave_out=None, save_resolution=True):
     """
     Simulate spectra from an input set of wavelength and flux and writes a FITS file in the Spectra format that can
     be used as input to the redshift fitter.
@@ -50,7 +50,10 @@ def sim_spectra(wave, flux, program, spectra_filename, obsconditions=None,
         meta : dictionnary, saved in primary fits header of the spectra file 
         fibermap_columns : add these columns to the fibermap
         fullsim : if True, write full simulation data in extra file per camera
-        use_poisson : if False, do not use numpy.random.poisson to simulate the Poisson noise. This is useful to get reproducible random realizations.
+        use_poisson : if False, do not use numpy.random.poisson to simulate the Poisson noise. This is useful to get reproducible random
+        realizations.
+        save_resolution : if True it will save the Resolution matrix for each spectra.
+        If False returns a resolution matrix (useful for mocks to save disk space).
     """ 
     log = get_logger()
 
@@ -193,6 +196,8 @@ def sim_spectra(wave, flux, program, spectra_filename, obsconditions=None,
     for camera in sim.instrument.cameras:
         R = Resolution(camera.get_output_resolution_matrix())
         resolution[camera.name] = np.tile(R.to_fits_array(), [nspec, 1, 1])
+        if not save_resolution :
+            resolution[camera.name] = R.to_fits_array()
 
     skyscale = skyerr * random_state.normal(size=sim.num_fibers)
 
@@ -245,10 +250,18 @@ def sim_spectra(wave, flux, program, spectra_filename, obsconditions=None,
             ivar = ivar / scale**2
             mask  = np.zeros(flux.shape).astype(int)
 
-            spec = Spectra([band], {band : wave}, {band : flux}, {band : ivar}, 
-                        resolution_data={band : resolution[band]}, 
-                        mask={band : mask}, 
-                        fibermap=spectra_fibermap, 
+            if not save_resolution :
+                spec = Spectra([band], {band : wave}, {band : flux}, {band : ivar},
+                        resolution_data=None,
+                        mask={band : mask},
+                        fibermap=spectra_fibermap,
+                        meta=meta,
+                        single=True)
+            else :
+                spec = Spectra([band], {band : wave}, {band : flux}, {band : ivar},
+                        resolution_data={band : resolution[band]},
+                        mask={band : mask},
+                        fibermap=spectra_fibermap,
                         meta=meta,
                         single=True)
 
@@ -265,6 +278,9 @@ def sim_spectra(wave, flux, program, spectra_filename, obsconditions=None,
     desisim.specsim._simulators.clear()
     desisim.specsim._simdefaults.clear()
 
+    if not save_resolution :
+        return resolution
+
 
 def parse(options=None):
     parser=argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter,