Merge pull request #266 from lsst/tickets/DM-39460

DM-39460: Bug-fix for CHARGE_SUM photodiode integration method

python/lsst/ip/isr/photodiode.py

@@ -297,10 +297,17 @@ def integrateChargeSum(self):
         # The current per interval to use for baseline subtraction
         # without assuming all of the dt values are the same:
         current = charge/dt
-        # For the baseline current level, select current values < 5%
-        # of the maximum, measured relative to the overall minimum.
+        # To determine the baseline current level, exclude points with
+        # signal levels > 5% of the maximum (measured relative to the
+        # overall minimum), and extend that selection 2 entries on
+        # either side to avoid otherwise low-valued points that sample
+        # the signal ramp and which should not be included in the
+        # baseline estimate.
         dy = np.max(current) - np.min(current)
-        index = np.where(current < dy/20. + np.min(current))
-        bg_current = np.sum(charge[index])/np.sum(dt[index])
+        signal, = np.where(current > dy/20. + np.min(current))
+        imin = signal[0] - 2
+        imax = signal[-1] + 2
+        bg = np.concatenate([np.arange(0, imin), np.arange(imax, len(current))])
+        bg_current = np.sum(charge[bg])/np.sum(dt[bg])
         # Return the background-subtracted total charge.
         return np.sum(charge - bg_current*dt)

tests/test_photodiode.py

@@ -22,6 +22,8 @@
 import unittest
 import tempfile
 
+import numpy as np
+
 import lsst.afw.cameraGeom as cameraGeom
 import lsst.geom as geom
 import lsst.utils.tests
@@ -63,7 +65,7 @@ def testDataOnly(self):
         self.assertFloatsAlmostEqual(calib.integrate(), 2.88414e-10, rtol=1e-14)
         self.assertFloatsAlmostEqual(calib.integrateDirectSum(), 2.88414e-10, rtol=1e-14)
         self.assertFloatsAlmostEqual(calib.integrateTrimmedSum(), 2.88720e-10, rtol=1e-14)
-        self.assertFloatsAlmostEqual(calib.integrateChargeSum(), 3.5978848e-09, rtol=1e-14)
+        self.assertFloatsAlmostEqual(calib.integrateChargeSum(), 3.5960649e-09, rtol=1e-14)
 
         self.assertEqual(calib.timeSamples.shape, (16, ))
         self.assertEqual(calib.currentSamples.shape, (16, ))
@@ -87,7 +89,7 @@ def testFullyPopulated(self):
         self.assertFloatsAlmostEqual(calib.integrate(), 2.88414e-10, rtol=1e-14)
         self.assertFloatsAlmostEqual(calib.integrateDirectSum(), 2.88414e-10, rtol=1e-14)
         self.assertFloatsAlmostEqual(calib.integrateTrimmedSum(), 2.88720e-10, rtol=1e-14)
-        self.assertFloatsAlmostEqual(calib.integrateChargeSum(), 3.5978848e-09, rtol=1e-14)
+        self.assertFloatsAlmostEqual(calib.integrateChargeSum(), 3.5960649e-09, rtol=1e-14)
 
         self.assertEqual(calib.timeSamples.shape, (16, ))
         self.assertEqual(calib.currentSamples.shape, (16, ))
@@ -102,6 +104,25 @@ def testFullyPopulated(self):
         newPhotodiode = PhotodiodeCalib().readFits(outPath)
         self.assertEqual(calib, newPhotodiode)
 
+    def testIntegrateChargeSum(self):
+        """
+        This tests the .integrateChargeSum method in the presence of
+        an above-baseline point which is still below the initial
+        threshold for identifying points to exclude in the baseline
+        current estimate.
+        """
+        # Compute a value to insert at time 0.47.
+        test_time = 0.47
+        test_current = ((max(self.currentSeries) - min(self.currentSeries))/40.
+                        + min(self.currentSeries))
+        new_times = np.array(self.timeSeries + [test_time])
+        new_currents = np.array(self.currentSeries + [test_current])
+        # Put these arrays in time order, though this isn't strictly needed.
+        index = np.argsort(new_times)
+        calib = PhotodiodeCalib(timeSamples=new_times[index],
+                                currentSamples=new_currents[index])
+        self.assertFloatsAlmostEqual(calib.integrateChargeSum(), 3.6049277e-09, rtol=1e-14)
+
 
 class MemoryTester(lsst.utils.tests.MemoryTestCase):
     pass