Don't err if SiliconSensor accumulates on image outside range of tree…

…ring function
GalSim-developers · Oct 31, 2023 · f6c2f70 · f6c2f70
1 parent 1b56bb4
commit f6c2f70
Show file tree

Hide file tree

Showing 2 changed files with 27 additions and 4 deletions.
diff --git a/src/Silicon.cpp b/src/Silicon.cpp
@@ -550,9 +550,9 @@ namespace galsim {
             double ty = (double)j + poly[n].y - _treeRingCenter.y + (double)orig_center.y;
             xdbg<<"tx,ty = "<<tx<<','<<ty<<std::endl;
             double r = sqrt(tx * tx + ty * ty);
-            shift = _tr_radial_table.lookup(r);
-            xdbg<<"r = "<<r<<", shift = "<<shift<<std::endl;
-            if (r > 0) {
+            if (r > 0 && r < _tr_radial_table.argMax()) {
+                shift = _tr_radial_table.lookup(r);
+                xdbg<<"r = "<<r<<", shift = "<<shift<<std::endl;
                 // Shifts are along the radial vector in direction of the doping gradient
                 double dx = shift * tx / r;
                 double dy = shift * ty / r;
@@ -583,7 +583,7 @@ namespace galsim {
                              double ty = (double)j + p.y - _treeRingCenter.y + (double)orig_center.y;
                              //xdbg<<"tx,ty = "<<tx<<','<<ty<<std::endl;
                              double r = sqrt(tx * tx + ty * ty);
-                             if (r > 0) {
+                             if (r > 0 && r < _tr_radial_table.argMax()) {
                                 double shift = _tr_radial_table.lookup(r);
                                 //xdbg<<"r = "<<r<<", shift = "<<shift<<std::endl;
                                 // Shifts are along the radial vector in direction of the doping gradient

diff --git a/tests/test_sensor.py b/tests/test_sensor.py
@@ -887,6 +887,29 @@ def test_treerings():
     # Mostly checking that there aren't any nan's here from division by 0.
     assert np.min(areas8.array) > 0
 
+    # Also check that things behave sensibly if the stamp is outside the arg range of
+    # the treering function
+    # tr6 has a max radius of 2000.
+    im.setCenter(2000,2000)
+    obj.drawImage(im, method='phot', sensor=sensor6, rng=rng6)
+    print('im.sum = ',im.array.sum(), im.added_flux)
+
+    areas6a = sensor6.calculate_pixel_areas(im, use_flux=True)
+    assert np.min(areas6a.array) > 0
+    print('areas6a = ',areas6a.array)
+    areas6b = sensor6.calculate_pixel_areas(im, use_flux=False)
+    assert np.min(areas6b.array) > 0
+    # When the stamp is outside the range of the treering function, the areas will be
+    # identical for the use_flux=False case.
+    print('areas6b = ',areas6b.array)
+    print('min, max = ',areas6b.array.min(), areas6b.array.max())
+    assert np.all(areas6b.array == areas6b.array[0,0])
+    # But not when use_flux=True
+    assert not np.all(areas6a.array == areas6a.array[0,0])
+    print('mean, std when use_flux=True: ', np.mean(areas6a.array), np.std(areas6a.array))
+    np.testing.assert_allclose(np.mean(areas6a.array), np.mean(areas6b.array),
+                               rtol=np.std(areas6a.array)/np.sqrt(np.prod(areas6b.array.shape))*3)
+
 
 @timer
 def test_resume():