diff --git a/src/Silicon.cpp b/src/Silicon.cpp index 1521c676ad..301cea32b2 100644 --- a/src/Silicon.cpp +++ b/src/Silicon.cpp @@ -596,7 +596,15 @@ namespace galsim { if (photons.hasAllocatedWavelengths()) { double lambda = photons.getWavelength(i); // in nm // Lookup the absorption length in the imported table - double abs_length = _abs_length_table.lookup(lambda); // in microns + double abs_length; + try { + abs_length = _abs_length_table.lookup(lambda); // in microns + } catch (std::runtime_error) { + if (lambda <= _abs_length_table.argMin()) + abs_length = _abs_length_table.lookup(_abs_length_table.argMin()); + else + abs_length = _abs_length_table.lookup(_abs_length_table.argMax()); + } si_length = -abs_length * log(1.0 - randomNumber); // in microns #ifdef DEBUGLOGGING if (i % 1000 == 0) { diff --git a/tests/test_sensor.py b/tests/test_sensor.py index c81b8fb69c..27f34d3229 100644 --- a/tests/test_sensor.py +++ b/tests/test_sensor.py @@ -643,6 +643,44 @@ def test_sensor_wavelengths_and_angles(): print('r1 = %f, r4 = %f, 2*sigma_r = %f'%(r1,r4,2.*sigma_r)) assert r4 > r1 +@timer +def test_bad_wavelengths(): + """If a wavelength is outside the valid range of our tabulated absorption length table, + make sure it does something reasonable, rather than abort. + """ + rng = galsim.BaseDeviate(1234) + + nphot = 7 + x = rng.np.uniform(10,20, size=nphot) + y = rng.np.uniform(10,20, size=nphot) + flux = np.ones(nphot) + dxdz = rng.np.uniform(0,0.5, size=nphot) + dydz = rng.np.uniform(0,0.5, size=nphot) + + # The original bug that triggered this test involved photons with wavelength=0. + # That's not physically possible, but now this works, using the min or max wavelength + # in the lookup table for any photons that are too blue or too red respectively. + # Note: the valid range is [255, 1450] + wave = [0.0, 55, 255., 800., 1450., 4000., np.inf] + photons = galsim.PhotonArray(nphot, x=x, y=y, flux=flux, dxdz=dxdz, dydz=dydz, wavelength=wave) + + image = galsim.Image(32,32) + sensor = galsim.SiliconSensor(name='lsst_itl_50_8', rng=rng) + sensor.accumulate(photons, image) + + # The real test is just that that didn't throw an exception. + # But check that all the photons were put somewhere. + assert np.sum(image.array) == nphot + + # The original code that triggered this was actually putting flux=0 for the bad wavelengths, + # so let's do that too. + flux = [0, 0, 1, 1, 1, 0, 0] + photons = galsim.PhotonArray(nphot, x=x, y=y, flux=flux, dxdz=dxdz, dydz=dydz, wavelength=wave) + image.setZero() + sensor.accumulate(photons, image) + assert np.sum(image.array) == np.sum(flux) + + @timer def test_bf_slopes(): """Test the brighter-fatter slopes