Update datamodel for templates metadata table

doc/changes.rst

@@ -7,8 +7,11 @@ desisim change log
 
 * Major refactor of newexp to add connection to upstream mocks, surveysims,
   and fiber assignment (`PR #250`_).
+* Support latest (>DR4) data model in the templates metadata table and also
+  scale simulated templates by 1e17 erg/s/cm2/Angstrom (`PR #252`_).
 
 .. _`PR #250`: https://github.com/desihub/desisim/pull/250
+.. _`PR #252`: https://github.com/desihub/desisim/pull/252
 
 0.20.0 (2017-07-12)
 -------------------

py/desisim/io.py

@@ -863,7 +863,7 @@ def write_templates(outfile, flux, wave, meta, objtype=None,
         CDELT1 = (wave[1]-wave[0], 'Wavelength step [Angstrom]'),
         LOGLAM = (0, 'linear wavelength steps, not log10'),
         AIRORVAC = ('vac', 'wavelengths in vacuum (vac) or air'),
-        BUNIT = ('erg/s/cm2/A', 'spectrum flux units')
+        BUNIT = ('1e-17 erg/(s cm2 Angstrom)', 'spectrum flux units')
         )
     hdr = fitsheader(header)
 
@@ -916,9 +916,17 @@ def empty_metatable(nmodel=1, objtype='ELG', subtype='', add_SNeIa=None):
     meta.add_column(Column(name='REDSHIFT', length=nmodel, dtype='f4',
                            data=np.zeros(nmodel)))
     meta.add_column(Column(name='MAG', length=nmodel, dtype='f4',
-                           data=np.zeros(nmodel)-1))
-    meta.add_column(Column(name='DECAM_FLUX', shape=(6,), length=nmodel, dtype='f4'))
-    meta.add_column(Column(name='WISE_FLUX', shape=(2,), length=nmodel, dtype='f4'))
+                           data=np.zeros(nmodel)-1, unit='mag'))
+    meta.add_column(Column(name='FLUX_G', length=nmodel, dtype='f4',
+                           unit='nanomaggies'))
+    meta.add_column(Column(name='FLUX_R', length=nmodel, dtype='f4',
+                           unit='nanomaggies'))
+    meta.add_column(Column(name='FLUX_Z', length=nmodel, dtype='f4',
+                           unit='nanomaggies'))
+    meta.add_column(Column(name='FLUX_W1', length=nmodel, dtype='f4',
+                           unit='nanomaggies'))
+    meta.add_column(Column(name='FLUX_W2', length=nmodel, dtype='f4',
+                           unit='nanomaggies'))
 
     meta.add_column(Column(name='OIIFLUX', length=nmodel, dtype='f4',
                            data=np.zeros(nmodel)-1, unit='erg/(s*cm2)'))

py/desisim/lya_spectra.py

@@ -95,6 +95,7 @@ def get_spectra(lyafile, nqso=None, wave=None, templateid=None, normfilter='sdss
         flux1, _, meta1 = qso.make_templates(nmodel=1, redshift=np.atleast_1d(zqso[ii]), 
                                              mag=np.atleast_1d(mag_g[ii]), seed=templateseed[ii],
                                              nocolorcuts=nocolorcuts)
+        flux1 *= 1e-17
         for col in meta1.colnames:
             meta[col][ii] = meta1[col][0]
 
@@ -114,8 +115,8 @@ def get_spectra(lyafile, nqso=None, wave=None, templateid=None, normfilter='sdss
                                                        mask_invalid=True)[normfilter])
         factor = 10**(-0.4 * mag_g[ii]) / normmaggies
         flux1 *= factor
-        meta['DECAM_FLUX'][ii] *= factor
-        meta['WISE_FLUX'][ii] *= factor
+        for key in ('FLUX_G', 'FLUX_R', 'FLUX_Z', 'FLUX_W1', 'FLUX_W2'):
+            meta[key][ii] *= factor
 
         flux[ii, :] = flux1[:]
 

py/desisim/quickcat.py

@@ -160,16 +160,26 @@ def get_redshift_efficiency(simtype, targets, truth, targets_in_tile, obsconditi
     targetid = targets['TARGETID']
     n = len(targetid)
 
+    try:
+        if 'DECAM_FLUX' in targets.colnames:
+            true_gflux = targets['DECAM_FLUX'][:, 1]
+            true_rflux = targets['DECAM_FLUX'][:, 2]
+        else:
+            true_gflux = targets['FLUX_G']
+            true_rflux = targets['FLUX_R']
+    except:
+        raise Exception('Missing photometry needed to estimate redshift efficiency!')
+
+    if (obsconditions is None) and (truth['OIIFLUX'] not in truth.colnames):
+        raise Exception('Missing obsconditions and flux information to estimate redshift efficiency')
+
     if (simtype == 'ELG'):
         # Read the model OII flux threshold (FDR fig 7.12 modified to fit redmonster efficiency on OAK)
         filename = resource_filename('desisim', 'data/quickcat_elg_oii_flux_threshold.txt')
         fdr_z, modified_fdr_oii_flux_threshold = np.loadtxt(filename, unpack=True)
 
         # Get OIIflux from truth
-        try:
-            true_oii_flux = truth['OIIFLUX']
-        except:
-            raise Exception('Missing OII flux information to estimate redshift efficiency for ELGs')
+        true_oii_flux = truth['OIIFLUX']
 
         # Compute OII flux thresholds for truez
         oii_flux_threshold = np.interp(truth['TRUEZ'],fdr_z,modified_fdr_oii_flux_threshold)
@@ -186,11 +196,6 @@ def get_redshift_efficiency(simtype, targets, truth, targets_in_tile, obsconditi
         magr, magr_eff = np.loadtxt(filename, unpack=True)
 
         # Get Rflux from truth
-        try:
-            true_rflux = targets['DECAM_FLUX'][:,2]
-        except:
-            raise Exception('Missing Rmag information to estimate redshift efficiency for LRGs')
-
         r_mag = 22.5 - 2.5*np.log10(true_rflux)
 
         mean_eff_mag=np.interp(r_mag,magr,magr_eff)
@@ -205,11 +210,8 @@ def get_redshift_efficiency(simtype, targets, truth, targets_in_tile, obsconditi
         zc, qso_gmag_threshold_vs_z = np.loadtxt(filename, unpack=True)
 
         # Get Gflux from truth
-        try:
-            true_gmag = targets['DECAM_FLUX'][:,1]
-        except:
-            raise Exception('Missing Gmag information to estimate redshift efficiency for QSOs')
-
+        true_gmag = 22.5 - 2.5 * np.log10(true_gflux) 
+
         # Computes QSO mag thresholds for truez
         qso_gmag_threshold=np.interp(truth['TRUEZ'],zc,qso_gmag_threshold_vs_z)
         assert (qso_gmag_threshold.size == true_gmag.size),"qso_gmag_threshold and true_gmag should have the same size"
@@ -233,9 +235,6 @@ def get_redshift_efficiency(simtype, targets, truth, targets_in_tile, obsconditi
         log.warning('using default redshift efficiency of {} for {}'.format(default_zeff, simtype))
         simulated_eff = default_zeff * np.ones(n)
 
-    if (obsconditions is None) and (truth['OIIFLUX'] is None) and (targets['DECAM_FLUX'] is None):
-        raise Exception('Missing obsconditions and flux information to estimate redshift efficiency')
-
     #- Get the corrections for observing conditions per tile, then
     #- correct targets on those tiles.  Parameterize in terms of failure
     #- rate instead of success rate to handle bookkeeping of targets that
@@ -343,6 +342,16 @@ def get_observed_redshifts(targets, truth, targets_in_tile, obsconditions):
     truez = truth['TRUEZ']
     targetid = truth['TARGETID']
 
+    try:
+        if 'DECAM_FLUX' in targets.colnames:
+            true_gflux = targets['DECAM_FLUX'][:, 1]
+            true_rflux = targets['DECAM_FLUX'][:, 2]
+        else:
+            true_gflux = targets['FLUX_G']
+            true_rflux = targets['FLUX_R']
+    except:
+        raise Exception('Missing photometry needed to estimate redshift efficiency!')
+
     zout = truez.copy()
     zerr = np.zeros(len(truez), dtype=np.float32)
     zwarn = np.zeros(len(truez), dtype=np.int32)
@@ -375,10 +384,8 @@ def get_observed_redshifts(targets, truth, targets_in_tile, obsconditions):
             elif (objtype == 'LRG'):
                 redbins=np.linspace(0.55,1.1,5)
                 mag = np.linspace(20.5,23.,50)
-                try:
-                    true_magr=targets['DECAM_FLUX'][ii,2]
-                except:
-                    raise Exception('Missing DECAM r flux information to estimate redshift error for LRGs')
+
+                true_magr = 22.5 - 2.5 * np.log10(true_rflux)
 
                 coefs = [[9.46882282e-05,-1.87022383e-03],[6.14601021e-05,-1.17406643e-03],\
                              [8.85342362e-05,-1.76079966e-03],[6.96202042e-05,-1.38632104e-03]]
@@ -419,10 +426,7 @@ def get_observed_redshifts(targets, truth, targets_in_tile, obsconditions):
                 redbins = np.linspace(0.5,3.5,7)
                 mag = np.linspace(21.,23.,50)
 
-                try:
-                    true_magg=targets['DECAM_FLUX'][ii,1]
-                except:
-                    raise Exception('Missing DECAM g flux information to estimate redshift error for QSOs')
+                true_magg = 22.5 - 2.5 * np.log10(true_gflux)
 
                 coefs = [[0.000156950059747,-0.00320719603886],[0.000461779391179,-0.00924485142818],\
                              [0.000458672517009,-0.0091254038977],[0.000461427968475,-0.00923812594293],\

py/desisim/scripts/quickgen.py

@@ -386,7 +386,7 @@ def main(args):
             jj.append(ii.tolist())
 
             # Sanity check on units; templates currently return ergs, not 1e-17 ergs...
-            assert (thisobj == 'SKY') or (np.max(truth['FLUX']) < 1e-6)
+            # assert (thisobj == 'SKY') or (np.max(truth['FLUX']) < 1e-6)
 
         # Sort the metadata table.
         jj = sum(jj,[])

py/desisim/targets.py

@@ -335,17 +335,15 @@ def get_targets(nspec, program, tileid=None, seed=None, specmin=0):
         else:
             raise ValueError('Unable to simulate OBJTYPE={}'.format(objtype))
 
-        flux[ii] = 1e17 * simflux
+        flux[ii] = simflux
         meta[ii] = meta1
 
-        ugrizy = 22.5-2.5*np.log10(meta1['DECAM_FLUX'].data)
-        wise = 22.5-2.5*np.log10(meta1['WISE_FLUX'].data)
         fibermap['FILTER'][ii, :6] = ['DECAM_G', 'DECAM_R', 'DECAM_Z', 'WISE_W1', 'WISE_W2']
-        fibermap['MAG'][ii, 0] = ugrizy[:, 1]
-        fibermap['MAG'][ii, 1] = ugrizy[:, 2]
-        fibermap['MAG'][ii, 2] = ugrizy[:, 4]
-        fibermap['MAG'][ii, 3] = wise[:, 0]
-        fibermap['MAG'][ii, 4] = wise[:, 1]
+        fibermap['MAG'][ii, 0] = 22.5 - 2.5 * np.log10(meta['FLUX_G'][ii])
+        fibermap['MAG'][ii, 1] = 22.5 - 2.5 * np.log10(meta['FLUX_R'][ii])
+        fibermap['MAG'][ii, 2] = 22.5 - 2.5 * np.log10(meta['FLUX_Z'][ii])
+        fibermap['MAG'][ii, 3] = 22.5 - 2.5 * np.log10(meta['FLUX_W1'][ii])
+        fibermap['MAG'][ii, 4] = 22.5 - 2.5 * np.log10(meta['FLUX_W2'][ii])
 
     #- Load fiber -> positioner mapping and tile information
     fiberpos = desimodel.io.load_fiberpos()