Use hardware bandpass for sky level

config/imsim-config-instcat.yaml

@@ -24,6 +24,7 @@ input.instance_catalog:
     # This should be overridden below or on the command line.
     file_name: default_catalog_file.txt
     sed_dir: $os.environ.get('SIMS_SED_LIBRARY_DIR')
+    pupil_area: $pupil_area
 
 input.opsim_data:
     # Read the visit meta data. By default, we use the same file as the above

config/imsim-config-skycat.yaml

@@ -25,6 +25,7 @@ input.sky_catalog:
     file_name: default_sky_cat.yaml
     band: $band
     mjd: { type: OpsimData, field: mjd }
+    pupil_area: $pupil_area
 
 input.opsim_data:
     file_name: default_opsim.db

config/imsim-config.yaml

@@ -39,6 +39,12 @@ eval_variables:
         type: Degrees
         theta: { type: OpsimData, field: rotTelPos }
 
+    fpupil_R_outer: 4.18
+    fpupil_R_inner: 2.55  # M1 inner diameter is 2.558, but we need a bit of slack for off-axis rays
+    fpupil_area:
+        type: Eval
+        str: "np.pi * (pupil_R_outer**2 - pupil_R_inner**2) * 100**2"
+
     sband: { type: OpsimData, field: band }
 
     fexptime: { type: OpsimData, field: exptime }
@@ -75,6 +81,7 @@ input:
         # background level.
         exptime: $exptime
         mjd: { type: OpsimData, field: mjd }
+        pupil_area: $pupil_area
 
     atm_psf:
         # This enables the AtmosphericPSF type for the PSF
@@ -260,8 +267,8 @@ stamp:
             exptime: $exptime
         -
             type: PupilAnnulusSampler
-            R_outer: 4.18
-            R_inner: 2.55  # M1 inner diameter is 2.558, but we need a bit of slack for off-axis rays
+            R_outer: "$pupil_R_outer"
+            R_inner: "$pupil_R_inner"
         -
             type: PhotonDCR
             base_wavelength: $bandpass.effective_wavelength

imsim/bandpass.py

@@ -96,6 +96,7 @@ def RubinBandpass(
         case _:
             camera = camera
     tp_path = Path(os.getenv("RUBIN_SIM_DATA_DIR")) / "throughputs"
+    bp_hardware = None
     if airmass is None and camera is None:
         file_name = tp_path / "baseline" / f"total_{band}.dat"
         if not file_name.is_file():
@@ -183,6 +184,12 @@ def RubinBandpass(
     bp = bp.truncate(relative_throughput=1.e-3)
     bp = bp.thin()
     bp = bp.withZeropoint('AB')
+    if bp_hardware is not None:
+        bp_hardware.truncate(relative_throughput=1.e-3)
+        bp_hardware.thin()
+        bp_hardware.withZeropoint('AB')
+        bp.bp_hardware = bp_hardware
+
     return bp
 
 

imsim/instcat.py

@@ -12,6 +12,7 @@
 from galsim import CelestialCoord
 import galsim
 import pickle
+from .utils import RUBIN_AREA
 
 
 def clarify_radec_limits(
@@ -171,11 +172,13 @@ class InstCatalog(object):
 
     def __init__(self, file_name, wcs, xsize=4096, ysize=4096, sed_dir=None,
                  edge_pix=100, sort_mag=True, flip_g2=True, approx_nobjects=None,
-                 min_source=None, skip_invalid=True, logger=None):
+                 pupil_area=RUBIN_AREA, min_source=None, skip_invalid=True,
+                 logger=None):
         logger = galsim.config.LoggerWrapper(logger)
         self.file_name = file_name
         self.flip_g2 = flip_g2
         self.approx_nobjects = approx_nobjects
+        self.pupil_area = pupil_area
         self._sed_cache = {}
 
         if sed_dir is None:
@@ -516,7 +519,7 @@ def getObj(self, index, gsparams=None, rng=None, exptime=30):
 
         # This gives the normalization in photons/cm^2/sec.
         # Multiply by area and exptime to get photons.
-        fAt = flux * self._rubin_area * exptime
+        fAt = flux * self.pupil_area * exptime
 
         sed = self.getSED(index)
         return obj.withFlux(fAt) * sed
@@ -599,6 +602,7 @@ def getKwargs(self, config, base, logger):
                 'sort_mag' : bool,
                 'flip_g2' : bool,
                 'approx_nobjects' : int,
+                'pupil_area' : float,
                 'min_source' : int,
                 'skip_invalid' : bool,
               }

imsim/sky_model.py

@@ -3,21 +3,19 @@
 import warnings
 import numpy as np
 import galsim
-import pickle
 from galsim.config import InputLoader, RegisterInputType, RegisterValueType
 from scipy.interpolate import RegularGridInterpolator, RectBivariateSpline
 import os
 from .meta_data import data_dir
-
-RUBIN_AREA = 0.25 * np.pi * 649**2  # cm^2
+from .utils import RUBIN_AREA
 
 
 __all__ = ['SkyModel', 'SkyGradient']
 
 
 class SkyModel:
     """Interface to rubin_sim.skybrightness model."""
-    def __init__(self, exptime, mjd, bandpass, eff_area=RUBIN_AREA):
+    def __init__(self, exptime, mjd, bandpass, pupil_area=RUBIN_AREA, logger=None):
         """
         Parameters
         ----------
@@ -27,15 +25,22 @@ def __init__(self, exptime, mjd, bandpass, eff_area=RUBIN_AREA):
             MJD of observation.
         bandpass : `galsim.Bandpass`
             Bandpass to use for flux calculation.
-        eff_area : `float`
-            Collecting area of telescope in cm^2. Default: Rubin value from
-            https://confluence.lsstcorp.org/display/LKB/LSST+Key+Numbers
+        pupil_area : `float`
+            Collecting area of telescope in cm^2.  The default value
+            uses R_outer=418 cm and R_inner=255 cm.
+        logger : Logger object.
         """
         from rubin_sim import skybrightness
+        logger = galsim.config.LoggerWrapper(logger)
         self.exptime = exptime
         self.mjd = mjd
-        self.eff_area = eff_area
-        self.bandpass = bandpass
+        self.pupil_area = pupil_area
+        if hasattr(bandpass, "bp_hardware"):
+            self.bandpass = bandpass.bp_hardware
+        else:
+            logger.info("A separate hardware bandpass is not available. "
+                        "Using the total bandpass for the sky level.")
+            self.bandpass = bandpass
         self._rubin_sim_sky_model = skybrightness.SkyModel()
 
     def get_sky_level(self, skyCoord):
@@ -76,7 +81,7 @@ def get_sky_level(self, skyCoord):
         flux = sed.calculateFlux(self.bandpass)
 
         # Return photons/arcsec^2
-        value = flux * self.eff_area * self.exptime
+        value = flux * self.pupil_area * self.exptime
         return value
 
 
@@ -112,6 +117,7 @@ def __call__(self, x, y):
         """
         return (self.a*x + self.b*y + self.c)/self.sky_level_center
 
+
 class CCD_Fringing:
     """
     Class generates normalized fringing map.
@@ -164,7 +170,6 @@ def generate_heightfield(self, fractal_dimension=2.5, n=4096):
 
         return np.fft.ifft2(A)
 
-
     def simulate_fringes(self, amp=0.002):
         '''
         # Generate random fringing pattern from a heightfield
@@ -183,7 +188,6 @@ def simulate_fringes(self, amp=0.002):
         #Z  += 1
         return(Z)
 
-
     def fringe_variation_level(self):
         '''
         Function implementing temporal and spatial variation of fringing.
@@ -237,18 +241,20 @@ def calculate_fringe_amplitude(self,x,y,amplitude = 0.002):
 
         return (interp_func((x,y)))
 
+
 class SkyModelLoader(InputLoader):
     """
     Class to load a SkyModel object.
     """
     def getKwargs(self, config, base, logger):
         req = {'exptime': float,
                'mjd': float}
-        opt = {'eff_area': float}
+        opt = {'pupil_area': float}
         kwargs, safe = galsim.config.GetAllParams(config, base, req=req, opt=opt)
         bandpass, safe1 = galsim.config.BuildBandpass(base['image'], 'bandpass', base, logger)
         safe = safe and safe1
         kwargs['bandpass'] = bandpass
+        kwargs['logger'] = logger
         return kwargs, safe
 
 
@@ -261,5 +267,6 @@ def SkyLevel(config, base, value_type):
     value = sky_model.get_sky_level(base['world_center'])
     return value, False
 
-RegisterInputType('sky_model', SkyModelLoader(SkyModel))
+
+RegisterInputType('sky_model', SkyModelLoader(SkyModel, takes_logger=True))
 RegisterValueType('SkyLevel', SkyLevel, [float], input_type='sky_model')

imsim/skycat.py

@@ -7,6 +7,7 @@
 from galsim.config import InputLoader, RegisterInputType, RegisterValueType, \
     RegisterObjectType
 from skycatalogs import skyCatalogs
+from .utils import RUBIN_AREA
 
 # Hot-fix to set interpolant for 500nm magnorm bandpass in skyCatalogs.
 # Once this is fixed in a skyCatalogs release, we can remove it here.
@@ -20,14 +21,10 @@
 
 class SkyCatalogInterface:
     """Interface to skyCatalogs package."""
-    # Rubin effective area computed using numbers at
-    # https://confluence.lsstcorp.org/display/LKB/LSST+Key+Numbers
-    _eff_area = 0.25 * np.pi * 649**2  # cm^2
-
     def __init__(self, file_name, wcs, band, mjd, xsize=4096, ysize=4096,
                  obj_types=None, skycatalog_root=None, edge_pix=100,
-                 max_flux=None, logger=None, apply_dc2_dilation=False,
-                 approx_nobjects=None):
+                 pupil_area=RUBIN_AREA, max_flux=None, logger=None,
+                 apply_dc2_dilation=False, approx_nobjects=None):
         """
         Parameters
         ----------
@@ -53,6 +50,9 @@ def __init__(self, file_name, wcs, band, mjd, xsize=4096, ysize=4096,
         edge_pix : float [100]
             Size in pixels of the buffer region around nominal image
             to consider objects.
+        pupil_area : float [RUBIN_AREA]
+            The area of the telescope pupil in cm**2.  The default value
+            uses R_outer=418 cm and R_inner=255 cm.
         max_flux : float [None]
             If object flux exceeds max_flux, the return None for that object.
             if max_flux == None, then don't apply a maximum flux cut.
@@ -81,6 +81,7 @@ def __init__(self, file_name, wcs, band, mjd, xsize=4096, ysize=4096,
         else:
             self.skycatalog_root = skycatalog_root
         self.edge_pix = edge_pix
+        self.pupil_area = pupil_area
         self.max_flux = max_flux
         self.logger = galsim.config.LoggerWrapper(logger)
         self.apply_dc2_dilation = apply_dc2_dilation
@@ -190,7 +191,7 @@ def getObj(self, index, gsparams=None, rng=None, exptime=30):
         for component in gsobjs:
             if component in seds:
                 gs_obj_list.append(gsobjs[component]*seds[component]
-                                   *exptime*self._eff_area)
+                                   *exptime*self.pupil_area)
 
         if not gs_obj_list:
             return None
@@ -202,7 +203,7 @@ def getObj(self, index, gsparams=None, rng=None, exptime=30):
 
         # Compute the flux or get the cached value.
         gs_object.flux \
-            = skycat_obj.get_LSST_flux(self.band, mjd=self.mjd)*exptime*self._eff_area
+            = skycat_obj.get_LSST_flux(self.band, mjd=self.mjd)*exptime*self.pupil_area
 
         if ((self.max_flux is not None and gs_object.flux > self.max_flux)
             or gs_object.flux < 0):
@@ -224,6 +225,7 @@ def getKwargs(self, config, base, logger):
                'apply_dc2_dilation': bool,
                'approx_nobjects': int,
                'mjd': float,
+               'pupil_area': float,
               }
         kwargs, safe = galsim.config.GetAllParams(config, base, req=req,
                                                   opt=opt)

imsim/utils.py

@@ -16,6 +16,8 @@
 iers.conf.auto_download = False
 iers.conf.iers_degraded_accuracy = 'ignore'
 
+RUBIN_AREA = np.pi * (418.**2 - 255.**2)  # cm^2
+
 
 def ignore_erfa_warnings(func):
     @wraps(func)

tests/data/reference_sky_levels.json

@@ -1 +1 @@
-{"u": 966.8740383242452, "g": 9565.16844049534, "r": 20335.64330578381, "i": 39770.05348346839, "z": 66820.68250072759, "y": 62387.30799288845}
+{"u": 1882.1390687502703, "g": 13292.23529652257, "r": 26707.726378502637, "i": 47153.44660256306, "z": 75328.82015750805, "y": 63810.809646841895}

tests/sky_level_reference_values.py

@@ -7,7 +7,7 @@
 import imsim
 
 
-RUBIN_AREA = 0.25 * np.pi * 649**2  # cm^2
+RUBIN_AREA = np.pi * (418.**2 - 255.**2)  # cm^2
 
 ra = 54.9348753510528
 dec = -35.8385705255579
@@ -24,7 +24,7 @@
 
 sky_levels = {}
 for band in 'ugrizy':
-    bandpass = imsim.RubinBandpass(band)
+    bandpass = imsim.RubinBandpass(band, camera='LsstCam', det_name='R22_S11').bp_hardware
     wave, spec = sky_model.return_wave_spec()
     lut = galsim.LookupTable(wave, spec[0])
     sed = galsim.SED(lut, wave_type='nm', flux_type='flambda')

tests/test_instcat_parser.py

@@ -275,7 +275,7 @@ def test_object_extraction_stars(self):
                 # Instead, this basically recapitulates the calculation in the InstCatalog class.
                 magnorm = cat.getMagNorm(i)
                 flux = np.exp(-0.9210340371976184 * magnorm)
-                rubin_area = 0.25 * np.pi * 649**2 # cm^2
+                rubin_area = np.pi * (418**2 - 255**2) # cm^2
                 exptime = 30
                 fAt = flux * rubin_area * exptime
                 sed = cat.getSED(i)
@@ -459,7 +459,7 @@ def test_object_extraction_galaxies(self):
                 # Instead, this basically recapitulates the calculation in the InstCatalog class.
                 magnorm = cat.getMagNorm(i)
                 flux = np.exp(-0.9210340371976184 * magnorm)
-                rubin_area = 0.25 * np.pi * 649**2 # cm^2
+                rubin_area = np.pi * (418**2 - 255**2) # cm^2
                 exptime = 30
                 fAt = flux * rubin_area * exptime
                 sed = cat.getSED(i) # This applies the redshift internally.

tests/test_sky_model.py

@@ -30,7 +30,7 @@ def test_sky_model():
         expected_sky_levels = json.load(fobj)
 
     for band in 'ugrizy':
-        bandpass = RubinBandpass(band)
+        bandpass = RubinBandpass(band, camera='LsstCam', det_name='R22_S11')
         sky_model = SkyModel(exptime, mjd, bandpass)
         sky_level = sky_model.get_sky_level(skyCoord)
         np.testing.assert_approx_equal(sky_level, expected_sky_levels[band],
@@ -41,7 +41,7 @@ def test_sky_model():
     # presumably innocuous reasons (makes a difference ~part per 10000).
     # Test still passes at significant=3.
     for band in 'ugrizy':
-        bandpass = RubinBandpass(band, airmass=1.2)
+        bandpass = RubinBandpass(band, airmass=1.2, camera='LsstCam', det_name='R22_S11')
         sky_model = SkyModel(exptime, mjd, bandpass)
         sky_level = sky_model.get_sky_level(skyCoord)
         np.testing.assert_approx_equal(sky_level, expected_sky_levels[band],

tests/test_stamp.py

@@ -198,6 +198,9 @@ def test_stamp_sizes():
         input.atm_psf.screen_size: 40.96
         input.checkpoint: ""
         input.sky_catalog.file_name: data/sky_cat_9683.yaml
+        input.sky_catalog.pupil_area:
+            type: Eval
+            str: "0.25 * np.pi * 649**2"
         input.opsim_data.file_name: data/small_opsim_9683.db
         input.opsim_data.visit: 449053
         input.tree_rings.only_dets: [R22_S11]
@@ -463,6 +466,9 @@ def test_stamp_bandpass_airmass():
         input.atm_psf.screen_size: 40.96
         input.checkpoint: ""
         input.sky_catalog.file_name: data/sky_cat_9683.yaml
+        input.sky_catalog.pupil_area:
+            type: Eval
+            str: "0.25 * np.pi * 649**2"
         input.opsim_data.file_name: data/small_opsim_9683.db
         input.opsim_data.visit: 449053
         input.tree_rings.only_dets: [R22_S11]