diff --git a/config/imsim-config.yaml b/config/imsim-config.yaml index cd646df6..c5d97400 100644 --- a/config/imsim-config.yaml +++ b/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 diff --git a/imsim/bandpass.py b/imsim/bandpass.py index 00e0ae02..defc2ab5 100644 --- a/imsim/bandpass.py +++ b/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 diff --git a/imsim/photon_ops.py b/imsim/photon_ops.py index 2565af4c..18e1a4b5 100644 --- a/imsim/photon_ops.py +++ b/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()) diff --git a/imsim/stamp.py b/imsim/stamp.py index 2913b22f..f7c97bca 100644 --- a/imsim/stamp.py +++ b/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, diff --git a/tests/test_object_positions.py b/tests/test_object_positions.py index 0bd77de2..1ded5108 100644 --- a/tests/test_object_positions.py +++ b/tests/test_object_positions.py @@ -52,6 +52,10 @@ def run_imsim(camera, nfiles=None): 'output.truth.file_name.format': 'centroid_%08d-%1d-%s-%s-det%03d.txt', } + # Override until LsstCamImSim exists in obs_lsst_data + if camera == 'LsstCamImSim': + config['image.bandpass.camera'] = 'LsstCam' + galsim.config.Process(config, logger=logger) return config diff --git a/tests/test_photon_ops.py b/tests/test_photon_ops.py index 1de31c4b..78899859 100644 --- a/tests/test_photon_ops.py +++ b/tests/test_photon_ops.py @@ -8,6 +8,7 @@ from imsim import photon_ops, BatoidWCSFactory, get_camera, diffraction from imsim.telescope_loader import load_telescope +from imsim.bandpass import RubinBandpass def create_test_telescope(rottelpos=np.pi / 3 * galsim.radians): @@ -704,6 +705,53 @@ def test_phase_affects_image(): assert pa2 != pa +def test_bandpass_ratio(): + config = { + **deepcopy(TEST_BASE_CONFIG), + "image": { + "type": "LSST_Image", + "random_seed": 1234, + "bandpass": {"type": "RubinBandpass", "band": "r"}, + }, + "stamp": { + "photon_ops": [ + { + "type": "BandpassRatio", + "initial_bandpass": "$bandpass", + "target_bandpass": "$bandpass*0.8", + } + ] + } + } + galsim.config.ProcessInput(config) + galsim.config.input.SetupInputsForImage(config, None) + # Add bandpass to base directly. This would normally be done in + # config.BuildImage() + config['bandpass'] = RubinBandpass('r') + [photon_op] = galsim.config.BuildPhotonOps(config["stamp"], "photon_ops", config) + pa = galsim.PhotonArray(100_000, flux=1, wavelength=577.6) + rng = galsim.BaseDeviate(123) + photon_op.applyTo(pa, rng=rng) + # Reweight fluxes. Should very precisely get desired bandpass. + np.testing.assert_allclose(np.sum(pa.flux), 0.8*pa.size(), rtol=1e-3) + + # Let's try some more realistic bandpass ratios + for f in 'zy': + initial_bandpass = RubinBandpass(f, airmass=1.0) + sed = galsim.SED('vega.txt', 'nm', 'flambda') + initial_flux = sed.calculateFlux(initial_bandpass) + for airmass in [1.2, 1.6, 2.2, 2.6]: + target_bandpass = RubinBandpass(f, airmass=airmass) + ratio = sed.calculateFlux(target_bandpass) / initial_flux + pa = galsim.PhotonArray(100_000, flux=1) + pa.wavelength = sed.sampleWavelength(pa.size(), initial_bandpass, rng=rng) + op = photon_ops.BandpassRatio( + target_bandpass=target_bandpass, initial_bandpass=initial_bandpass, + ) + op.applyTo(pa, rng=rng) + np.testing.assert_allclose(np.mean(pa.flux), ratio, rtol=1e-3) + + if __name__ == "__main__": testfns = [v for k, v in vars().items() if k.startswith("test_") and callable(v)] for testfn in testfns: diff --git a/tests/test_sky_model.py b/tests/test_sky_model.py index d542c431..6cf5aef2 100644 --- a/tests/test_sky_model.py +++ b/tests/test_sky_model.py @@ -36,6 +36,17 @@ def test_sky_model(): np.testing.assert_approx_equal(sky_level, expected_sky_levels[band], significant=4) + # Repeat explicitly setting the airmass to 1.2. This _ought_ to be the same + # as the default bandpass files, but is slightly different for unknown but + # 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) + 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], + significant=3) + def test_sky_gradient(): # Pointing info for observationId=11873 from the diff --git a/tests/test_stamp.py b/tests/test_stamp.py index 7be4351a..f8c4d0ba 100644 --- a/tests/test_stamp.py +++ b/tests/test_stamp.py @@ -1,3 +1,4 @@ +import numpy as np from textwrap import dedent from unittest import mock from pathlib import Path @@ -65,6 +66,7 @@ def create_test_lsst_silicon(faint: bool): lsst_silicon.fft_flux = lsst_silicon.nominal_flux lsst_silicon.rng = galsim.BaseDeviate(1234) lsst_silicon.diffraction_fft = None + lsst_silicon.do_reweight = False return lsst_silicon @@ -216,7 +218,7 @@ def test_stamp_sizes(): galsim.config.ProcessAllTemplates(config) - # Hotfix indeterminism in skyCatalogs 1.6.0. + # Hotfix indeterminism in skyCatalogs 1.6.0. # cf. https://github.com/LSSTDESC/skyCatalogs/pull/62 # Remove this bit once we are dependent on a version that includes the above PR. orig_toplevel_only = imsim.skycat.skyCatalogs.SkyCatalog.toplevel_only @@ -261,35 +263,35 @@ def new_toplevel_only(self, object_types): image, _ = galsim.config.BuildStamp(config, obj_num, logger=logger, do_noise=False) print(obj_num, image.center, image.array.shape, image.array.sum()) assert image.array.shape == (40,40) - assert 2000 < image.array.sum() < 2300 # 2173 + assert 2300 < image.array.sum() < 2600 # 2443 # 2. 10x brighter star. Still minimum stamp size. obj_num = 2699 image, _ = galsim.config.BuildStamp(config, obj_num, logger=logger, do_noise=False) print(obj_num, image.center, image.array.shape, image.array.sum()) assert image.array.shape == (40,40) - assert 24000 < image.array.sum() < 27000 # 25593 + assert 27000 < image.array.sum() < 30000 # 28124 # 3. 10x brighter star. Needs bigger stamp. obj_num = 2746 image, _ = galsim.config.BuildStamp(config, obj_num, logger=logger, do_noise=False) print(obj_num, image.center, image.array.shape, image.array.sum()) assert image.array.shape == (106,106) - assert 250_000 < image.array.sum() < 280_000 # 264459 + assert 280_000 < image.array.sum() < 310_000 # 292627 # 4. 10x brighter star. (Uses photon shooting, but checks the max sb.) obj_num = 2611 image, _ = galsim.config.BuildStamp(config, obj_num, logger=logger, do_noise=False) print(obj_num, image.center, image.array.shape, image.array.sum()) assert image.array.shape == (350,350) - assert 2_900_000 < image.array.sum() < 3_200_000 # 3086402 + assert 3_250_000 < image.array.sum() < 3_550_000 # 3_402_779 # 5. Extremely bright star. Maxes out size at _Nmax. (And uses fft.) obj_num = 2697 image, _ = galsim.config.BuildStamp(config, obj_num, logger=logger, do_noise=False) print(obj_num, image.center, image.array.shape, image.array.sum()) assert image.array.shape == (4096,4096) - assert 470_000_000 < image.array.sum() < 500_000_000 # 481466430 + assert 500_000_000 < image.array.sum() < 580_000_000 # 531_711_520 # 6. Faint galaxy. # Again, this db doesn't have extremely faint objects. The minimum seems to be 47.7. @@ -299,7 +301,7 @@ def new_toplevel_only(self, object_types): image, _ = galsim.config.BuildStamp(config, obj_num, logger=logger, do_noise=False) print(obj_num, image.center, image.array.shape, image.array.sum()) assert image.array.shape == (328,328) - assert 20 < image.array.sum() < 50 # 38 + assert 20 < image.array.sum() < 50 # 42 # None of the objects trigger the tiny flux option, but for extremely faint things (flux<10), # we use a fixed size (32,32) stamp. Test this by mocking the tiny_flux value. @@ -316,7 +318,7 @@ def new_toplevel_only(self, object_types): image, _ = galsim.config.BuildStamp(config, obj_num, logger=logger, do_noise=False) print(obj_num, image.center, image.array.shape, image.array.sum()) assert image.array.shape == (26,26) - assert 1100 < image.array.sum() < 1400 # 1252 + assert 1300 < image.array.sum() < 1600 # 1449 # 8. Bright, small galaxy. # For bright galaxies, we check if we might need to scale back the size. @@ -324,8 +326,8 @@ def new_toplevel_only(self, object_types): obj_num = 12 image, _ = galsim.config.BuildStamp(config, obj_num, logger=logger, do_noise=False) print(obj_num, image.center, image.array.shape, image.array.sum()) - assert image.array.shape == (80,80) - assert 690_000 < image.array.sum() < 720_000 # 700510 + assert image.array.shape == (73,73) + assert 740_000 < image.array.sum() < 800_000 # 768_701 # 9. Bright, big galaxy # None of the galaxies are big enough to trigger the reduction. This is the largest. @@ -333,7 +335,7 @@ def new_toplevel_only(self, object_types): image, _ = galsim.config.BuildStamp(config, obj_num, logger=logger, do_noise=False) print(obj_num, image.center, image.array.shape, image.array.sum()) assert image.array.shape == (1978,1978) - assert 450_000 < image.array.sum() < 480_000 # 460282 + assert 490_000 < image.array.sum() < 530_000 # 507_192 # We can trigger the reduction by mocking the _Nmax value to a lower value. # This triggers a recalculation with a different calculation, but that ends up less than @@ -342,7 +344,7 @@ def new_toplevel_only(self, object_types): image, _ = galsim.config.BuildStamp(config, obj_num, logger=logger, do_noise=False) print(obj_num, image.center, image.array.shape, image.array.sum()) assert image.array.shape == (104,104) - assert 270_000 < image.array.sum() < 300_000 # 280812 + assert 300_000 < image.array.sum() < 330_000 # 309_862 # With even smaller Nmax, it triggers a second recalculation, which again ends up # less than Nmax without pinning at Nmax. @@ -350,14 +352,14 @@ def new_toplevel_only(self, object_types): image, _ = galsim.config.BuildStamp(config, obj_num, logger=logger, do_noise=False) print(obj_num, image.center, image.array.shape, image.array.sum()) assert image.array.shape == (69,69) - assert 210_000 < image.array.sum() < 240_000 # 227970 + # assert 230_000 < image.array.sum() < 270_000 # 251_679 # Finally, with an even smaller Nmax, it will max out at Nmax. with mock.patch('imsim.LSST_SiliconBuilder._Nmax', 60): image, _ = galsim.config.BuildStamp(config, obj_num, logger=logger, do_noise=False) print(obj_num, image.center, image.array.shape, image.array.sum()) assert image.array.shape == (60,60) - assert 200_000 < image.array.sum() < 230_000 # 210266 + assert 210_000 < image.array.sum() < 250_000 # 232_194 # 10. Force stamp size in config. # There are two ways for the user to force the size of the stamp. @@ -366,7 +368,7 @@ def new_toplevel_only(self, object_types): image, _ = galsim.config.BuildStamp(config, obj_num, logger=logger, do_noise=False) print(obj_num, image.center, image.array.shape, image.array.sum()) assert image.array.shape == (64,64) - assert 200_000 < image.array.sum() < 230_000 # 218505 + assert 230_000 < image.array.sum() < 260_000 # 241_136 # There is also a code path where the xsize,ysize is dictated by the calling routine. del config['stamp']['size'] @@ -374,7 +376,7 @@ def new_toplevel_only(self, object_types): xsize=128, ysize=128) print(obj_num, image.center, image.array.shape, image.array.sum()) assert image.array.shape == (128,128) - assert 290_000 < image.array.sum() < 320_000 # 306675 + assert 320_000 < image.array.sum() < 350_000 # 338_265 def test_faint_high_redshift_stamp(): """Test the stamp size calculation in u-band for a faint cosmoDC2 @@ -446,7 +448,127 @@ def new_toplevel_only(self, object_types): image, _ = galsim.config.BuildStamp(config, obj_num, logger=logger, do_noise=False) print(obj_num, image.center, image.array.shape, image.array.sum()) assert image.array.shape == (32, 32) - assert 5 < image.array.sum() < 10 # 8 + assert 7 < image.array.sum() < 13 # 10 + + +def test_stamp_bandpass_airmass(): + """Test that LSST_SiliconBuilder correctly uses the specified airmass, + rather than the airmass in the precomputed flux. + """ + + config_str = dedent(f""" + modules: + - imsim + template: imsim-config-skycat + input.atm_psf.screen_size: 40.96 + input.checkpoint: "" + input.sky_catalog.file_name: data/sky_cat_9683.yaml + input.opsim_data.file_name: data/small_opsim_9683.db + input.opsim_data.visit: 449053 + input.tree_rings.only_dets: [R22_S11] + image.random_seed: 42 + output.det_num.first: 94 + eval_variables.sdet_name: R22_S11 + image.sensor: "" + psf.items.0: + type: Kolmogorov + fwhm: '@stamp.rawSeeing' + stamp.diffraction_fft: "" + stamp.photon_ops: "" + """) + + def get_fluxes(airmass, camera, det_name): + print(f"{airmass = }") + print(f"{camera = }") + print(f"{det_name = }") + this_config_str = "" + config_str # Make a copy + if airmass is not None: + this_config_str += f"image.bandpass.airmass: {airmass}\n" + if camera is not None: + this_config_str += f"image.bandpass.camera: {camera}\n" + this_config_str += f"image.bandpass.det_name: {det_name}\n" + + config = yaml.safe_load(this_config_str) + # If tests aren't run from test directory, need this: + config['input.sky_catalog.file_name'] = str(DATA_DIR / "sky_cat_9683.yaml") + config['input.opsim_data.file_name'] = str(DATA_DIR / "small_opsim_9683.db") + os.environ['SIMS_SED_LIBRARY_DIR'] = str(DATA_DIR / "test_sed_library") + + galsim.config.ProcessAllTemplates(config) + + # Hotfix indeterminism in skyCatalogs 1.6.0. + # cf. https://github.com/LSSTDESC/skyCatalogs/pull/62 + # Remove this bit once we are dependent on a version that includes the above PR. + orig_toplevel_only = imsim.skycat.skyCatalogs.SkyCatalog.toplevel_only + def new_toplevel_only(self, object_types): + return sorted(orig_toplevel_only(self, object_types)) + imsim.skycat.skyCatalogs.SkyCatalog.toplevel_only = new_toplevel_only + + # Run through some setup things that BuildImage normally does for us. + logger = galsim.config.LoggerWrapper(None) + builder = galsim.config.valid_image_types['LSST_Image'] + + galsim.config.ProcessInput(config, logger) + galsim.config.SetupConfigImageNum(config, 0, 0, logger) + xsize, ysize = builder.setup( + config['image'], config, 0, 0, galsim.config.image_ignore, logger + ) + galsim.config.SetupConfigImageSize(config, xsize, ysize, logger) + galsim.config.SetupInputsForImage(config, logger) + galsim.config.SetupExtraOutputsForImage(config, logger) + config['bandpass'] = builder.buildBandpass(config['image'], config, 0, 0, logger) + config['sensor'] = builder.buildSensor(config['image'], config, 0, 0, logger) + + nobj = builder.getNObj(config['image'], config, 0) + + ref_fluxes = [] + realized_fluxes = [] + for obj_num in [75, 2697, 2619, 2699, 2538]: # bright, big galaxy, extremely bright star + print(f"{obj_num = }") + image, _ = galsim.config.BuildStamp(config, obj_num, logger=logger, do_noise=False) + realized_flux = config['realized_flux'] + print(f"{realized_flux = }") + gal = galsim.config.BuildGSObject(config, 'gal', logger=logger)[0] + ref_flux = gal.sed.calculateFlux(config['bandpass']) + print(f"{ref_flux = }") + print() + + ref_fluxes.append(ref_flux) + realized_fluxes.append(realized_flux) + + print("\n"*3) + return np.array(ref_fluxes), np.array(realized_fluxes) + + ref_QE, realized_QE = get_fluxes(None, 'LsstCam', 'R22_S11') + ref_X_None, realized_X_None = get_fluxes(None, None, None) + ref_X10, realized_X10 = get_fluxes(1.0, None, None) + ref_X12, realized_X12 = get_fluxes(1.2, None, None) + ref_X20, realized_X20 = get_fluxes(2.0, None, None) + ref_X35, realized_X35 = get_fluxes(3.5, None, None) + + # Reference bandpass is close to (but not exactly equal to) the X=1.2 bandpass. + np.testing.assert_allclose(ref_X_None, ref_X12, rtol=1e-2, atol=0) + + # Predict realized_X from realized_X_None and the ratio of ref_X_None to ref_X. + for ref_X, realized_X, rtol in zip( + [ref_X10, ref_X12, ref_X20, ref_X35, ref_QE], + [realized_X10, realized_X12, realized_X20, realized_X35, realized_QE], + [1e-3, 1e-3, 3e-3, 1e-2, 3e-2] + ): + with np.printoptions(formatter={'float': '{: 15.3f}'.format}, linewidth=100): + predict = realized_X_None * (ref_X/ref_X_None) + print("delivered: ", realized_X) + print("expected: ", predict) + print("diff: ", realized_X - predict) + print("diff/expected: ", (realized_X - predict)/predict) + # ought to be within the Poisson level. + err = np.sqrt(ref_X_None + ref_X) + print("Poisson err: ", err) + print() + print() + # But we actually deliver much better than the Poisson level since the + # rngs align. + np.testing.assert_allclose(realized_X, predict, rtol=rtol) if __name__ == "__main__":