diff --git a/py/desisim/scripts/quickquasars.py b/py/desisim/scripts/quickquasars.py index 07ca9770d..52f1d3312 100644 --- a/py/desisim/scripts/quickquasars.py +++ b/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() diff --git a/py/desisim/scripts/quickspectra.py b/py/desisim/scripts/quickspectra.py index 50e74b4ea..4bb89d177 100644 --- a/py/desisim/scripts/quickspectra.py +++ b/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,