Merge f670268 into 040a4d5

bin/desi_average_flux_calibration

@@ -27,10 +27,48 @@ import yaml
 import matplotlib.pyplot as plt
 
 
+def get_ffracflux_wave(seeing, ffracflux, wave, fac_wave_power):
+    """
+    Compute the fiber acceptance for a given seeing, and normalize it at ffracflux at 6500A
+
+    Args:
+        seeing : in the r-band from the GFA [arcsec]
+        ffracflux : FIBER_FRACFLUX from GFA
+        wave : wavelength for output [A]
+        fac_wave_power : shape of the wavelength-dependence of seeing (wave ** fac_wave_power)
+        fdiam_um : physical fiber diameter [um]
+        fdiam_asec : fiber diameter [arcsec]
+
+    Returns:
+        Fiber acceptance for seeing, normalized to ffracflux at 6500A
+
+    
+    Comment:
+        GFA now provide ffrac for fiber_diameter=1.52", so no need to correct for different diameters
+        we do the normalization using fac_wave, and then to be interpolated to wave    
+    """
+    #
+    desi = desimodel.io.load_desiparams()
+    fdiam_um = desi['fibers']['diameter_um'] # AR 107 um: physical fiber diameter
+    fdiam_avg_asec = desi['fibers']['diameter_arcsec'] # AR 1.52 arcsec ; average fiber diameter
+    #
+    fac = FastFiberAcceptance("{}/data/throughput/galsim-fiber-acceptance.fits".format(os.getenv("DESIMODEL")))
+    fac_wave = np.arange(3500, 10000) # AR to do the fac calculation, and then to be interpolated to wave
+    seeing_wave = seeing * (fac_wave / 6500) ** args.fac_wave_power
+    sigma = seeing_wave / 2.35  * (fdiam_um / fdiam_avg_asec)
+    offset = np.zeros(sigma.shape)
+    ffracflux_wave = fac.value("POINT", sigma, offset) # AR fac for seeing[i] and fdiam_avg_asec
+    ffracflux_wave /= ffracflux_wave[np.abs(fac_wave-6500).argmin()] # AR normalizing to 1 at 6500A
+    ffracflux_wave *= ffracflux # AR normalising to ffracflux[i] at 6500A
+    return np.interp(wave, fac_wave, ffracflux_wave) # AR interpolating on wave
+
+
+
 if __name__ == '__main__':
 
     parser = argparse.ArgumentParser(description="Compute the average calibration for a DESI spectrograph camera using precomputed flux calibration vectors.")
 
+    parser.add_argument('--gfa_matched_coadd', type = str, default = None, required=True, help = 'full path to the GFA matched_coadd file, e.g. $DESI_ROOT/survey/GFA/offline_matched_coadd_ccds_SV1-thru_20210422.fits')
     parser.add_argument('-i','--infiles', type = str, default = None, required=True, help = 'path to ASCII file with full path to DESI frame calib fits files (1 per line)')
     parser.add_argument('-o','--outfile', type = str, default = None, required=True, help = 'output calibration file')
     parser.add_argument('--plot', action = 'store_true', help = 'plot the result')
@@ -39,6 +77,7 @@ if __name__ == '__main__':
     parser.add_argument('--corr_transp', action = 'store_true', help = 'correct for transparency')
     parser.add_argument('--seeing_key', type = str, default = 'RADPROF_FWHM_ASEC', choices = ['FWHM_ASEC', 'RADPROF_FWHM_ASEC'], help = 'key in GFA catalog for seeing (default=RADPROF_FWHM_ASEC)')
     parser.add_argument('--fac_wave_power', type = float, default = -0.25, help = 'wavelength dependence of fiber acceptance (default=-0.25)')
+    parser.add_argument('--first_night', type = int, default = None, help = 'first night for which this calibration is usable (default=None)')
 
     args = parser.parse_args()
     log=get_logger()
@@ -61,7 +100,9 @@ if __name__ == '__main__':
 
     # AR GFA information: airmass, seeing, transparency, fiber_fracflux
     # AR working only with airmass and fiber_fracflux
-    # AR using the latest GFA file (daily updated)
+    # AR now the filename is provided as input argument
+    # AR we restrict to matched_coadd, as exposures not in matched_coadd are a few
+    # AR and likely problematic
     # AR as of early Mar. 2021: GFA now reports fiber_fracflux for a 1.52" diameter fiber
     filenames = np.loadtxt(args.infiles, dtype=np.str)
     log.info("{} files in {}".format(len(filenames), args.infiles))
@@ -70,21 +111,17 @@ if __name__ == '__main__':
     seeing = np.nan + np.zeros(len(filenames))
     transparency = np.nan + np.zeros(len(filenames))
     ffracflux = np.nan + np.zeros(len(filenames))
-    # AR first checking GFA acq, then matched_coadd
-    # AR hence, retained values are matched_coadd if present, then acq if no matched_coadd
-    # AR if neither present in matched_coadd nor in acq, we discard the exposure
-    acq_fn = np.sort(glob("{}/survey/GFA/offline_acq_ccds_SV1-thru_20??????.fits".format(os.getenv("DESI_ROOT"))))[-1]
-    matched_coadd_fn = np.sort(glob("{}/survey/GFA/offline_matched_coadd_ccds_SV1-thru_20??????.fits".format(os.getenv("DESI_ROOT"))))[-1]
+    # AR matched_coadd
+    # AR we discard the exposure if not present in matched_coadd
     unq_expids = np.unique(expids)
-    for fn in [acq_fn, matched_coadd_fn]:
-        gfa = fitsio.read(fn, ext=3)
-        _, ii, gfa_ii = np.intersect1d(unq_expids, gfa["EXPID"], return_indices=True)
-        for i, gfa_i in zip(ii, gfa_ii):
-            jj = np.where(expids == unq_expids[i])[0]
-            airmass[jj] = gfa["AIRMASS"][gfa_i]
-            seeing[jj] = gfa[args.seeing_key][gfa_i]
-            transparency[jj] = gfa["TRANSPARENCY"][gfa_i]
-            ffracflux[jj] = gfa["FIBER_FRACFLUX"][gfa_i]
+    gfa = fitsio.read(args.gfa_matched_coadd, ext=3)
+    _, ii, gfa_ii = np.intersect1d(unq_expids, gfa["EXPID"], return_indices=True)
+    for i, gfa_i in zip(ii, gfa_ii):
+        jj = np.where(expids == unq_expids[i])[0]
+        airmass[jj] = gfa["AIRMASS"][gfa_i]
+        seeing[jj] = gfa[args.seeing_key][gfa_i]
+        transparency[jj] = gfa["TRANSPARENCY"][gfa_i]
+        ffracflux[jj] = gfa["FIBER_FRACFLUX"][gfa_i]
     # AR cutting the input file list on exposures having - meaningful - GFA measurements
     keep = (airmass >= 1.0) & (seeing >= 0) & (ffracflux >= 0)
     if args.corr_transp:
@@ -93,13 +130,6 @@ if __name__ == '__main__':
     log.info("discarding the {} following exposures (no GFA): {}".format((~keep).sum(), ",".join(filenames[~keep])))
     filenames = filenames[keep]
 
-    # AR fiber acceptance, fiber diameters
-    fac = FastFiberAcceptance("{}/data/throughput/galsim-fiber-acceptance.fits".format(os.getenv("DESIMODEL")))
-    fac_wave = np.arange(3500, 10000) # AR to do the fac calculation, and then to be interpolated to wave
-    desi = desimodel.io.load_desiparams()
-    fdiam_um = desi['fibers']['diameter_um'] # AR 107 um: physical fiber diameter
-    fdiam_avg_asec = desi['fibers']['diameter_arcsec'] # AR 1.52 arcsec ; average fiber diameter
-
     # AR storing indexes kept in the loop
     ii = []
     # AR looping on filenames
@@ -148,15 +178,9 @@ if __name__ == '__main__':
 
         # AR GFA now provide ffrac for fiber_diameter=1.52", so no need to correct for different diameters
         # AR we do the normalization using fac_wave
-        seeing_wave = seeing[i] * (fac_wave / 6500) ** args.fac_wave_power
-        sigma = seeing_wave / 2.35  * (fdiam_um / fdiam_avg_asec)
-        offset = np.zeros(sigma.shape)
-        ffracflux_wave = fac.value("POINT", sigma, offset) # AR fac for seeing[i] and fdiam_avg_asec
-        ffracflux_wave /= ffracflux_wave[np.abs(fac_wave-6500).argmin()] # AR normalizing to 1 at 6500A
-        ffracflux_wave *= ffracflux[i] # AR normalising to ffracflux[i] at 6500A
-        ffracflux_wave_interp = np.interp(wave, fac_wave, ffracflux_wave) # AR interpolating on wave
+        ffracflux_wave = get_ffracflux_wave(seeing[i], ffracflux[i], wave, args.fac_wave_power)
         # AR normalizing
-        norm = exptime * ffracflux_wave_interp
+        norm = exptime * ffracflux_wave
         text = "applying calibs.append(mcalib / exptime / ffracflux_wave_interp"
         if args.corr_transp:
             norm *= transparency[i]
@@ -167,10 +191,19 @@ if __name__ == '__main__':
         #
         ii += [i]
 
+    # AR cutting airmass, seeing, ffracflux on kept exposures
+    airmass, seeing, ffracflux = airmass[ii], seeing[ii], ffracflux[ii]
+
+    # AR we multiply by the ffracflux_wave from median conditions
+    # AR to have a curve including the fiber acceptance
+    # AR in the initial spirit of this script
+    # AR not sure one needs to normalize by transparency if args.corr_transp; not done here
+    median_ffracflux_wave = get_ffracflux_wave(np.median(seeing), np.median(ffracflux), wave, args.fac_wave_power)
+    log.info("multiplying by the ffracflux_wave from median conditions")
+    calibs = [calib * median_ffracflux_wave for calib in calibs]
+    #
     calibs=np.array(calibs)
     nexp=calibs.shape[0]
-    # AR cutting airmass, seeing on kept exposures
-    airmass, seeing = airmass[ii], seeing[ii]
 
     # compute an average calibration vector
     ###########################################################
@@ -345,7 +378,12 @@ if __name__ == '__main__':
                              pivot_airmass=pivot_airmass,
                              pivot_seeing=pivot_seeing,
                              atmospheric_extinction_uncertainty=airmass_term_uncertainty,
-                             seeing_term_uncertainty=seeing_term_uncertainty)
+                             seeing_term_uncertainty=seeing_term_uncertainty,
+                             median_seeing=np.median(seeing),
+                             median_ffracflux=np.median(ffracflux),
+                             fac_wave_power=args.fac_wave_power,
+                             ffracflux_wave=median_ffracflux_wave,
+                             first_night=args.first_night)
 
     # write the result ...
 

bin/desi_compute_skymags

@@ -19,6 +19,8 @@ from astropy.table import Table
 from desiutil.log import get_logger
 from desispec.io import specprod_root
 from desispec.skymag import compute_skymag
+from desispec.util import parse_int_args
+
 
 
 def parse(options=None):
@@ -33,7 +35,7 @@ def parse(options=None):
     parser.add_argument('-e','--expids', type = str, default = None, required=False,
                         help = 'Comma separated list of exp ids to process')
     parser.add_argument('-n','--nights', type = str, default = None, required=False,
-                        help = 'Comma separated list of nights to process')
+                        help = 'Comma, or colon separated list of nights to process. ex: 20210501,20210502 or 20210501:20210531')
     parser.add_argument('--nproc', type = int, default = 1,
                         help = 'Multiprocessing.')
 
@@ -93,7 +95,7 @@ def main():
             except ValueError as e :
                 log.warning("ignore {}".format(dirname))
     else :
-        nights=[int(val) for val in args.nights.split(",")]
+        nights=parse_int_args(args.nights)
 
     log.info("nights = {}".format(nights))
     if expids is not None : log.info('expids = {}'.format(expids))

py/desispec/averagefluxcalibration.py

@@ -1,11 +1,14 @@
 import numpy as np
 import scipy.constants
 from desimodel.io import load_desiparams
+from desiutil.log import get_logger
 
 class AverageFluxCalib(object):
     def __init__(self, wave, average_calib, atmospheric_extinction, seeing_term, \
                  pivot_airmass, pivot_seeing,\
-                 atmospheric_extinction_uncertainty = None, seeing_term_uncertainty = None):
+                 atmospheric_extinction_uncertainty = None, seeing_term_uncertainty = None,\
+                 median_seeing = None, median_ffracflux = None, fac_wave_power = None, ffracflux_wave = None,
+                 first_night = None):
         """Lightweight wrapper object for average flux calibration data
 
         Args:
@@ -17,6 +20,11 @@ def __init__(self, wave, average_calib, atmospheric_extinction, seeing_term, \
             pivot_seeing : float, seeing value for average_calib (same definition and unit as SEEING keyword in images)
             atmospheric_extinction_uncertainty : 1D[nwave] uncertainty on extinction term, magnitude
             seeing_term_uncertainty : 1D[nwave], uncertainty on seeing term magnitude
+            median_seeing : float, median seeing [arcsec]
+            median_ffracflux : float, median GFA FIBER_FRACFLUX
+            fac_wave_power : float, wavelength-dependence of the fiber acceptance (wave ** fac_wave_power)
+            ffracflux_wave : 1D[nwave], fiber acceptance for {median_seeing, median_ffrac} at 6500A, following wave ** fac_wave_power
+            first_night : first night for which this calibration is usable
 
         All arguments become attributes,
         the calib vector should be in units of [electrons]/[1e-17 erg/cm^2].
@@ -38,6 +46,11 @@ def __init__(self, wave, average_calib, atmospheric_extinction, seeing_term, \
         self.pivot_seeing = pivot_seeing
         self.atmospheric_extinction_uncertainty = atmospheric_extinction_uncertainty
         self.seeing_term_uncertainty = seeing_term_uncertainty
+        self.median_seeing = median_seeing
+        self.median_ffracflux = median_ffracflux
+        self.fac_wave_power = fac_wave_power
+        self.ffracflux_wave = ffracflux_wave
+        self.first_night = first_night
         self.meta = dict(units='electrons/(1e-17 erg/cm^2)')
         self.desiparams = load_desiparams()
 
@@ -80,4 +93,14 @@ def throughput(self,airmass=1,seeing=1.1) :
         cspeed *= 1e10 # A/s
         energy = hplanck*cspeed/self.wave # (erg/photon)
         cal_scale *= energy # scale*value in (electron/photon)
+        # AR dividing by ffracflux_wave if present
+        # AR (if present, it means the average_calib has been multiplied by it)
+        if self.ffracflux_wave is not None:
+            cal_scale /= self.ffracflux_wave
+        else :
+            default_fiber_fracflux = 0.6
+            log = get_logger()
+            log.warning("Considering a default fiber frac flux of {} to convert calib to throughput".format(default_fiber_fracflux))
+            cal_scale /= default_fiber_fracflux
+
         return cal_scale * self.value(airmass,seeing)

py/desispec/efftime.py

@@ -64,8 +64,8 @@ def compute_efftime(table,
     fiber_area_arcsec2 = np.pi*(fiber_diameter_arcsec/2)**2
 
     # flux in fiber artificially divided by fiber_area_arcsec2  because the sky flux is per arcsec2
-    fflux_bright = flux_bright_nom * fiber_fracflux_bgs / airfac / ebvfac / fiber_area_arcsec2
-    fflux_backup = flux_backup_nom * fiber_fracflux_psf / airfac / ebvfac / fiber_area_arcsec2
+    fflux_bright = flux_bright_nom * transparency * fiber_fracflux_bgs / airfac / ebvfac / fiber_area_arcsec2
+    fflux_backup = flux_backup_nom * transparency * fiber_fracflux_psf / airfac / ebvfac / fiber_area_arcsec2
 
     # AR effective sky
     effsky_dark = (sky + sky_rdn * exptime_nom / exptime) / (1.0 + sky_rdn / sky_nom)

py/desispec/io/fluxcalibration.py

@@ -212,7 +212,14 @@ def write_average_flux_calibration(outfile, averagefluxcalib):
       hx.append( fits.ImageHDU(averagefluxcalib.atmospheric_extinction_uncertainty.astype('f4'), name='ATERM_ERR') )
     if averagefluxcalib.seeing_term_uncertainty is not None :
         hx.append( fits.ImageHDU(averagefluxcalib.seeing_term_uncertainty.astype('f4'), name='STERM_ERR') )
-
+    # AR wave-dependent FIBER_FRACFLUX curve for the median seeing and median FIBER_FRACFLUX of the used exposures
+    if averagefluxcalib.ffracflux_wave is not None:
+        hx.append( fits.ImageHDU(averagefluxcalib.ffracflux_wave.astype('f4'), name='FFRACFLUX_WAVE') )
+        hx[-1].header['MDSEEING'] = averagefluxcalib.median_seeing
+        hx[-1].header['MDFFRACF'] = averagefluxcalib.median_ffracflux
+        hx[-1].header['FACWPOW'] = averagefluxcalib.fac_wave_power
+        hx[-1].header['FSTNIGHT'] = averagefluxcalib.first_night
+
     t0 = time.time()
     hx.writeto(outfile+'.tmp', overwrite=True, checksum=True)
     os.rename(outfile+'.tmp', outfile)
@@ -245,6 +252,16 @@ def read_average_flux_calibration(filename):
             seeing_term_uncertainty = native_endian(fx["STERM_ERR"].data.astype('f8'))
         else :
             seeing_term_uncertainty = None
+        if "FFRACFLUX_WAVE" in fx:
+            median_seeing          = fx["FFRACFLUX_WAVE"].header["MDSEEING"]
+            median_ffracflux       = fx["FFRACFLUX_WAVE"].header["MDFFRACF"]
+            fac_wave_power         = fx["FFRACFLUX_WAVE"].header["FACWPOW"]
+            ffracflux_wave         = native_endian(fx["FFRACFLUX_WAVE"].data.astype('f8'))
+            first_night            = fx["FFRACFLUX_WAVE"].header["FSTNIGHT"]
+        else:
+            median_seeing, median_ffracflux, fac_wave_power, ffracflux_wave = None, None, None, None
+            first_night = None
+
 
     duration = time.time() - t0
     log.info(iotime.format('read', filename, duration))
@@ -256,7 +273,12 @@ def read_average_flux_calibration(filename):
                                   pivot_airmass=pivot_airmass,
                                   pivot_seeing=pivot_seeing,
                                   atmospheric_extinction_uncertainty=atmospheric_extinction_uncertainty,
-                                  seeing_term_uncertainty=seeing_term_uncertainty)
+                                  seeing_term_uncertainty=seeing_term_uncertainty,
+                                  median_seeing=median_seeing,
+                                  median_ffracflux=median_ffracflux,
+                                  fac_wave_power=fac_wave_power,
+                                  ffracflux_wave=ffracflux_wave,
+                                  first_night=first_night)
 
     return afluxcalib
 

py/desispec/skymag.py

@@ -96,20 +96,30 @@ def compute_skymag(night, expid, specprod_dir=None):
             header=fitsio.read_header(filename)
             exptime=header["EXPTIME"]
 
-            # for now we use a fixed calibration as used in DESI-6043 for which we know what was the fiber aperture loss
-            cal_filename="{}/spec/fluxcalib/fluxcalibnight-{}-20201216.fits".format(os.environ["DESI_SPECTRO_CALIB"],camera)
-            # apply the correction from
-            fiber_acceptance_for_point_sources = 0.60 # see DESI-6043
-            mean_fiber_diameter_arcsec = 1.52 # see DESI-6043
-            fiber_area_arcsec = np.pi*(mean_fiber_diameter_arcsec/2)**2
+            # use fixed calibrations
+            if night < 20210318 : # before mirror cleaning
+                cal_filename="{}/spec/fluxcalib/fluxcalibaverage-{}-20201214.fits".format(os.environ["DESI_SPECTRO_CALIB"],camera)
+            else :
+                cal_filename="{}/spec/fluxcalib/fluxcalibaverage-{}-20210318.fits".format(os.environ["DESI_SPECTRO_CALIB"],camera)
 
             acal = _get_average_calibration(cal_filename)
-
             begin, end = istitch[camera]
             flux = np.interp(fullwave[begin:end], skywave, skyflux)
-            acal_val = np.interp(fullwave[begin:end], acal.wave, acal.value())
+            acal_val = acal.value()
+
+            if acal.ffracflux_wave is not None :
+                acal_val /= acal.ffracflux_wave
+            else :
+                default_ffracflux = 0.6 # see DESI-6043
+                log.warning("use a default fiber acceptance correction = {}".format(default_ffracflux))
+                acal_val /= default_ffracflux
+
+            acal_val = np.interp(fullwave[begin:end], acal.wave, acal_val)
+
+            mean_fiber_diameter_arcsec = 1.52 # see DESI-6043
+            fiber_area_arcsec = np.pi*(mean_fiber_diameter_arcsec/2)**2
 
-            sky[begin:end] = flux / exptime / acal_val * fiber_acceptance_for_point_sources / fiber_area_arcsec * 1e-17 # ergs/s/cm2/A/arcsec2
+            sky[begin:end] = flux / exptime / acal_val / fiber_area_arcsec * 1e-17 # ergs/s/cm2/A/arcsec2
 
         if not ok : continue # to next spectrograph
         sky_spectra.append(sky)