error propagation for kgain multip. and photon/shot noise using ptc/ENF

corgidrp/detector.py

@@ -816,4 +816,18 @@ def create_onsky_flatfield(dataset, planet=None,band=None,up_radius=55,im_size=N
     onsky_flatfield.err=raster_com[1]
 
     return(onsky_flatfield)
+
+def ENF(g, Nem):
+    """Returns the extra-noise function (ENF).
+    
+    Args:
+        g : float 
+        EM gain.  >= 1.
+    Nem : int
+        Number of gain register cells.
     
+    Returns
+        ENF : float, extra-noise function
+    """
+    return np.sqrt(2*(g-1)*g**(-(Nem+1)/Nem) + 1/g) 
+

corgidrp/l2a_to_l2b.py

@@ -1,22 +1,33 @@
 # A file that holds the functions that transmogrify l2a data to l2b data
 import numpy as np
+from scipy.interpolate import interp1d
+
 import corgidrp.data as data
 from corgidrp.darks import build_synthesized_dark
 from corgidrp.detector import detector_areas
+from corgidrp.detector import ENF
 
-def add_photon_noise(input_dataset):
+def add_photon_noise(input_dataset, kgain, detector_params):
     """
     Propagate the photon/shot noise determined from the image signal to the error map. 
 
     Args:
        input_dataset (corgidrp.data.Dataset): a dataset of Images (L2a-level)
+       kgain (corgidrp.data.KGain): kgain calibration object
+       detector_params (corgidrp.data.DetectorParams): detector parameters calibration object
 
     Returns:
         corgidrp.data.Dataset: photon noise propagated to the image error extensions of the input dataset
     """
     # you should make a copy the dataset to start
-    phot_noise_dataset = input_dataset.copy() # necessary at all?
-
+    phot_noise_dataset = input_dataset.copy()
+
+    #get the noise from the ptc curve
+    ptc = kgain.ptc
+    #better say shot noise?
+    #should we also add the read noise as error layer?
+    nem = detector_params.params['Nem']
+
     for i, frame in enumerate(phot_noise_dataset.frames):
         try: # use measured gain if available TODO change hdr name if necessary
             em_gain = phot_noise_dataset[i].ext_hdr["EMGAIN_M"]
@@ -25,10 +36,15 @@ def add_photon_noise(input_dataset):
                 em_gain = phot_noise_dataset[i].ext_hdr["EMGAIN_A"]
             except: # otherwise use commanded EM gain
                 em_gain = phot_noise_dataset[i].ext_hdr["CMDGAIN"]
-        phot_err = np.sqrt(frame.data)
+        #estimate of photon/shot noise by interpolation of the photon transfer curve
+        interp_func = interp1d(ptc[:,0], ptc[:,1], kind='linear', fill_value='extrapolate')
+
+        phot_err = interp_func(frame.data)
+        #phot_err = np.sqrt(frame.data), simple first estimate of photon noise, which is only true for high signals
         #add excess noise in case of em_gain
         if em_gain > 1:
-            phot_err *= np.sqrt(2)           
+            nem = detector_params.params['Nem']
+            phot_err *= ENF(em_gain, nem)           
         frame.add_error_term(phot_err, "photnoise_error")
 
     new_all_err = np.array([frame.err for frame in phot_noise_dataset.frames])        
@@ -250,11 +266,11 @@ def convert_to_electrons(input_dataset, k_gain):
     kgain_cube = kgain_dataset.all_data
 
     kgain = k_gain.value #extract from caldb
-    error_frame = kgain_dataset[0].data * k_gain.error
+    error_frame = kgain_dataset[0].data * k_gain.error #kgain_cube * error?
     kgain_cube *= kgain
 
-    #scale also the old error with kgain and propagate the error of 
-    kgain_dataset.rescale_error(kgain, "kgain")
+    #scale also the old error with kgain and propagate the error 
+    kgain_dataset.rescale_error(kgain, "kgain") #should be 2/3 dim?
     kgain_dataset.add_error_term(error_frame, "kgain_error")
 
     history_msg = "data converted to detected EM electrons by kgain {0}".format(str(kgain))
@@ -432,4 +448,4 @@ def update_to_l2b(input_dataset):
     history_msg = "Updated Data Level to L2b"
     updated_dataset.update_after_processing_step(history_msg)
 
-    return updated_dataset
+    return updated_dataset

tests/test_add_phot_noise.py

@@ -2,8 +2,9 @@
 import numpy as np
 import corgidrp
 from corgidrp.mocks import create_default_headers
-from corgidrp.data import Image, Dataset
+from corgidrp.data import Image, Dataset, DetectorParams, KGain
 from corgidrp.l2a_to_l2b import add_photon_noise
+from corgidrp.detector import ENF
 import pytest
 
 old_err_tracking = corgidrp.track_individual_errors
@@ -21,19 +22,27 @@ def test_add_phot_noise():
     image2 = Image(data,pri_hdr = prhd, ext_hdr = exthd, err = err, dq = dq)
 
     dataset = Dataset([image1, image2])
-    dataset_add = add_photon_noise(dataset)
+
+    detector_params = DetectorParams({})
+    gain_value = np.array([[8.7]])
+    signal_array = np.linspace(0, 50)
+    noise_array = np.sqrt(signal_array)
+    ptc = np.column_stack([signal_array, noise_array])
+    kgain = KGain(gain_value, ptc = ptc, pri_hdr = prhd, ext_hdr = exthd, input_dataset = dataset)
+
+    dataset_add = add_photon_noise(dataset, kgain, detector_params)
     all_err = dataset.all_err
     all_err1 = dataset_add.all_err
     assert not np.array_equal(all_err, all_err1)
-    assert np.array_equal(all_err1[0,1], np.sqrt(data))
-    assert np.allclose(all_err1[0,0], np.sqrt(data + np.square(err)))
+    assert np.allclose(all_err1[0,1], np.sqrt(data), rtol = 0.01)
+    assert np.allclose(all_err1[0,0], np.sqrt(data + np.square(err)), rtol = 0.01)
     assert "noise" in str(dataset_add.frames[0].ext_hdr["HISTORY"])
     assert "photnoise_error" == dataset_add.frames[0].err_hdr["Layer_2"]
     #check that excess noise is applied
     dataset[0].ext_hdr["CMDGAIN"] = 3000
-    dataset_add1 = add_photon_noise(dataset)
+    dataset_add1 = add_photon_noise(dataset, kgain, detector_params)
     all_err2 = dataset_add1.all_err
-    assert np.array_equal(all_err2[0,1], np.sqrt(data)*np.sqrt(2))
+    assert np.allclose(all_err2[0,1], np.sqrt(data)*ENF(3000, 604), rtol =0.01)
 
     corgidrp.track_individual_errors = old_err_tracking