Merge pull request #195 from LSSTDESC/u/jchiang/dark_current_read_noi…

…se_fixes

dark current and read noise fixes

python/desc/imsim/camera_readout.py

@@ -80,6 +80,11 @@ def __init__(self, image_array, exptime, sensor_id, visit=42, logger=None):
             logging.Logger object to use. If None, then a logger with level
             INFO will be used.
         """
+        if logger is None:
+            self.logger = get_logger('INFO')
+        else:
+            self.logger = logger
+
         self.eimage = fits.HDUList()
         self.eimage.append(fits.PrimaryHDU(image_array))
         self.eimage_data = self.eimage[0].data.transpose()
@@ -88,17 +93,13 @@ def __init__(self, image_array, exptime, sensor_id, visit=42, logger=None):
         self.sensor_id = sensor_id
         self.visit = visit
 
-        self._seed = None
+        self._rng = None
+        self.seed = None
 
         self.camera_info = CameraInfo()
 
         self._make_amp_images()
 
-        if logger is None:
-            self.logger = get_logger('INFO')
-        else:
-            self.logger = logger
-
         self.ratel = 0
         self.dectel = 0
         self.rotangle = 0
@@ -290,12 +291,14 @@ def _make_amp_image(self, amp_name):
         imaging_segment.getArray()[:] = data
 
         # Add dark current.
-        dark_current = config['electronics_readout']['dark_current']
-        imaging_arr = imaging_segment.getArray()
-        rng = galsim.PoissonDeviate(self.seed, dark_current*self.exptime)
-        dc_data = np.zeros(np.prod(imaging_arr.shape))
-        rng.generate(dc_data)
-        imaging_arr += dc_data.reshape(imaging_arr.shape)
+        if self.exptime > 0:
+            dark_current = config['electronics_readout']['dark_current']
+            imaging_arr = imaging_segment.getArray()
+            rng = galsim.PoissonDeviate(seed=self.rng,
+                                        mean=dark_current*self.exptime)
+            dc_data = np.zeros(np.prod(imaging_arr.shape))
+            rng.generate(dc_data)
+            imaging_arr += dc_data.reshape(imaging_arr.shape)
 
         # Add defects.
 
@@ -328,22 +331,25 @@ def _add_read_noise_and_bias(self, amp_name):
         """
         amp_info = self.camera_info.get_amp_info(amp_name)
         full_arr = self.amp_images[amp_name].getArray()
-        rng = galsim.GaussianDeviate(self.seed, amp_info.getReadNoise())
+        rng = galsim.GaussianDeviate(seed=self.rng,
+                                     sigma=amp_info.getReadNoise())
         rn_data = np.zeros(np.prod(full_arr.shape))
         rng.generate(rn_data)
         full_arr += rn_data.reshape(full_arr.shape)
         full_arr += config['electronics_readout']['bias_level']
 
     @property
-    def seed(self):
+    def rng(self):
         """
-        Random seed derived from visit and sensor id.  This is used as
-        the seed for both the read noise and dark current
-        calculations.
+        galsim.BaseDeviate seeded with an integer derived from a hash
+        of the visit and sensor id.  This is used as the underlying
+        random number generator for both the read noise and dark
+        current calculations.
         """
-        if self._seed is None:
-            self._seed = CosmicRays.generate_seed(self.visit, self.sensor_id)
-        return self._seed
+        if self._rng is None:
+            self.seed = CosmicRays.generate_seed(self.visit, self.sensor_id)
+            self._rng = galsim.BaseDeviate(self.seed)
+        return self._rng
 
     def _apply_crosstalk(self):
         """

tests/test_instcat_parser.py

@@ -202,9 +202,9 @@ def test_object_extraction_stars(self):
             fnorm = sed.calcFluxNorm(truth_data['magNorm'][i_obj], imsim_bp)
             sed.multiplyFluxNorm(fnorm)
             sed.resampleSED(wavelen_match=bp_dict.wavelenMatch)
-            a_x, b_x = sed.setupCCMab()
-            sed.addCCMDust(a_x, b_x, A_v=truth_data['Av'][i_obj],
-                           R_v=truth_data['Rv'][i_obj])
+            a_x, b_x = sed.setupCCM_ab()
+            sed.addDust(a_x, b_x, A_v=truth_data['Av'][i_obj],
+                        R_v=truth_data['Rv'][i_obj])
 
             for bp in ('u', 'g', 'r', 'i', 'z', 'y'):
                 flux = sed.calcADU(bp_dict[bp], phot_params)*phot_params.gain
@@ -331,15 +331,15 @@ def test_object_extraction_galaxies(self):
             fnorm = sed.calcFluxNorm(truth_data['magNorm'][i_obj], imsim_bp)
             sed.multiplyFluxNorm(fnorm)
 
-            a_x, b_x = sed.setupCCMab()
-            sed.addCCMDust(a_x, b_x, A_v=truth_data['internalAv'][i_obj],
-                           R_v=truth_data['internalRv'][i_obj])
+            a_x, b_x = sed.setupCCM_ab()
+            sed.addDust(a_x, b_x, A_v=truth_data['internalAv'][i_obj],
+                        R_v=truth_data['internalRv'][i_obj])
 
             sed.redshiftSED(truth_data['redshift'][i_obj], dimming=True)
             sed.resampleSED(wavelen_match=bp_dict.wavelenMatch)
-            a_x, b_x = sed.setupCCMab()
-            sed.addCCMDust(a_x, b_x, A_v=truth_data['galacticAv'][i_obj],
-                           R_v=truth_data['galacticRv'][i_obj])
+            a_x, b_x = sed.setupCCM_ab()
+            sed.addDust(a_x, b_x, A_v=truth_data['galacticAv'][i_obj],
+                        R_v=truth_data['galacticRv'][i_obj])
 
             for bp in ('u', 'g', 'r', 'i', 'z', 'y'):
                 flux = sed.calcADU(bp_dict[bp], phot_params)*phot_params.gain