Add airmass dependent Bandpasses

config/imsim-config.yaml

@@ -133,7 +133,14 @@ image:
 
     det_name: "$det_name"
 
-    bandpass: { type: RubinBandpass, band: "$band" }
+    bandpass:
+        type: RubinBandpass
+        band: $band
+        # You can optionally add an airmass here to incorporate into the bandpass.  If you leave
+        # this out, then the fiducial bandpass with airmass = 1.2 will be used.
+        airmass: { type: OpsimData, field: airmass }
+        camera: "@output.camera"
+        det_name: $det_name
 
     wcs:
         type: Batoid

imsim/bandpass.py

@@ -1,12 +1,71 @@
-
+import numpy as np
 import os
+from pathlib import Path
 import galsim
 from galsim.config import RegisterBandpassType
+from astropy.table import Table
 
 __all__ = ['RubinBandpass']
 
 
-def RubinBandpass(band, logger=None):
+class AtmInterpolator:
+    """ Interpolate atmospheric transmission curves from throughputs repo.
+    Linear interpolation of log(transmission) independently for every wavelength
+    does a good job.  Extrapolation is done by assuming a constant slope in
+    log(transmission).
+
+    Parameters
+    ----------
+    Xs : `np.array`
+        Airmass values at which the transmission curves are tabulated.
+    arr : `np.array`
+        Transmission curves at the airmass values.  First index is the airmass,
+        second index is the wavelength.
+    """
+    def __init__(self, Xs, arr):
+        self.Xs = Xs
+        self.arr = arr
+        with np.errstate(all='ignore'):
+            self.logarr = np.log(arr)
+            self.log_extrapolation_slope = (
+                (self.logarr[-1] - self.logarr[-2])/(self.Xs[-1]-self.Xs[-2])
+            )
+
+    def __call__(self, X):
+        """ Evaluate atmospheric transmission curve at airmass X.
+
+        Parameters
+        ----------
+        X : `float`
+            Airmass at which to evaluate the transmission curve.
+
+        Returns
+        -------
+        out : `np.array`
+            Transmission curve at the requested airmass.
+        """
+        assert X >= 1.0
+        idx = np.searchsorted(self.Xs, X, side='right')
+        if idx == len(self.Xs):  # extrapolate
+            frac = (X - self.Xs[idx-1])
+            out = np.array(self.logarr[idx-1])
+            out += frac*self.log_extrapolation_slope
+        else:
+            frac = (X - self.Xs[idx-1]) / (self.Xs[idx] - self.Xs[idx-1])
+            out = (1-frac)*np.array(self.logarr[idx-1])
+            out += frac*self.logarr[idx]
+        out = np.exp(out)
+        out[~np.isfinite(out)] = 0.0  # Clean up exp(log(0)) => 0
+        return out
+
+
+def RubinBandpass(
+    band,
+    airmass=None,
+    camera=None,
+    det_name=None,
+    logger=None
+):
     """Return one of the Rubin bandpasses, specified by the single-letter name.
 
     The zeropoint is automatically set to the AB zeropoint normalization.
@@ -15,24 +74,118 @@ def RubinBandpass(band, logger=None):
     ----------
     band : `str`
         The name of the bandpass.  One of u,g,r,i,z,y.
+    airmass : `float`, optional
+        The airmass at which to evaluate the bandpass.  If None, use the
+        standard X=1.2 bandpass.
+    camera : `str`, optional
+        Name of the camera for which to incorporate the detector-specific
+        quantum efficiency.  If None, use the standard hardware bandpass.
+    det_name : `str`, optional
+        Name of the detector for which to incorporate the quantum efficiency.  If
+        None, use the standard hardware bandpass.
     logger : logging.Logger
         If provided, a logger for logging debug statements.
     """
-    # This uses the baseline throughput files from lsst.sims
-    sims_dir = os.getenv("RUBIN_SIM_DATA_DIR")
-    file_name = os.path.join(sims_dir, "throughputs", "baseline", f"total_{band}.dat")
-    if not os.path.isfile(file_name):
-        # If the user doesn't have the RUBIN_SIM_DATA_DIR defined, or if they don't have
-        # the correct files installed, revert to the GalSim files.
-        logger = galsim.config.LoggerWrapper(logger)
-        logger.warning("Warning: Using the old bandpass files from GalSim, not lsst.sims")
-        file_name = f"LSST_{band}.dat"
-    bp = galsim.Bandpass(file_name, wave_type='nm')
+    if (camera is None) != (det_name is None):
+        raise ValueError("Must provide both camera and det_name if using one.")
+    match camera:
+        case 'LsstCam':
+            camera = 'lsstCam'
+        case 'LsstComCam':
+            camera = 'comCamSim'
+        case _:
+            camera = camera
+    tp_path = Path(os.getenv("RUBIN_SIM_DATA_DIR")) / "throughputs"
+    if airmass is None and camera is None:
+        file_name = tp_path / "baseline" / f"total_{band}.dat"
+        if not file_name.is_file():
+            logger = galsim.config.LoggerWrapper(logger)
+            logger.warning("Warning: Using the old bandpass files from GalSim, not rubin_sim")
+            file_name = f"LSST_{band}.dat"
+        bp = galsim.Bandpass(str(file_name), wave_type='nm')
+    else:
+        if airmass is None:
+            airmass = 1.2
+        # Could be more efficient by only reading in the bracketing airmasses,
+        # but probably doesn't matter much here.
+        atmos = {}
+        for f in (tp_path / "atmos").glob("atmos_??_aerosol.dat"):
+            X = float(f.name[-14:-12])/10.0
+            wave, tput = np.genfromtxt(f).T
+            atmos[X] = wave, tput
+        Xs = np.array(sorted(atmos.keys()))
+        arr = np.array([atmos[X][1] for X in Xs])
+
+        interpolator = AtmInterpolator(Xs, arr)
+        tput = interpolator(airmass)
+        bp_atm = galsim.Bandpass(
+            galsim.LookupTable(
+                wave, tput, interpolant='linear'
+            ),
+            wave_type='nm'
+        )
+
+        if camera is not None:
+            try:
+                old_path = Path(os.getenv("OBS_LSST_DATA_DIR"))
+            except TypeError:
+                raise ValueError(
+                    "Unable to find OBS_LSST_DATA; required if using camera or det_name for bandpass."
+                )
+            old_path = old_path / camera / "transmission_sensor" / det_name.lower()
+            det_files = list(old_path.glob("*.ecsv"))
+            if len(det_files) != 1:
+                raise ValueError(f"Expected 1 detector file, found {len(det_files)}")
+            det_file = det_files[0]
+            table = Table.read(det_file)
+            # Average over amplifiers
+            amps = np.unique(table['amp_name'])
+            det_wave = table[table['amp_name'] == amps[0]]['wavelength']
+            det_tput = table[table['amp_name'] == amps[0]]['efficiency']/100
+            for amp in amps[1:]:
+                assert np.all(det_wave == table[table['amp_name'] == amp]['wavelength'])
+                det_tput += table[table['amp_name'] == amp]['efficiency']/100
+            det_tput /= len(amps)
+            bp_det = galsim.Bandpass(
+                galsim.LookupTable(
+                    det_wave, det_tput, interpolant='linear'
+                ),
+                wave_type='nm'
+            )
+
+            # Get the rest of the hardware throughput
+            optics_wave, optics_tput = np.genfromtxt(
+                tp_path / "baseline" / f"filter_{band}.dat"
+            ).T
+            for f in ["m1.dat", "m2.dat", "m3.dat", "lens1.dat", "lens2.dat", "lens3.dat"]:
+                wave, tput = np.genfromtxt(
+                    tp_path / "baseline" / f
+                ).T
+                assert np.all(wave == optics_wave)
+                optics_tput *= tput
+            bp_optics = galsim.Bandpass(
+                galsim.LookupTable(
+                    optics_wave, optics_tput, interpolant='linear'
+                ),
+                wave_type='nm'
+            )
+            bp_hardware = bp_det * bp_optics
+        else:
+            file_name = tp_path / "baseline" / f"hardware_{band}.dat"
+            wave, hardware_tput = np.genfromtxt(file_name).T
+            bp_hardware = galsim.Bandpass(
+                galsim.LookupTable(
+                    wave, hardware_tput, interpolant='linear'
+                ),
+                wave_type='nm'
+            )
+        bp = bp_atm * bp_hardware
     bp = bp.truncate(relative_throughput=1.e-3)
     bp = bp.thin()
     bp = bp.withZeropoint('AB')
     return bp
 
+
 class RubinBandpassBuilder(galsim.config.BandpassBuilder):
     """A class for building a RubinBandpass in the config file
     """
@@ -48,9 +201,19 @@ def buildBandpass(self, config, base, logger):
             the constructed Bandpass object.
         """
         req = { 'band' : str }
-        kwargs, safe = galsim.config.GetAllParams(config, base, req=req)
+        opt = {
+            'airmass' : float,
+            'camera' : str,
+            'det_name' : str
+        }
+        kwargs, safe = galsim.config.GetAllParams(
+            config, base, req=req, opt=opt
+        )
         kwargs['logger'] = logger
         bp = RubinBandpass(**kwargs)
+
+        # Also, store the kwargs=None version in the base config.
+        base['fiducial_bandpass'] = RubinBandpass(band=kwargs['band'], logger=logger)
         logger.debug('bandpass = %s', bp)
         return bp, safe
 

imsim/photon_ops.py

@@ -3,10 +3,10 @@
 from functools import lru_cache
 import numpy as np
 
-import galsim
 import batoid
-from galsim import PhotonArray, PhotonOp, GaussianDeviate
+from galsim import Bandpass, PhotonArray, PhotonOp, GaussianDeviate
 from galsim.config import RegisterPhotonOpType, PhotonOpBuilder, GetAllParams
+from galsim.config import BuildBandpass, GalSimConfigError
 from galsim.celestial import CelestialCoord
 from galsim.config.util import get_cls_params
 from coord import Angle
@@ -481,3 +481,36 @@ def ray_vector_to_photon_array(
     out.dxdz, out.dydz = (out.dxdz.ravel(), out.dydz.ravel())
     out.flux[ray_vector.vignetted] = 0.0
     return out
+
+
+class BandpassRatio(PhotonOp):
+    """Photon operator that reweights photon fluxes to effect
+    a specified bandpass from photons initially sampled from
+    a different bandpass.
+    """
+    def __init__(
+        self,
+        target_bandpass: Bandpass,
+        initial_bandpass: Bandpass
+    ):
+        self.target = target_bandpass
+        self.initial = initial_bandpass
+        self.ratio = self.target / self.initial
+
+    def applyTo(self, photon_array, local_wcs=None, rng=None):
+        photon_array.flux *= self.ratio(photon_array.wavelength)
+
+
+class BandpassRatioBuilder(PhotonOpBuilder):
+    def buildPhotonOp(self, config, base, logger):
+        if 'target_bandpass' not in config:
+            raise GalSimConfigError("target_bandpass is required for BandpassRatio")
+        if 'initial_bandpass' not in config:
+            raise GalSimConfigError("initial_bandpass is required for BandpassRatio")
+        kwargs = {}
+        kwargs['target_bandpass'] = BuildBandpass(config, 'target_bandpass', base, logger)[0]
+        kwargs['initial_bandpass'] = BuildBandpass(config, 'initial_bandpass', base, logger)[0]
+        return BandpassRatio(**kwargs)
+
+
+RegisterPhotonOpType('BandpassRatio', BandpassRatioBuilder())

imsim/stamp.py

@@ -8,7 +8,7 @@
 
 from .diffraction_fft import apply_diffraction_psf
 from .camera import get_camera
-
+from .photon_ops import BandpassRatio
 
 @dataclass
 class DiffractionFFT:
@@ -119,10 +119,20 @@ def setup(self, config, base, xsize, ysize, ignore, logger):
         camera = get_camera(params.get('camera', 'LsstCam'))
         if self.vignetting:
             self.det = camera[params['det_name']]
-        if not hasattr(gal, 'flux'):
-            # In this case, the object flux has not been precomputed
-            # or cached by the skyCatalogs code.
+        if hasattr(gal, 'flux'):
+            # In this case, the object flux has been precomputed.  If our
+            # realized bandpass is different from the fiducial bandpass used
+            # in the precomputation, then we'll need to reweight the flux.
+            # We'll only do this if the bandpass was obtained from
+            # RubinBandpass though.
+            self.fiducial_bandpass = base.get('fiducial_bandpass', None)
+            self.do_reweight = (
+                self.fiducial_bandpass is not None
+                and self.fiducial_bandpass != bandpass
+            )
+        else:
             gal.flux = gal.calculateFlux(bandpass)
+            self.do_reweight = False
         self.nominal_flux = gal.flux
 
         # For photon shooting rendering, precompute the realization of the Poisson variate.
@@ -636,6 +646,12 @@ def draw(self, prof, image, method, offset, config, base, logger):
         max_flux_simple = config.get('max_flux_simple', 100)
         faint = self.nominal_flux < max_flux_simple
         bandpass = base['bandpass']
+
+        if self.do_reweight:
+            initial_flux_bandpass = self.fiducial_bandpass
+        else:
+            initial_flux_bandpass = base['bandpass']
+
         if faint:
             logger.info("Flux = %.0f  Using trivial sed", self.nominal_flux)
             # cosmoDC2 galaxies with z > 2.71 and some SNANA objects
@@ -650,7 +666,7 @@ def draw(self, prof, image, method, offset, config, base, logger):
                 if sed_value != 0:
                     break
             if sed_value == 0:
-                # We can't evalue the profile for this object, so skip it.
+                # We can't evaluate the profile for this object, so skip it.
                 obj_num = base.get('obj_num')
                 object_id = base.get('object_id')
                 logger.warning("Zero-valued SED for faint object %d, "
@@ -675,10 +691,8 @@ def draw(self, prof, image, method, offset, config, base, logger):
             maxN = galsim.config.ParseValue(config, 'maxN', base, int)[0]
 
         if method == 'fft':
-            # For FFT, we may have adjusted the flux to account for vignetting.
-            # So update the flux to self.fft_flux if it's different.
             if self.fft_flux != self.nominal_flux:
-                gal = gal.withFlux(self.fft_flux, bandpass)
+                gal = gal.withFlux(self.fft_flux, initial_flux_bandpass)
 
             fft_image = image.copy()
             fft_offset = offset
@@ -724,12 +738,26 @@ def draw(self, prof, image, method, offset, config, base, logger):
         else:
             # For photon shooting, use the poisson-realization of the flux
             # and tell GalSim not to redo the Poisson realization.
-            gal = gal.withFlux(self.phot_flux, bandpass)
+            # Use the initial_flux_bandpass here, and use a photon op to get to the realized
+            # bandpass below.
+            gal = gal.withFlux(self.phot_flux, initial_flux_bandpass)
 
             if not faint and 'photon_ops' in config:
                 photon_ops = galsim.config.BuildPhotonOps(config, 'photon_ops', base, logger)
             else:
                 photon_ops = []
+
+            if self.do_reweight:
+                photon_ops.append(
+                    BandpassRatio(
+                        target_bandpass=bandpass,
+                        initial_bandpass=self.fiducial_bandpass,
+                    )
+                )
+                bp_for_drawImage = self.fiducial_bandpass
+            else:
+                bp_for_drawImage = bandpass
+
             # Put the psfs at the start of the photon_ops.
             # Probably a little better to put them a bit later than the start in some cases
             # (e.g. after TimeSampler, PupilAnnulusSampler), but leave that as a todo for now.
@@ -742,7 +770,7 @@ def draw(self, prof, image, method, offset, config, base, logger):
                 if sensor is not None:
                     sensor.updateRNG(self.rng)
 
-            image = gal.drawImage(bandpass,
+            image = gal.drawImage(bp_for_drawImage,
                                   method='phot',
                                   offset=offset,
                                   rng=self.rng,