diff --git a/.github/pull_request_template.md b/.github/pull_request_template.md
index 09cdcbf7..8f85a7ac 100644
--- a/.github/pull_request_template.md
+++ b/.github/pull_request_template.md
@@ -16,3 +16,4 @@ Please delete options that are not relevant (and this line).
 - [ ] I have linted my code
 - [ ] I have verified that all unit tests pass in a clean environment and added new unit tests, as needed
 - [ ] I have verified that all docstrings are properly formatted and added new documentation, as needed
+- [ ] I have filled out the Unit Test Definition Table on confluence, as needed
\ No newline at end of file
diff --git a/corgidrp/caldb.py b/corgidrp/caldb.py
index f60b8d84..a38c7f9f 100644
--- a/corgidrp/caldb.py
+++ b/corgidrp/caldb.py
@@ -168,7 +168,7 @@ def _get_values_from_entry(self, entry, is_calib=True):
             drp_version,
             obsid,
             naxis1,
-            naxis2,
+            naxis2
         ]
 
         # rest are ext_hdr keys we can copy
@@ -243,7 +243,7 @@ def remove_entry(self, entry, to_disk=True):
 
     def get_calib(self, frame, dtype, to_disk=True):
         """
-        Outputs the best calibration file of the given type for the input sciene frame.
+        Outputs the best calibration file of the given type for the input science frame.
 
         Args:
             frame (corgidrp.data.Image): an image frame to request a calibration for. If None is passed in, looks for the 
@@ -288,14 +288,12 @@ def get_calib(self, frame, dtype, to_disk=True):
             options = calibdf.loc[
                 (
                     (calibdf["EXPTIME"] == frame_dict["EXPTIME"])
-                    & (calibdf["NAXIS1"] == frame_dict["NAXIS1"])
-                    & (calibdf["NAXIS2"] == frame_dict["NAXIS2"])
                 )
             ]
 
             if len(options) == 0:
-                raise ValueError("No valid Dark with EXPTIME={0} and dimension ({1},{2})"
-                                 .format(frame_dict["EXPTIME"], frame_dict["NAXIS1"], frame_dict["NAXIS2"]))
+                raise ValueError("No valid Dark with EXPTIME={0})"
+                                 .format(frame_dict["EXPTIME"]))
 
             # select the one closest in time
             result_index = np.abs(options["MJD"] - frame_dict["MJD"]).argmin()
diff --git a/corgidrp/calibrate_kgain.py b/corgidrp/calibrate_kgain.py
index 405ba47d..d4212164 100644
--- a/corgidrp/calibrate_kgain.py
+++ b/corgidrp/calibrate_kgain.py
@@ -12,7 +12,7 @@
 from corgidrp.detector import slice_section, detector_areas
 
 # Dictionary with constant kgain calibration parameters
-kgain_params= {
+kgain_params_default= {
 # ROI constants
 'rowroi1': 9,
 'rowroi2': 1000,
@@ -30,34 +30,30 @@
 'signal_bins_N': 400,
 }
 
-def check_kgain_params(
-    ):
-    """ Checks integrity of kgain parameters in the dictionary kgain_params. """
-    if 'offset_colroi1' not in kgain_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
-    if 'offset_colroi2' not in kgain_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
+def check_kgain_params(kgain_params):
+    """ Checks integrity of kgain parameters in the dictionary kgain_params. 
+    
+    Args:
+        kgain_params (dict):  Dictionary of parameters used for calibrating the k gain.
+    """
+    
     if 'rowroi1' not in kgain_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
+        raise ValueError('Missing parameter:  rowroi1.')
     if 'rowroi2' not in kgain_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
+        raise ValueError('Missing parameter:  rowroi2.')
     if 'colroi1' not in kgain_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
+        raise ValueError('Missing parameter:  colroi1.')
     if 'colroi2' not in kgain_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
+        raise ValueError('Missing parameter:  colroi2.')
     if 'rn_bins1' not in kgain_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
+        raise ValueError('Missing parameter:  rn_bins1.')
     if 'rn_bins2' not in kgain_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
+        raise ValueError('Missing parameter:  rn_bins2.')
     if 'max_DN_val' not in kgain_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
+        raise ValueError('Missing parameter:  max_DN_val.')
     if 'signal_bins_N' not in kgain_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
+        raise ValueError('Missing parameter:  signal_bins_N.')
 
-    if not isinstance(kgain_params['offset_colroi1'], (float, int)):
-        raise TypeError('offset_colroi1 is not a number')
-    if not isinstance(kgain_params['offset_colroi2'], (float, int)):
-        raise TypeError('offset_colroi2 is not a number')
     if not isinstance(kgain_params['rowroi1'], (float, int)):
         raise TypeError('rowroi1 is not a number')
     if not isinstance(kgain_params['rowroi2'], (float, int)):
@@ -281,7 +277,7 @@ def calibrate_kgain(dataset_kgain,
                     n_cal=10, n_mean=30, min_val=800, max_val=3000, binwidth=68,
                     make_plot=True,plot_outdir='figures', show_plot=False,
                     logspace_start=-1, logspace_stop=4, logspace_num=200,
-                    verbose=False, detector_regions=None):
+                    verbose=False, detector_regions=None, kgain_params=None):
     """
     kgain (e-/DN) is calculated from the means and variances
     within the defined minimum and maximum mean values. A photon transfer curve
@@ -344,12 +340,29 @@ def calibrate_kgain(dataset_kgain,
       detector_regions (dict): a dictionary of detector geometry properties.
         Keys should be as found in detector_areas in detector.py.  Defaults to
         that dictionary.
+      kgain_params (dict): (Optional) Dictionary containing row and col specifications
+        for the region of interest (indicated by 'rowroi1','rowroi2','colroi1',and 'colroi2').
+        The 'roi' needs one square region specified, and 'back' needs two square regions, 
+        where a '1' ending indicates the smaller of two values, and a '2' ending indicates the larger 
+        of two values.  The coordinates of the square region are specified by matching 
+        up as follows: (rowroi1, colroi1), (rowroi2, colroi1), etc. 
+        Also must contain:
+        'rn_bins1': lower bound of counts histogram for fitting or read noise
+        'rn_bins2': upper bound of counts histogram for fitting or read noise 
+        'max_DN_val': maximum DN value to be included in photon transfer curve (PTC)
+        'signal_bins_N': number of bins in the signal variables of PTC curve
+        Defaults to kgain_params_default included in this file.
     
     Returns:
       corgidrp.data.KGain: kgain estimate from the least-squares fit to the photon
         transfer curve (in e-/DN). The expected value of kgain for EXCAM with
         flight readout sequence should be between 8 and 9 e-/DN
     """
+    if kgain_params is None:
+        kgain_params = kgain_params_default
+        
+    check_kgain_params(kgain_params)
+
     if detector_regions is None:
         detector_regions = detector_areas
 
diff --git a/corgidrp/calibrate_nonlin.py b/corgidrp/calibrate_nonlin.py
index 199ed826..c5a6b2dc 100644
--- a/corgidrp/calibrate_nonlin.py
+++ b/corgidrp/calibrate_nonlin.py
@@ -16,7 +16,7 @@
 from corgidrp.calibrate_kgain import CalKgainException
 
 # Dictionary with constant non-linearity calibration parameters
-nonlin_params = {
+nonlin_params_default = {
     # ROI constants
     'rowroi1': 305,
     'rowroi2': 736,
@@ -44,41 +44,44 @@
     'min_mask_factor': 1.1,
     }
  
-def check_nonlin_params(
-    ):
-    """ Checks integrity of kgain parameters in the dictionary nonlin_params. """
+def check_nonlin_params(nonlin_params):
+    """ Checks integrity of kgain parameters in the dictionary nonlin_params. 
+
+    Args:
+        nonlin_params (dict):  Dictionary of parameters used for calibrating nonlinearity.
+    """
     if 'rowroi1' not in nonlin_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
+        raise ValueError('Missing parameter:  rowroi1.')
     if 'rowroi2' not in nonlin_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
+        raise ValueError('Missing parameter:  rowroi2.')
     if 'colroi1' not in nonlin_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
+        raise ValueError('Missing parameter:  colroi1.')
     if 'colroi2' not in nonlin_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
+        raise ValueError('Missing parameter:  colroi2.')
     if 'rowback11' not in nonlin_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
+        raise ValueError('Missing parameter:  rowback11.')
     if 'rowback12' not in nonlin_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
+        raise ValueError('Missing parameter:  rowback12.')
     if 'rowback21' not in nonlin_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
+        raise ValueError('Missing parameter:  rowback21.')
     if 'rowback22' not in nonlin_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
+        raise ValueError('Missing parameter:  rowback22.')
     if 'colback11' not in nonlin_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
+        raise ValueError('Missing parameter:  colback11.')
     if 'colback12' not in nonlin_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
+        raise ValueError('Missing parameter:  colback12.')
     if 'colback21' not in nonlin_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
+        raise ValueError('Missing parameter:  colback21.')
     if 'colback22' not in nonlin_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
+        raise ValueError('Missing parameter:  colback22.')
     if 'min_exp_time' not in nonlin_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
+        raise ValueError('Missing parameter:  min_exp_time.')
     if 'num_bins' not in nonlin_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
+        raise ValueError('Missing parameter:  num_bins.')
     if 'min_bin' not in nonlin_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
+        raise ValueError('Missing parameter:  min_bin.')
     if 'min_mask_factor' not in nonlin_params:
-        raise ValueError('Missing parameter in directory pointer YAML file.')
+        raise ValueError('Missing parameter:  min_mask_factor.')
     
     if not isinstance(nonlin_params['rowroi1'], (float, int)):
         raise TypeError('rowroi1 is not a number')
@@ -124,7 +127,7 @@ def calibrate_nonlin(dataset_nl,
                      pfit_upp_cutoff1 = -2, pfit_upp_cutoff2 = -3,
                      pfit_low_cutoff1 = 2, pfit_low_cutoff2 = 1,
                      make_plot=True, plot_outdir='figures', show_plot=False,
-                     verbose=False):
+                     verbose=False, nonlin_params=None):
     """
     Function that derives the non-linearity calibration table for a set of DN
     and EM values.
@@ -210,6 +213,16 @@ def calibrate_nonlin(dataset_nl,
       show_plot (bool): (Optional) display the plots. Default is False.
       verbose (bool): (Optional) display various diagnostic print messages.
         Default is False.
+      nonlin_params (dict): (Optional) Dictionary of row and col specifications
+        for the region of interest (indicated by 'roi') where the frame is illuminated and for 
+        two background regions (indicated by 'back1' and 'back2') where the frame is not illuminated.  
+        Must contain 'rowroi1','rowroi2','colroi1','colroi2','rowback11','rowback12',
+        'rowback21','rowback22','colback11','colback12','colback21',and 'colback22'.
+        The 'roi' needs one square region specified, and 'back' needs two square regions, 
+        where a '1' ending indicates the smaller of two values, and a '2' ending indicates the larger 
+        of two values.  The coordinates of each square are specified by matching 
+        up as follows: (rowroi1, colroi1), (rowroi1, colroi2), (rowback11, colback11), 
+        (rowback11, colback12), etc. Defaults to nonlin_params_default specified in this file.
     
     Returns:
       nonlin_arr (NonLinearityCalibration): 2-D array with nonlinearity values
@@ -217,6 +230,11 @@ def calibrate_nonlin(dataset_nl,
         input signal in DN is the first column. Signal values start with min_write
         and run through max_write in steps of 20 DN.
     """
+    if nonlin_params is None:
+        nonlin_params = nonlin_params_default
+        
+    check_nonlin_params(nonlin_params)
+
     # dataset_nl.all_data must be 3-D 
     if np.ndim(dataset_nl.all_data) != 3:
         raise Exception('dataset_nl.all_data must be 3-D')
diff --git a/corgidrp/combine.py b/corgidrp/combine.py
index 3235db97..7b2663f2 100644
--- a/corgidrp/combine.py
+++ b/corgidrp/combine.py
@@ -43,7 +43,7 @@ def combine_images(data_subset, err_subset, dq_subset, collapse, num_frames_scal
     
     # dq collpase: keep all flags on
     dq_collapse = np.bitwise_or.reduce(dq_subset, axis=0)
-    # except those pixels that have been replaced 
+    # except for those pixels that have been replaced with good values
     dq_collapse[np.where((dq_collapse > 0) & (~np.isnan(data_collapse)))] = 0
 
     return data_collapse, err_collapse, dq_collapse
diff --git a/corgidrp/darks.py b/corgidrp/darks.py
index 93c53f82..0abec4aa 100644
--- a/corgidrp/darks.py
+++ b/corgidrp/darks.py
@@ -105,14 +105,14 @@ def mean_combine(image_list, bpmap_list, err=False):
             sum_im += masked
         map_im += (im_m.mask == False).astype(int)
 
-    if err: # sqrt of sum of sigma**2 terms
-        sum_im = np.sqrt(sum_im)
-
     # Divide sum_im by map_im only where map_im is not equal to 0 (i.e.,
     # not masked).
     # Where map_im is equal to 0, set combined_im to zero
     comb_image = np.divide(sum_im, map_im, out=np.zeros_like(sum_im),
                             where=map_im != 0)
+   
+    if err: # (sqrt of sum of sigma**2 terms)/sqrt(n)
+        comb_image = np.sqrt(comb_image)
 
     # Mask any value that was never mapped (aka masked in every frame)
     comb_bpmap = (map_im == 0).astype(int)
@@ -133,12 +133,13 @@ def build_trad_dark(dataset, detector_params, detector_regions=None, full_frame=
     - have had masks made for cosmic rays
     - have been corrected for nonlinearity
     - have been converted from DN to e-
-    - have been desmeared if desmearing is appropriate.  Under normal
-    circumstances, darks should not be desmeared.  The only time desmearing
-    would be useful is in the unexpected case that, for example,
-    dark current is so high that it stands far above other noise that is
-    not smeared upon readout, such as clock-induced charge, 
-    fixed-pattern noise, and read noise.
+    - have NOT been desmeared. Darks should not be desmeared.  The only component 
+    of dark frames that would be subject to a smearing effect is dark current 
+    since it linearly increases with time, so the extra row read time affects 
+    the dark current per pixel.  However, illuminated images
+    would also contain this smeared dark current, so dark subtraction should 
+    remove this smeared dark current (and then desmearing may be applied to the 
+    processed image if appropriate).  
 
     Also, add_photon_noise() should NOT have been applied to the frames in
     dataset.  And note that creation of the
@@ -268,12 +269,13 @@ def calibrate_darks_lsq(dataset, detector_params, detector_regions=None):
     - have had masks made for cosmic rays
     - have been corrected for nonlinearity
     - have been converted from DN to e-
-    - have been desmeared if desmearing is appropriate.  Under normal
-    circumstances, darks should not be desmeared.  The only time desmearing
-    would be useful is in the unexpected case that, for example,
-    dark current is so high that it stands far above other noise that is
-    not smeared upon readout, such as clock-induced charge, 
-    fixed-pattern noise, and read noise.
+    - have NOT been desmeared. Darks should not be desmeared.  The only component 
+    of dark frames that would be subject to a smearing effect is dark current 
+    since it linearly increases with time, so the extra row read time affects 
+    the dark current per pixel.  However, illuminated images
+    would also contain this smeared dark current, so dark subtraction should 
+    remove this smeared dark current (and then desmearing may be applied to the 
+    processed image if appropriate).  
 
     Also, add_photon_noise() should NOT have been applied to the frames in
     dataset.  And note that creation of the
@@ -329,7 +331,7 @@ def calibrate_darks_lsq(dataset, detector_params, detector_regions=None):
         output Dark's dq after assigning these pixels a flag value of 256.
         They should have large err values.
         The pixels that are masked for EVERY frame in all sub-stacks
-        but 4 (or less) are assigned a flag value of
+        but 3 (or less) are assigned a flag value of
         1, which falls under the category of "Bad pixel - unspecified reason".
         These pixels would have no reliability for dark subtraction.
 
@@ -415,7 +417,7 @@ def calibrate_darks_lsq(dataset, detector_params, detector_regions=None):
         output Dark's dq after assigning these pixels a flag value of 256.
         They should have large err values.
         The pixels that are masked for EVERY frame in all sub-stacks
-        but 4 (or less) are assigned a flag value of
+        but 3 (or less) are assigned a flag value of
         1, which falls under the category of "Bad pixel - unspecified reason".
         These pixels would have no reliability for dark subtraction.
     FPN_std_map : array-like (full frame)
@@ -522,7 +524,7 @@ def calibrate_darks_lsq(dataset, detector_params, detector_regions=None):
     # flag value of 256; unreliable pixels, large err
     output_dq = (unreliable_pix_map >= len(datasets)-3).astype(int)*256
     # flag value of 1 for those that are masked all the way through for all
-    # but 4 (or less) stacks; this overwrites the flag value of 256 that was assigned to
+    # but 3 (or less) stacks; this overwrites the flag value of 256 that was assigned to
     # these pixels in previous line
     unfittable_ind = np.where(unfittable_pix_map >= len(datasets)-3)
     output_dq[unfittable_ind] = 1
@@ -577,7 +579,7 @@ def calibrate_darks_lsq(dataset, detector_params, detector_regions=None):
     # input data error comes from .err arrays; could use this for error bars
     # in input data for weighted least squares, but we'll just simply get the
     # std error and add it in quadrature to least squares fit standard dev
-    stacks_err = np.sqrt(np.sum(mean_err_stack**2, axis=0))/len(mean_err_stack)
+    stacks_err = np.sqrt(np.sum(mean_err_stack**2, axis=0)/np.sqrt(len(mean_err_stack)))
 
     # matrix to be used for least squares and covariance matrix
     X = np.array([np.ones([len(EMgain_arr)]), EMgain_arr, EMgain_arr*exptime_arr]).T
@@ -865,5 +867,5 @@ def build_synthesized_dark(dataset, noisemaps, detector_regions=None, full_frame
 
         master_dark = Dark(md_data, prihdr, exthdr, input_data, md_noise, FDCdq,
                         errhdr)
-
+        
         return master_dark
\ No newline at end of file
diff --git a/corgidrp/data.py b/corgidrp/data.py
index 97adc970..f1afd7c9 100644
--- a/corgidrp/data.py
+++ b/corgidrp/data.py
@@ -104,6 +104,8 @@ def update_after_processing_step(self, history_entry, new_all_data=None, new_all
             if new_all_err.shape[-2:] != self.all_err.shape[-2:] or new_all_err.shape[0] != self.all_err.shape[0]:
                 raise ValueError("The shape of new_all_err is {0}, whereas we are expecting {1}".format(new_all_err.shape, self.all_err.shape))
             self.all_err = new_all_err
+            for i in range(len(self.frames)):
+                self.frames[i].err = self.all_err[i]
         if new_all_dq is not None:
             if new_all_dq.shape != self.all_dq.shape:
                 raise ValueError("The shape of new_all_dq is {0}, whereas we are expecting {1}".format(new_all_dq.shape, self.all_dq.shape))
@@ -608,7 +610,7 @@ class Dark(Image):
         data_or_filepath (str or np.array): either the filepath to the FITS file to read in OR the 2D image data
         pri_hdr (astropy.io.fits.Header): the primary header (required only if raw 2D data is passed in)
         ext_hdr (astropy.io.fits.Header): the image extension header (required only if raw 2D data is passed in)
-        input_dataset (corgidrp.data.Dataset): the Image files combined together to make this noise map (required only if raw 2D data is passed in and if raw data filenames not already archived in ext_hdr)
+        input_dataset (corgidrp.data.Dataset): the Image files combined together to make this dark (required only if raw 2D data is passed in and if raw data filenames not already archived in ext_hdr)
         err (np.array): the error array (required only if raw data is passed in)
         err_hdr (astropy.io.fits.Header): the error header (required only if raw data is passed in)
         dq (np.array): the DQ array (required only if raw data is passed in)
@@ -625,7 +627,8 @@ def __init__(self, data_or_filepath, pri_hdr=None, ext_hdr=None, input_dataset=N
                 raise ValueError("This appears to be a new dark. The dataset of input files needs to be passed in to the input_dataset keyword to record history of this dark.")
             self.ext_hdr['DATATYPE'] = 'Dark' # corgidrp specific keyword for saving to disk
             self.ext_hdr['BUNIT'] = 'detected electrons'
-
+            if 'PC_STAT' not in ext_hdr:
+                self.ext_hdr['PC_STAT'] = 'analog master dark'
             # log all the data that went into making this calibration file
             if 'DRPNFILE' not in ext_hdr.keys() and input_dataset is not None:
                 self._record_parent_filenames(input_dataset)
@@ -636,8 +639,15 @@ def __init__(self, data_or_filepath, pri_hdr=None, ext_hdr=None, input_dataset=N
             # give it a default filename using the first input file as the base
             # strip off everything starting at .fits
             if input_dataset is not None:
-                orig_input_filename = input_dataset[0].filename.split(".fits")[0]
-                self.filename = "{0}_dark.fits".format(orig_input_filename)
+                if self.ext_hdr['PC_STAT'] != 'photon-counted master dark':
+                    orig_input_filename = input_dataset[0].filename.split(".fits")[0]
+                    self.filename = "{0}_dark.fits".format(orig_input_filename)
+                else:
+                    orig_input_filename = input_dataset[0].filename.split(".fits")[0]
+                    self.filename = "{0}_pc_dark.fits".format(orig_input_filename)
+        
+        if 'PC_STAT' not in self.ext_hdr:
+            self.ext_hdr['PC_STAT'] = 'analog master dark'
 
         if err_hdr is not None:
             self.err_hdr['BUNIT'] = 'detected electrons'
@@ -784,7 +794,8 @@ def __init__(self, data_or_filepath, pri_hdr=None, ext_hdr=None,
             raise ValueError("File that was loaded was not a NonLinearityCalibration file.")
         if self.ext_hdr['DATATYPE'] != 'NonLinearityCalibration':
             raise ValueError("File that was loaded was not a NonLinearityCalibration file.")
-
+        self.dq_hdr['COMMENT'] = 'DQ not meaningful for this calibration; just present for class consistency' 
+        
 class KGain(Image):
     """
     Class for KGain calibration file. Until further insights it is just one float value.
@@ -808,6 +819,13 @@ def __init__(self, data_or_filepath, err = None, ptc = None, pri_hdr=None, ext_h
        # run the image class contructor
         super().__init__(data_or_filepath, err=err, pri_hdr=pri_hdr, ext_hdr=ext_hdr, err_hdr=err_hdr)
 
+        # initialize these headers that have been recently added so that older calib files still contain this keyword when initialized and allow for tests that don't require 
+        # these values to run smoothly; if these values are actually required for 
+        # a particular process, the user would be alerted since these values below would result in an error as they aren't numerical
+        if 'RN' not in self.ext_hdr:
+            self.ext_hdr['RN'] = ''
+        if 'RN_ERR' not in self.ext_hdr:
+            self.ext_hdr['RN_ERR'] = ''
         # File format checks
         if self.data.shape != (1,1):
             raise ValueError('The KGain calibration data should be just one float value')
@@ -944,6 +962,8 @@ def __init__(self, data_or_filepath, pri_hdr=None, ext_hdr=None, input_dataset=N
             raise ValueError("File that was loaded was not a BadPixelMap file.")
         if self.ext_hdr['DATATYPE'] != 'BadPixelMap':
             raise ValueError("File that was loaded was not a BadPixelMap file.")
+        self.dq_hdr['COMMENT'] = 'DQ not meaningful for this calibration; just present for class consistency' 
+        self.err_hdr['COMMENT'] = 'err not meaningful for this calibration; just present for class consistency' 
 
 
 class DetectorNoiseMaps(Image):
@@ -1214,7 +1234,8 @@ def __init__(self, data_or_filepath, pri_hdr=None, ext_hdr=None, input_dataset=N
         # check that this is actually an AstrometricCalibration file that was read in
         if 'DATATYPE' not in self.ext_hdr or self.ext_hdr['DATATYPE'] != 'AstrometricCalibration':
             raise ValueError("File that was loaded was not an AstrometricCalibration file.")    
-        
+        self.dq_hdr['COMMENT'] = 'DQ not meaningful for this calibration; just present for class consistency'     
+    
 class TrapCalibration(Image):
     """
 
@@ -1261,6 +1282,7 @@ def __init__(self,data_or_filepath, pri_hdr=None,ext_hdr=None, input_dataset=Non
         # since if only a filepath was passed in, any file could have been read in
         if 'DATATYPE' not in self.ext_hdr or self.ext_hdr['DATATYPE'] != 'TrapCalibration':
             raise ValueError("File that was loaded was not a TrapCalibration file.")
+        self.dq_hdr['COMMENT'] = 'DQ not meaningful for this calibration; just present for class consistency' 
 
 class FluxcalFactor(Image):
     """
diff --git a/corgidrp/l2a_to_l2b.py b/corgidrp/l2a_to_l2b.py
index e294ef1b..2c0bd943 100644
--- a/corgidrp/l2a_to_l2b.py
+++ b/corgidrp/l2a_to_l2b.py
@@ -82,10 +82,14 @@ def dark_subtraction(input_dataset, dark, detector_regions=None, outputdir=None)
             outputdir = '.' #current directory
         dark.save(filedir=outputdir)
     elif type(dark) is data.Dark:
-       # In this case, the Dark loaded in should already match the arry dimensions
-       # of input_dataset, specified by full_frame argument of build_trad_dark
-       # when this Dark was built
-       pass
+        if 'PC_STAT' in dark.ext_hdr:
+            if dark.ext_hdr['PC_STAT'] != 'analog master dark':
+                raise Exception('The input \'dark\' is a photon-counted dark and cannot be used in the dark_subtraction step function, which is intended for analog frames.')
+        # In this case, the Dark loaded in should already match the arry dimensions
+        # of input_dataset, specified by full_frame argument of build_trad_dark
+        # when this Dark was built
+        if (dark.ext_hdr['EXPTIME'], dark.ext_hdr['CMDGAIN'], dark.ext_hdr['KGAIN']) != unique_vals[0]:
+            raise Exception('Dark should have the same EXPTIME, CMDGAIN, and KGAIN as input_dataset.')
     else:
         raise Exception('dark type should be either corgidrp.data.Dark or corgidrp.data.DetectorNoiseMaps.')
 
@@ -255,16 +259,20 @@ def convert_to_electrons(input_dataset, k_gain):
    # you should make a copy the dataset to start
     kgain_dataset = input_dataset.copy()
     kgain_cube = kgain_dataset.all_data
-    kgain_error = kgain_dataset.all_err
 
     kgain = k_gain.value #extract from caldb
+    kgainerr = k_gain.error[0]
+    # x = c*kgain, where c (counts beforehand) and kgain both have error, so do propogation of error due to the product of 2 independent sources
+    kgain_error = np.zeros_like(input_dataset.all_err)
+    kgain_tot_err = (np.sqrt((kgain*kgain_dataset.all_err[:,0,:,:])**2 + (kgain_cube*kgainerr)**2))
+    kgain_error[:,0,:,:] = kgain_tot_err
     kgain_cube *= kgain
-    kgain_error *= kgain
-
+    
     history_msg = "data converted to detected EM electrons by kgain {0}".format(str(kgain))
 
     # update the output dataset with this converted data and update the history
-    kgain_dataset.update_after_processing_step(history_msg, new_all_data=kgain_cube, new_all_err=kgain_error, header_entries = {"BUNIT":"detected EM electrons", "KGAIN":kgain})
+    kgain_dataset.update_after_processing_step(history_msg, new_all_data=kgain_cube, new_all_err=kgain_error, header_entries = {"BUNIT":"detected EM electrons", "KGAIN":kgain, 
+                                                                                    "KGAIN_ER": k_gain.error[0], "RN":k_gain.ext_hdr['RN'], "RN_ERR":k_gain.ext_hdr["RN_ERR"]})
     return kgain_dataset
 
 def em_gain_division(input_dataset):
diff --git a/corgidrp/mocks.py b/corgidrp/mocks.py
index 12d82613..546e2f6c 100644
--- a/corgidrp/mocks.py
+++ b/corgidrp/mocks.py
@@ -1,6 +1,7 @@
 import os
 from pathlib import Path
 import numpy as np
+import warnings
 import scipy.ndimage
 import pandas as pd
 import astropy.io.fits as fits
@@ -17,6 +18,7 @@
 import corgidrp.detector as detector
 from corgidrp.detector import imaging_area_geom, unpack_geom
 from corgidrp.pump_trap_calibration import (P1, P1_P1, P1_P2, P2, P2_P2, P3, P2_P3, P3_P3, tau_temp)
+from corgidrp.data import DetectorParams
 
 from emccd_detect.emccd_detect import EMCCDDetect
 from emccd_detect.util.read_metadata_wrapper import MetadataWrapper
@@ -167,7 +169,7 @@ def create_synthesized_master_dark_calib(detector_areas):
     # image area, including "shielded" rows and cols:
     imrows, imcols, imr0c0 = imaging_area_geom('SCI', detector_areas)
     prerows, precols, prer0c0 = unpack_geom('SCI', 'prescan', detector_areas)
-
+    
     frame_list = []
     for i in range(len(EMgain_arr)):
         for l in range(N): #number of frames to produce
@@ -235,7 +237,8 @@ def create_dark_calib_files(filedir=None, numfiles=10):
     for i in range(numfiles):
         prihdr, exthdr = create_default_headers()
         exthdr['KGAIN'] = 7
-        np.random.seed(456+i); sim_data = np.random.poisson(lam=150., size=(1200, 2200)).astype(np.float64)
+        #np.random.seed(456+i); 
+        sim_data = np.random.poisson(lam=150., size=(1200, 2200)).astype(np.float64)
         frame = data.Image(sim_data, pri_hdr=prihdr, ext_hdr=exthdr)
         if filedir is not None:
             frame.save(filedir=filedir, filename=filepattern.format(i))
@@ -264,7 +267,8 @@ def create_simflat_dataset(filedir=None, numfiles=10):
     for i in range(numfiles):
         prihdr, exthdr = create_default_headers()
         # generate images in normal distribution with mean 1 and std 0.01
-        np.random.seed(456+i); sim_data = np.random.poisson(lam=150., size=(1024, 1024)).astype(np.float64)
+        #np.random.seed(456+i); 
+        sim_data = np.random.poisson(lam=150., size=(1024, 1024)).astype(np.float64)
         frame = data.Image(sim_data, pri_hdr=prihdr, ext_hdr=exthdr)
         if filedir is not None:
             frame.save(filedir=filedir, filename=filepattern.format(i))
@@ -451,7 +455,8 @@ def create_flatfield_dummy(filedir=None, numfiles=2):
     frames=[]
     for i in range(numfiles):
         prihdr, exthdr = create_default_headers()
-        np.random.seed(456+i); sim_data = np.random.normal(loc=1.0, scale=0.01, size=(1024, 1024))
+        #np.random.seed(456+i); 
+        sim_data = np.random.normal(loc=1.0, scale=0.01, size=(1024, 1024))
         frame = data.Image(sim_data, pri_hdr=prihdr, ext_hdr=exthdr)
         if filedir is not None:
             frame.save(filedir=filedir, filename=filepattern.format(i))
@@ -490,7 +495,8 @@ def create_nonlinear_dataset(nonlin_filepath, filedir=None, numfiles=2,em_gain=2
         data_range = np.linspace(800,65536,size)
         # Generate data for each row, where the mean increase from 10 to 65536
         for x in range(size):
-            np.random.seed(120+x); sim_data[:, x] = np.random.poisson(data_range[x], size).astype(np.float64)
+            #np.random.seed(120+x); 
+            sim_data[:, x] = np.random.poisson(data_range[x], size).astype(np.float64)
 
         non_linearity_correction = data.NonLinearityCalibration(nonlin_filepath)
 
@@ -543,7 +549,7 @@ def create_cr_dataset(nonlin_filepath, filedir=None, datetime=None, numfiles=2,
     im_width = dataset.all_data.shape[-1]
 
     # Overwrite dataset with a poisson distribution
-    np.random.seed(123)
+    #np.random.seed(123)
     dataset.all_data[:,:,:] = np.random.poisson(lam=150,size=dataset.all_data.shape).astype(np.float64)
 
     # Loop over images in dataset
@@ -554,7 +560,7 @@ def create_cr_dataset(nonlin_filepath, filedir=None, datetime=None, numfiles=2,
 
         # Pick random locations to add a cosmic ray
         for x in range(numCRs):
-            np.random.seed(123+x)
+            #np.random.seed(123+x)
             loc = np.round(np.random.uniform(0,im_width-1, size=2)).astype(int)
 
             # Add the CR plateau
@@ -1859,6 +1865,116 @@ def add_defect(sch_imgs, prob, ori, temp):
                 else:
                     hdul.writeto(filename, overwrite = True)
 
+def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, exptime=0.05, cosmic_rate=0, full_frame=True, smear=True, flux=1):
+    '''This creates mock L1 Dataset containing frames with large gain and short exposure time, illuminated and dark frames.
+    Used for unit tests for photon counting.  
+    
+    Args:
+        Nbrights (int):  number of illuminated frames to simulate
+        Ndarks (int):  number of dark frames to simulate
+        EMgain (float): EM gain
+        kgain (float): k gain (e-/DN)
+        exptime (float): exposure time (in s)
+        cosmic_rate: (float) simulated cosmic rays incidence, hits/cm^2/s
+        full_frame: (bool) If True, simulated frames are SCI full frames.  If False, 50x50 images are simulated.  Defaults to True.
+        smear: (bool) If True, smear is simulated.  Defaults to True.
+        flux: (float) Number of photons/s per pixel desired.  Defaults to 1.
+    
+    Returns:
+        ill_dataset (corgidrp.data.Dataset): Dataset containing the illuminated frames
+        dark_dataset (corgidrp.data.Dataset): Dataset containing the dark frames
+        ill_mean (float): mean electron count value simulated in the illuminated frames
+        dark_mean (float): mean electron count value simulated in the dark frames
+    '''
+    pix_row = 1024 #number of rows and number of columns
+    fluxmap = flux*np.ones((pix_row,pix_row)) #photon flux map, photons/s
+
+    emccd = EMCCDDetect(
+        em_gain=EMgain,
+        full_well_image=60000.,  # e-
+        full_well_serial=100000.,  # e-
+        dark_current=8.33e-4,  # e-/pix/s
+        cic=0.01,  # e-/pix/frame
+        read_noise=100.,  # e-/pix/frame
+        bias=20000,  # e-
+        qe=0.9,  # quantum efficiency, e-/photon
+        cr_rate=cosmic_rate,  # cosmic rays incidence, hits/cm^2/s
+        pixel_pitch=13e-6,  # m
+        eperdn=kgain,  
+        nbits=64, # number of ADU bits
+        numel_gain_register=604 #number of gain register elements
+        )
+
+    thresh = emccd.em_gain/10 # threshold
+
+    if np.average(exptime*fluxmap) > 0.1:
+        warnings.warn('average # of photons/pixel is > 0.1.  Decrease frame '
+        'time to get lower average # of photons/pixel.')
+
+    if emccd.read_noise <=0:
+        warnings.warn('read noise should be greater than 0 for effective '
+        'photon counting')
+    if thresh < 4*emccd.read_noise:
+        warnings.warn('thresh should be at least 4 or 5 times read_noise for '
+        'accurate photon counting')
+
+    avg_ph_flux = np.mean(fluxmap)
+    # theoretical electron flux for brights
+    ill_mean = avg_ph_flux*emccd.qe*exptime + emccd.dark_current*exptime + emccd.cic
+    # theoretical electron flux for darks
+    dark_mean = emccd.dark_current*exptime + emccd.cic
+
+    if smear:
+        #simulate smear to fluxmap
+        detector_params = DetectorParams({})
+        rowreadtime = detector_params.params['rowreadtime']
+        smear = np.zeros_like(fluxmap)
+        m = len(smear)
+        for r in range(m):
+            columnsum = 0
+            for i in range(r+1):
+                columnsum = columnsum + rowreadtime*fluxmap[i,:]
+            smear[r,:] = columnsum
+        
+        fluxmap = fluxmap + smear/exptime
+    
+    frame_e_list = []
+    frame_e_dark_list = []
+    prihdr, exthdr = create_default_headers()
+    for i in range(Nbrights):
+        # Simulate bright
+        if full_frame:
+            frame_dn = emccd.sim_full_frame(fluxmap, exptime)
+        else:
+            frame_dn = emccd.sim_sub_frame(fluxmap[:50,:50], exptime)
+        frame = data.Image(frame_dn, pri_hdr=prihdr, ext_hdr=exthdr)
+        frame.ext_hdr['CMDGAIN'] = EMgain
+        frame.ext_hdr['EXPTIME'] = exptime
+        frame.ext_hdr['KGAIN'] = kgain
+        frame.ext_hdr['ISPC'] = True
+        frame.pri_hdr["VISTYPE"] = "TDEMO"
+        frame.filename = 'L1_for_pc_ill_{0}.fits'.format(i)
+        frame_e_list.append(frame)
+
+    for i in range(Ndarks):
+        # Simulate dark
+        if full_frame:
+            frame_dn_dark = emccd.sim_full_frame(np.zeros_like(fluxmap), exptime)
+        else:
+            frame_dn_dark = emccd.sim_sub_frame(np.zeros_like(fluxmap[:50,:50]), exptime)
+        frame_dark = data.Image(frame_dn_dark, pri_hdr=prihdr.copy(), ext_hdr=exthdr.copy())
+        frame_dark.ext_hdr['CMDGAIN'] = EMgain
+        frame_dark.ext_hdr['EXPTIME'] = exptime
+        frame_dark.ext_hdr['KGAIN'] = kgain
+        frame_dark.ext_hdr['ISPC'] = True
+        frame_dark.pri_hdr["VISTYPE"] = "DARK"
+        frame.filename = 'L1_for_pc_dark_{0}.fits'.format(i)
+        frame_e_dark_list.append(frame_dark)
+
+    ill_dataset = data.Dataset(frame_e_list)
+    dark_dataset = data.Dataset(frame_e_dark_list)
+
+    return ill_dataset, dark_dataset, ill_mean, dark_mean
 
 def create_flux_image(star_flux, fwhm, cal_factor, filedir=None, color_cor = 1., platescale=21.8, add_gauss_noise=True, noise_scale=1., background = 0., file_save=False):
     """
@@ -1965,4 +2081,4 @@ def create_flux_image(star_flux, fwhm, cal_factor, filedir=None, color_cor = 1.,
     if filedir is not None and file_save:
         frame.save(filedir=filedir, filename=filename)
 
-    return frame
\ No newline at end of file
+    return frame
diff --git a/corgidrp/photon_counting.py b/corgidrp/photon_counting.py
new file mode 100644
index 00000000..132b556e
--- /dev/null
+++ b/corgidrp/photon_counting.py
@@ -0,0 +1,438 @@
+"""Photon count a stack of analog images and return a mean expected electron 
+count array with photometric corrections.
+
+"""
+import warnings
+import numpy as np
+from astropy.io import fits
+import corgidrp.data as data
+
+def varL23(g, L, T, N):
+    '''Expected variance after photometric correction.  
+    See https://doi.org/10.1117/1.JATIS.9.4.048006 for details.
+    
+    Args:
+        g (scalar): EM gain
+        L (2-D array): mean expected number of electrons
+        T (scalar): threshold 
+        N (2-D array): number of frames
+
+    Returns:
+        (float): variance from photon counting and the photometric correction
+
+    '''
+    Const = 6/(6 + L*(6 + L*(3 + L)))
+    eThresh = (Const*(np.e**(-T/g)*L*(2*g**2*(6 + L*(3 + L)) +
+            2*g*L*(3 + L)*T + L**2*T**2))/(12*g**2))
+    std_dev = np.sqrt(N * eThresh * (1-eThresh))
+
+    return (std_dev)**2*(((np.e**((T/g)))/N) +
+    2*((np.e**((2*T)/g)*(g - T))/(2*g*N**2))*(N*eThresh) +
+    3*(((np.e**(((3*T)/g)))*(4*g**2 - 8*g*T + 5*T**2))/(
+    12*g**2*N**3))*(N*eThresh)**2)**2
+
+class PhotonCountException(Exception):
+    """Exception class for photon_count module."""
+
+
+def photon_count(e_image, thresh):
+    """Convert an analog image into a photon-counted image.
+
+    Args: 
+        e_image (array_like, float): Analog image (e-).
+        thresh (float): Photon counting threshold (e-). Values > thresh will be assigned a 1,
+            values <= thresh will be assigned a 0.
+
+    Returns:
+        pc_image (array_like, float): Output digital image in units of photons.
+
+    B Nemati and S Miller - UAH - 06-Aug-2018
+
+    """
+    # Check if input is an array/castable to one
+    e_image = np.array(e_image).astype(float)
+    if len(e_image.shape) == 0:
+        raise PhotonCountException('e_image must have length > 0')
+
+    pc_image = np.zeros(e_image.shape, dtype=int)
+    pc_image[e_image > thresh] = 1
+
+    return pc_image
+
+
+def get_pc_mean(input_dataset, pc_master_dark=None, T_factor=None, pc_ecount_max=None, niter=2, mask_filepath=None, safemode=True, inputmode='illuminated'):
+    """Take a stack of images, frames of the same exposure 
+    time, k gain, read noise, and EM gain, and return the mean expected value per 
+    pixel. The frames are taken in photon-counting mode, which means high EM 
+    gain and very low exposure time.  
+
+    The frames in each stack should be SCI image (1024x1024) frames that have 
+    had some of the L2b steps applied:
+
+    - have had their bias subtracted
+    - have had masks made for cosmic rays
+    - have been corrected for nonlinearity (These first 3 steps make the frames L2a.)
+    - have been frame-selected (to weed out bad frames)
+    - have been converted from DN to e-
+
+    This algorithm will photon count each frame individually,
+    then co-add the photon-counted frames. The co-added frame is then averaged
+    and corrected for thresholding and coincidence loss, returning the mean
+    expected array in units of photoelectrons if dark-subtracted (detected electrons 
+    if not dark-subtracted).  The threshold is determined by the input "T_factor", 
+    and the value stored in DetectorParams is used if this input is None. 
+    
+    This function can be used for photon-counting illuminated and dark datasets 
+    (see Args below). 
+
+    Args:
+        input_dataset (corgidrp.data.Dataset): This is an instance of corgidrp.data.Dataset containing the 
+            frames to photon-count. All the frames must have the same 
+            exposure time, EM gain, k gain, and read noise.  If the input dataset's header key 'VISTYPE' is equal to 'DARK', a 
+            photon-counted master dark calibration product will be produced.  Otherwise, the input dataset is assumed to consist of illuminated 
+            frames intended for photon-counting.
+        pc_master_dark (corgidrp.data.Dark): Photon-counted master dark to be used for dark subtraction.  
+            If None, no dark subtraction is done.
+        T_factor (float): The number of read noise standard deviations at which to set the threshold for photon counting.  If None, the value is drawn from corgidrp.data.DetectorParams. 
+            Defaults to None.
+        pc_ecount_max (float): Maximum allowed electrons/pixel/frame for photon counting.  If None, the value is drawn from corgidrp.data.DetectorParams. 
+            Defaults to None.
+        niter (int, optional): Number of Newton's method iterations (used for photometric 
+            correction). Defaults to 2.
+        mask_filepath (str): Filepath to a .fits file with a pixel mask in the default HDU.  The mask should be an array of the same shape as that found in each frame of the input_dataset, and the mask 
+            should be 0 where the region of interest is and 1 for pixels to be masked (not considered).
+        safemode (bool): If False, the function does not halt due to an exception (useful for the iterative process of digging the dark hole).  
+            If True, the function halts with an exception if the mean intensity of the unmasked pixels (or all pixels if no mask provided) is greater than pc_ecount_max or if the minimum photon-counted pixel value is negative.
+            If False, the function gives a warning for these instead.
+            Defaults to True.
+        inputmode (str):  If 'illuminated', the frames are assumed to be illuminated frames.  If 'darks', frames are assumed to be dark frames input for creation of a photon-counted master dark. 
+            This flag shows the user's intention with the input, and this input is checked against the file type of the dataset for compatibility (e.g., if this input is 'darks' while 'VISTYPE' is not equal 
+            to 'DARK', an exception is raised).
+
+    Returns:
+        corgidrp.data.Dataset or corgidrp.data.Dark: If If the input dataset's header key 'VISTYPE' is not equal to 'DARK', 
+            corgidrp.data.Dataset is the output type, and the output is the processed illuminated set, whether 
+            dark subtraction happened or not.  Contains mean expected array (detected electrons if not dark-subtracted, 
+            photoelectrons if dark-subtracted). 
+            If the input dataset's header key 'VISTYPE' is equal to 'DARK', corgidrp.data.Dark is the output type, and the output
+            is the processed dark set.  Contains mean expected array (detected electrons).
+
+    References
+    ----------
+    [1] https://www.spiedigitallibrary.org/conference-proceedings-of-spie/11443/114435F/Photon-counting-and-precision-photometry-for-the-Roman-Space-Telescope/10.1117/12.2575983.full
+    [2] https://doi.org/10.1117/1.JATIS.9.4.048006
+
+    B Nemati, S Miller - UAH - 13-Dec-2020
+    Kevin Ludwick - UAH - 2023
+
+    """
+    if not isinstance(niter, (int, np.integer)) or niter < 1:
+            raise PhotonCountException('niter must be an integer greater than '
+                                        '0')
+    test_dataset, _ = input_dataset.copy().split_dataset(exthdr_keywords=['EXPTIME', 'CMDGAIN', 'KGAIN', 'RN'])
+    if len(test_dataset) > 1:
+        raise PhotonCountException('All frames must have the same exposure time, '
+                                   'commanded EM gain, and k gain.')
+    datasets, val = test_dataset[0].split_dataset(prihdr_keywords=['VISTYPE'])
+    if len(val) != 1:
+        raise PhotonCountException('There should only be 1 \'VISTYPE\' value for the dataset.')
+    if val[0] == 'DARK':
+        if inputmode != 'darks':
+            raise PhotonCountException('Inputmode is not \'darks\', but the input dataset has \'VISTYPE\' = \'DARK\'.')
+        if pc_master_dark is not None:
+            raise PhotonCountException('The input frames are \'VISTYPE\'=\'DARK\' frames, so no pc_master_dark should be provided.')
+    if val[0] != 'DARK':
+        if inputmode != 'illuminated':
+            raise PhotonCountException('Inputmode is not \'illuminated\', but the input dataset has \'VISTYPE\' not equal to \'DARK\'.')
+    if 'ISPC' in datasets[0].frames[0].ext_hdr:
+        if datasets[0].frames[0].ext_hdr['ISPC'] != True:
+            raise PhotonCountException('\'ISPC\' header value must be True if these frames are to be photon-counted.')
+
+    dataset = datasets[0]
+    
+    pc_means = []
+    errs = []
+    dqs = []
+    # getting number of read noise standard deviations at which to set the
+    # photon-counting threshold
+    if T_factor is None:
+        detector_params = data.DetectorParams({})
+        T_factor = detector_params.params['T_factor']
+    # getting maximum allowed electrons/pixel/frame for photon counting
+    if pc_ecount_max is None:
+        detector_params = data.DetectorParams({})
+        pc_ecount_max = detector_params.params['pc_ecount_max']
+    if mask_filepath is None:
+        mask = np.zeros_like(dataset.frames[0].data)
+    else:
+        mask = fits.getdata(mask_filepath)
+    
+    considered_indices = (mask==0)
+
+    # now get threshold to use for photon-counting
+    read_noise = test_dataset[0].frames[0].ext_hdr['RN']
+    thresh = T_factor*read_noise
+    if thresh < 0:
+        raise PhotonCountException('thresh must be nonnegative')
+    
+    if 'EMGAIN_M' in dataset.frames[0].ext_hdr: # if EM gain measured directly from frame TODO change hdr name if necessary
+        em_gain = dataset.frames[0].ext_hdr['EMGAIN_M']
+    elif 'EMGAIN_A' in dataset.frames[0].ext_hdr: # use applied EM gain if available
+        em_gain = dataset.frames[0].ext_hdr['EMGAIN_A']
+    elif 'CMDGAIN' in dataset.frames[0].ext_hdr: # use commanded gain otherwise
+        em_gain = dataset.frames[0].ext_hdr['CMDGAIN']
+
+    if thresh >= em_gain:
+        if safemode:
+            raise Exception('thresh should be less than em_gain for effective '
+            'photon counting')
+        else:
+            warnings.warn('thresh should be less than em_gain for effective '
+            'photon counting')
+    if np.nanmean(dataset.all_data[:, considered_indices]/em_gain) > pc_ecount_max:
+        if safemode:
+            raise Exception('average # of electrons/pixel is > pc_ecount_max, which means '
+            'the average # of photons/pixel may be > pc_ecount_max (depending on the QE).  Can decrease frame '
+            'time to get lower average # of photons/pixel.')
+        else:
+            warnings.warn('average # of electrons/pixel is > pc_ecount_max, which means '
+            'the average # of photons/pixel may be > pc_ecount_max (depending on the QE).  Can decrease frame '
+            'time to get lower average # of photons/pixel.')
+    if read_noise <=0:
+        raise Exception('read noise should be greater than 0 for effective '
+        'photon counting')
+    if thresh < 4*read_noise: # leave as just a warning
+        warnings.warn('thresh should be at least 4 or 5 times read_noise for '
+        'accurate photon counting')
+
+    # Photon count stack of frames
+    frames_pc = photon_count(dataset.all_data, thresh)
+    bool_map = dataset.all_dq.astype(bool).astype(float)
+    bool_map[bool_map > 0] = np.nan
+    bool_map[bool_map == 0] = 1
+    nframes = np.nansum(bool_map, axis=0)
+    # upper and lower bounds for PC (for getting accurate err)
+    frames_pc_up = photon_count(dataset.all_data+dataset.all_err[:,0], thresh)
+    frames_pc_low = photon_count(dataset.all_data-dataset.all_err[:,0], thresh)
+    frames_pc_masked = frames_pc * bool_map
+    frames_pc_masked_up = frames_pc_up * bool_map
+    frames_pc_masked_low = frames_pc_low * bool_map
+    # Co-add frames
+    frame_pc_coadded = np.nansum(frames_pc_masked, axis=0)
+    frame_pc_coadded_up = np.nansum(frames_pc_masked_up, axis=0)
+    frame_pc_coadded_low = np.nansum(frames_pc_masked_low, axis=0)
+    
+    # Correct for thresholding and coincidence loss; any pixel masked all the 
+    # way through the stack may give NaN, but it should be masked in lam_newton_fit(); 
+    # and it doesn't matter anyways since its DQ value will be 1 (it will be masked when the 
+    # bad pixel correction is run, which comes after this photon-counting step)
+    mean_expected = corr_photon_count(frame_pc_coadded, nframes, thresh,
+                                        em_gain, considered_indices, niter)
+    mean_expected_up = corr_photon_count(frame_pc_coadded_up, nframes, thresh,
+                                        em_gain, considered_indices, niter)
+    mean_expected_low = corr_photon_count(frame_pc_coadded_low, nframes, thresh,
+                                        em_gain, considered_indices, niter)
+    ##### error calculation: accounts for err coming from input dataset and 
+    # statistical error from the photon-counting and photometric correction process. 
+    # expected error from photon counting (biggest source from the actual values, not 
+    # mean_expected_up or mean_expected_low):
+    pc_variance = varL23(em_gain,mean_expected,thresh,nframes)
+    up = mean_expected_up +  pc_variance
+    low = mean_expected_low -  pc_variance
+    errs.append(np.max([up - mean_expected, mean_expected - low], axis=0))
+    dq = (nframes == 0).astype(int) 
+    pc_means.append(mean_expected)
+    dqs.append(dq)
+    
+    if pc_master_dark is not None:
+        if type(pc_master_dark) is not data.Dark:
+            raise Exception('Input type for pc_master_dark must be corgidrp.data.Dark.')
+        if 'PC_STAT' not in pc_master_dark.ext_hdr:
+            raise PhotonCountException('\'PC_STAT\' must be a key in the extension header of pc_master_dark.')
+        if pc_master_dark.ext_hdr['PC_STAT'] != 'photon-counted master dark':
+            raise PhotonCountException('pc_master_dark must be a photon-counted master dark (i.e., '
+                                       'the extension header key \'PC_STAT\' must be \'photon-counted master dark\').')
+        if 'PCTHRESH' not in pc_master_dark.ext_hdr:
+            raise PhotonCountException('Threshold should be stored under the header \'PCTHRESH\'.')
+        if pc_master_dark.ext_hdr['PCTHRESH'] != thresh:
+            raise PhotonCountException('Threshold used for photon-counted master dark should match the threshold to be used for the illuminated frames.')
+        pc_means.append(pc_master_dark.data)
+        dqs.append(pc_master_dark.dq)
+        errs.append(pc_master_dark.err)
+        dark_sub = "yes"
+        filename_end = ''
+    else:
+        pc_means.append(np.zeros_like(pc_means[0]))
+        dqs.append(np.zeros_like(pc_means[0]).astype(int))
+        errs.append(np.zeros_like(pc_means[0]))
+        dark_sub = "no"
+        filename_end = '_no_ds'
+
+    # now subtract the dark PC mean
+    combined_pc_mean = pc_means[0] - pc_means[1]
+    combined_pc_mean[combined_pc_mean<0] = 0
+    combined_err = np.sqrt(errs[0]**2 + errs[1]**2)
+    combined_dq = np.bitwise_or(dqs[0], dqs[1])
+    pri_hdr = dataset[0].pri_hdr.copy()
+    ext_hdr = dataset[0].ext_hdr.copy()
+    err_hdr = dataset[0].err_hdr.copy()
+    dq_hdr = dataset[0].dq_hdr.copy()
+    hdulist = dataset[0].hdu_list.copy()
+
+    if val[0] != "DARK":  
+        new_image = data.Image(combined_pc_mean, pri_hdr=pri_hdr, ext_hdr=ext_hdr, err=combined_err, dq=combined_dq, err_hdr=err_hdr, 
+                            dq_hdr=dq_hdr, input_hdulist=hdulist) 
+        new_image.filename = dataset[0].filename.split('.')[0]+'_pc'+filename_end+'.fits'
+        new_image.ext_hdr['PCTHRESH'] = thresh
+        new_image._record_parent_filenames(input_dataset) 
+        pc_ill_dataset = data.Dataset([new_image])
+        pc_ill_dataset.update_after_processing_step("Photon-counted {0} frames using T_factor={1} and niter={2}. Dark-subtracted: {3}.".format(len(input_dataset), T_factor, niter, dark_sub))
+        return pc_ill_dataset
+    else:
+        ext_hdr['PC_STAT'] = 'photon-counted master dark'
+        ext_hdr['NAXIS1'] = combined_pc_mean.shape[0]
+        ext_hdr['NAXIS2'] = combined_pc_mean.shape[1]
+        ext_hdr['PCTHRESH'] = thresh
+        pc_dark = data.Dark(combined_pc_mean, pri_hdr=pri_hdr, ext_hdr=ext_hdr, err=combined_err, dq=combined_dq, err_hdr=err_hdr, input_dataset=input_dataset)
+        return pc_dark
+
+
+def corr_photon_count(nobs, nfr, t, g, mask_indices, niter=2):
+    """Correct photon counted images.
+
+    Args:
+        nobs (array_like): Number of observations (Co-added photon-counted frame).
+        nfr (int): Number of coadded frames, accounting for masked pixels in the frames.
+        t (float): Photon-counting threshold.
+        g (float): EM gain.
+        mask_indices (array-like): indices of pixel positions to use.
+        niter (int, optional): Number of Newton's method iterations. Defaults to 2.
+
+    Returns:
+        lam (array_like): Mean expeted electrons per pixel (lambda).
+
+    """
+    # Get an approximate value of lambda for the first guess of the Newton fit
+    lam0 = calc_lam_approx(nobs, nfr, t, g)
+
+    # Use Newton's method to converge at a value for lambda
+    lam = lam_newton_fit(nobs, nfr, t, g, lam0, niter, mask_indices)
+
+    return lam
+
+
+def calc_lam_approx(nobs, nfr, t, g):
+    """Approximate lambda calculation.
+
+    This will calculate the first order approximation of lambda, and for values
+    that are out of bounds (e.g. from statistical fluctuations) it will revert
+    to the zeroth order.
+
+    Args:
+        nobs (array_like): Number of observations (Co-added photon counted frame).
+        nfr (int): Number of coadded frames.
+        t (float): Photon counting threshold.
+        g (float): EM gain used when taking images.
+
+    Returns:
+        lam1 (array_like): Mean expected (lambda).
+
+    """
+    # First step of equation (before taking log)
+    init = 1 - (nobs/nfr) * np.exp(t/g)
+    # Mask out all values less than or equal to 0
+    lam_m = np.zeros_like(init).astype(bool)
+    lam_m[init > 0] = True
+
+    # Use the first order approximation on all values greater than zero
+    lam1 = np.zeros_like(init)
+    lam1[lam_m] = -np.log(init[lam_m])
+
+    # For any values less than zero, revert to the zeroth order approximation
+    lam0 = nobs / nfr
+    lam1[~lam_m] = lam0[~lam_m]
+
+    return lam1
+
+
+def lam_newton_fit(nobs, nfr, t, g, lam0, niter, mask_indices):
+    """Newton fit for finding lambda.
+
+    Args:
+        nobs (array_like): Number of observations (Co-added photon counted frame).
+        nfr (int): Number of coadded frames.
+        t (float): Photon counting threshold.
+        g (float): EM gain used when taking images.
+        lam0 (array_like): Initial guess for lambda.
+        niter (int): Number of Newton's fit iterations to take.
+        mask_indices (array-like): indices of pixel positions to use.
+
+    Returns:
+        lam (array_like): Mean expected (lambda).
+
+    """
+    # Mask out zero values to avoid divide by zero
+    lam_est_m = np.ma.masked_array(lam0, mask=(lam0==0))
+    nobs_m = np.ma.masked_array(nobs, mask=(nobs==0))
+
+    # Iterate Newton's method
+    for i in range(niter):
+        func = _calc_func(nobs_m, nfr, t, g, lam_est_m)
+        dfunc = _calc_dfunc(nfr, t, g, lam_est_m)
+        lam_est_m -= func / dfunc
+
+    if np.nanmin(lam_est_m.data[mask_indices]) < 0:
+        raise PhotonCountException('negative number of photon counts; '
+        'try decreasing the frametime')
+
+    # Fill zero values back in
+    lam = lam_est_m.filled(0)
+
+    return lam
+
+def _calc_func(nobs, nfr, t, g, lam):
+    """Objective function for lambda for Newton's method for all applying photometric correction.
+    
+    Args:
+        nobs (array-like): number of frames per pixel that passed the threshold
+        nfr (array-like): number of unmasked frames per pixel total
+        t (float): threshold for photon counting
+        g (float): EM gain
+        lam (array-like): estimated mean expected electron count
+
+    Returns:
+        func (array-like): objective function
+
+    """
+    epsilon_prime = (lam*(2*g**2*(6 + lam*(3 + lam)) + 2*g*lam*(3 + lam)*t + 
+            lam**2*t**2))/(2.*np.e**(t/g)*g**2*(6 + lam*(6 + lam*(3 + lam))))
+
+    #if (nfr * epsilon_prime).any() > nobs.any():
+    #    warnings.warn('Input photon flux is too high; decrease frametime')
+    # This warning isn't necessary; could have a negative func but still
+    # close enough to 0 for Newton's method
+    func = nfr * epsilon_prime - nobs
+
+    return func
+
+
+def _calc_dfunc(nfr, t, g, lam):
+    """Derivative with respect to lambda of objective function.
+    
+    Args:
+        nfr (array-like): number of unmasked frames per pixel total
+        t (float): threshold for photon counting
+        g (float): EM gain
+        lam (array-like): estimated mean expected electron count
+    
+    Returns:
+        dfunc (array-like): derivative with respect to lambda of objective function
+    """
+    dfunc = (lam*nfr*(2*g**2*(3 + 2*lam) + 2*g*lam*t + 2*g*(3 + lam)*t + 
+            2*lam*t**2))/(2.*np.e**(t/g)*g**2*(6 + lam*(6 + lam*(3 + lam)))) - (lam*(6 + 
+            lam*(3 + lam) + lam*(3 + 2*lam))*nfr*
+         (2*g**2*(6 + lam*(3 + lam)) + 2*g*lam*(3 + lam)*t + lam**2*t**2))/(2.*np.e**(t/g)*g**2*(6 + 
+        lam*(6 + lam*(3 + lam)))**2) + (nfr*(2*g**2*(6 + lam*(3 + lam)) + 2*g*lam*(3 + lam)*t + 
+        lam**2*t**2))/(2.*np.e**(t/g)*g**2*(6 + lam*(6 + lam*(3 + lam))))
+
+    return dfunc
\ No newline at end of file
diff --git a/corgidrp/recipe_templates/l1_to_l2b_pc.json b/corgidrp/recipe_templates/l1_to_l2b_pc.json
new file mode 100644
index 00000000..7067704b
--- /dev/null
+++ b/corgidrp/recipe_templates/l1_to_l2b_pc.json
@@ -0,0 +1,82 @@
+{
+    "name" : "l1_to_l2b_pc",
+    "template" : true,
+    "drpconfig" : {
+        "track_individual_errors" : false
+    },
+    "inputs" : [],
+    "outputdir" : "",
+    "steps" : [
+        {
+            "name" : "prescan_biassub",
+            "calibs" : {
+                "DetectorNoiseMaps" : "AUTOMATIC"
+            },
+            "keywords" : {
+                "return_full_frame" : false
+            }
+        },
+        {
+            "name" : "detect_cosmic_rays",
+            "calibs" : {
+                "DetectorParams" : "AUTOMATIC",
+                "KGain" : "AUTOMATIC"
+            }
+        },
+        {
+            "name" : "correct_nonlinearity",
+            "calibs" : {
+                "NonLinearityCalibration" : "AUTOMATIC"
+            }
+        },
+        {
+            "name" : "update_to_l2a"
+        },
+
+        {
+            "name" : "frame_select"
+        },
+        {
+            "name" : "convert_to_electrons",
+            "calibs" : {
+                "KGain" : "AUTOMATIC"
+            }
+        },
+        {
+            "name": "get_pc_mean",
+            "calibs" : {
+                "Dark" : "AUTOMATIC"
+            }
+        },
+        {
+            "name" : "desmear",
+            "calibs" : {
+                "DetectorParams" : "AUTOMATIC"
+            }
+        },
+        {
+            "name" : "cti_correction",
+            "calibs" : {
+                "TrapCalibration" : "AUTOMATIC,OPTIONAL"
+            }
+        },
+        {
+            "name" : "flat_division",
+            "calibs" : {
+                "FlatField" : "AUTOMATIC"
+            }
+        },
+        {
+            "name" : "correct_bad_pixels",
+            "calibs" : {
+                "BadPixelMap" : "AUTOMATIC"
+            }
+        },
+        {
+            "name" : "update_to_l2b"
+        },
+        {
+            "name" : "save"
+        }
+    ]
+}
\ No newline at end of file
diff --git a/corgidrp/recipe_templates/l1_to_l2b_pc_dark.json b/corgidrp/recipe_templates/l1_to_l2b_pc_dark.json
new file mode 100644
index 00000000..02e97622
--- /dev/null
+++ b/corgidrp/recipe_templates/l1_to_l2b_pc_dark.json
@@ -0,0 +1,55 @@
+{
+    "name" : "l1_to_l2b_pc_dark",
+    "template" : true,
+    "drpconfig" : {
+        "track_individual_errors" : false
+    },
+    "inputs" : [],
+    "outputdir" : "",
+    "steps" : [
+        {
+            "name" : "prescan_biassub",
+            "calibs" : {
+                "DetectorNoiseMaps" : "AUTOMATIC"
+            },
+            "keywords" : {
+                "return_full_frame" : false
+            }
+        },
+        {
+            "name" : "detect_cosmic_rays",
+            "calibs" : {
+                "DetectorParams" : "AUTOMATIC",
+                "KGain" : "AUTOMATIC"
+            }
+        },
+        {
+            "name" : "correct_nonlinearity",
+            "calibs" : {
+                "NonLinearityCalibration" : "AUTOMATIC"
+            }
+        },
+        {
+            "name" : "update_to_l2a"
+        },
+
+        {
+            "name" : "frame_select"
+        },
+        {
+            "name" : "convert_to_electrons",
+            "calibs" : {
+                "KGain" : "AUTOMATIC"
+            }
+        },
+        {
+            "name": "get_pc_mean",
+            "keywords" : {
+                "inputmode" : "darks"
+            }
+        },
+        {
+            "name" : "save"
+        }
+    ]
+}
\ No newline at end of file
diff --git a/corgidrp/recipe_templates/l2a_build_trad_dark_image.json b/corgidrp/recipe_templates/l2a_build_trad_dark_image.json
new file mode 100644
index 00000000..7ba7146b
--- /dev/null
+++ b/corgidrp/recipe_templates/l2a_build_trad_dark_image.json
@@ -0,0 +1,26 @@
+{
+    "name" : "build_trad_dark_l2a",
+    "template" : true,
+    "drpconfig" : {
+        "track_individual_errors" : false
+    },
+    "inputs" : [],
+    "outputdir" : "",
+    "steps" : [
+        {
+            "name" : "convert_to_electrons",
+            "calibs" : {
+                "KGain" : "AUTOMATIC"
+            }
+        },
+        {
+            "name" : "build_trad_dark",
+            "calibs" : {
+                "DetectorParams" : "AUTOMATIC"
+            }
+        },
+        {
+            "name" : "save"
+        }
+    ]
+}
\ No newline at end of file
diff --git a/corgidrp/recipe_templates/l2a_to_l2b_pc.json b/corgidrp/recipe_templates/l2a_to_l2b_pc.json
new file mode 100644
index 00000000..0d69f670
--- /dev/null
+++ b/corgidrp/recipe_templates/l2a_to_l2b_pc.json
@@ -0,0 +1,56 @@
+{
+    "name" : "l2a_to_l2b_pc",
+    "template" : true,
+    "drpconfig" : {
+        "track_individual_errors" : false
+    },
+    "inputs" : [],
+    "outputdir" : "",
+    "steps" : [
+        {
+            "name" : "frame_select"
+        },
+        {
+            "name" : "convert_to_electrons",
+            "calibs" : {
+                "KGain" : "AUTOMATIC"
+            }
+        },
+        {
+            "name": "get_pc_mean",
+            "calibs" : {
+                "Dark" : "AUTOMATIC"
+            }
+        },
+        {
+            "name" : "desmear",
+            "calibs" : {
+                "DetectorParams" : "AUTOMATIC"
+            }
+        },
+        {
+            "name" : "cti_correction",
+            "calibs" : {
+                "TrapCalibration" : "AUTOMATIC,OPTIONAL"
+            }
+        },
+        {
+            "name" : "flat_division",
+            "calibs" : {
+                "FlatField" : "AUTOMATIC"
+            }
+        },
+        {
+            "name" : "correct_bad_pixels",
+            "calibs" : {
+                "BadPixelMap" : "AUTOMATIC"
+            }
+        },
+        {
+            "name" : "update_to_l2b"
+        },
+        {
+            "name" : "save"
+        }
+    ]
+}
\ No newline at end of file
diff --git a/corgidrp/recipe_templates/l2a_to_l2b_pc_dark.json b/corgidrp/recipe_templates/l2a_to_l2b_pc_dark.json
new file mode 100644
index 00000000..e53ac9d6
--- /dev/null
+++ b/corgidrp/recipe_templates/l2a_to_l2b_pc_dark.json
@@ -0,0 +1,29 @@
+{
+    "name" : "l2a_to_l2b_pc_dark",
+    "template" : true,
+    "drpconfig" : {
+        "track_individual_errors" : false
+    },
+    "inputs" : [],
+    "outputdir" : "",
+    "steps" : [
+        {
+            "name" : "frame_select"
+        },
+        {
+            "name" : "convert_to_electrons",
+            "calibs" : {
+                "KGain" : "AUTOMATIC"
+            }
+        },
+        {
+            "name": "get_pc_mean",
+            "keywords" : {
+                "inputmode" : "darks"
+            }
+        },
+        {
+            "name" : "save"
+        }
+    ]
+}
\ No newline at end of file
diff --git a/corgidrp/walker.py b/corgidrp/walker.py
index 25f9b441..6ca02469 100644
--- a/corgidrp/walker.py
+++ b/corgidrp/walker.py
@@ -11,6 +11,7 @@
 import corgidrp.caldb as caldb
 import corgidrp.l1_to_l2a
 import corgidrp.l2a_to_l2b
+import corgidrp.photon_counting
 import corgidrp.pump_trap_calibration
 import corgidrp.calibrate_nonlin
 import corgidrp.detector
@@ -41,7 +42,8 @@
     "create_onsky_flatfield" : corgidrp.detector.create_onsky_flatfield,
     "combine_subexposures" : corgidrp.combine.combine_subexposures,
     "build_trad_dark" : corgidrp.darks.build_trad_dark,
-    "sort_pupilimg_frames" : corgidrp.sorting.sort_pupilimg_frames
+    "sort_pupilimg_frames" : corgidrp.sorting.sort_pupilimg_frames,
+    "get_pc_mean" : corgidrp.photon_counting.get_pc_mean
 }
 
 recipe_dir = os.path.join(os.path.dirname(__file__), "recipe_templates")
@@ -235,16 +237,52 @@ def guess_template(dataset):
         elif image.pri_hdr['VISTYPE'] == "FFIELD":
             recipe_filename = "l1_flat_and_bp.json"
         elif image.pri_hdr['VISTYPE'] == "DARK":
-            recipe_filename = "l1_to_l2a_noisemap.json"
+            _, unique_vals = dataset.split_dataset(exthdr_keywords=['EXPTIME', 'CMDGAIN', 'KGAIN'])
+            if 'ISPC' in image.ext_hdr: # to be added eventually, so check if it's present first
+                if image.ext_hdr['ISPC']:
+                    recipe_filename = "l1_to_l2b_pc_dark.json"
+                elif len(unique_vals) > 1: # darks for noisemap creation
+                    recipe_filename = "l1_to_l2a_noisemap.json"
+                else: # then len(unique_vals) is 1 and not PC: traditional darks
+                    recipe_filename = "build_trad_dark_image.json"
+            else: # can't distinguish PC or analog then; default to analog
+                if len(unique_vals) > 1: # darks for noisemap creation
+                    recipe_filename = "l1_to_l2a_noisemap.json"
+                else: # then len(unique_vals) is 1: traditional darks
+                    recipe_filename = "build_trad_dark_image.json"
         elif image.pri_hdr['VISTYPE'] == "PUPILIMG":
             recipe_filename = ["l1_to_l2a_nonlin.json", "l1_to_kgain.json"]
         else:
-            recipe_filename = "l1_to_l2b.json"
+            if 'ISPC' in image.ext_hdr: 
+                if image.ext_hdr['ISPC']:
+                    recipe_filename = "l1_to_l2b_pc.json"
+                else:
+                    recipe_filename = "l1_to_l2b.json"   
+            else:  
+                recipe_filename = "l1_to_l2b.json"    
     elif image.ext_hdr['DATA_LEVEL'] == "L2a":
         if image.pri_hdr['VISTYPE'] == "DARK":
-            recipe_filename = "l2a_to_l2a_noisemap.json"
+            _, unique_vals = dataset.split_dataset(exthdr_keywords=['EXPTIME', 'CMDGAIN', 'KGAIN'])
+            if 'ISPC' in image.ext_hdr:
+                if image.ext_hdr['ISPC']:
+                    recipe_filename = "l2a_to_l2b_pc_dark.json"
+                elif len(unique_vals) > 1: # darks for noisemap creation
+                    recipe_filename = "l2a_to_l2a_noisemap.json"
+                else: # then len(unique_vals) is 1 and not PC: traditional darks
+                    recipe_filename = "l2a_build_trad_dark_image.json"
+            else: # can't distinguish PC or analog then; default to analog
+                if len(unique_vals) > 1: # darks for noisemap creation
+                    recipe_filename = "l2a_to_l2a_noisemap.json"
+                else: # then len(unique_vals) is 1: traditional darks
+                    recipe_filename = "l2a_build_trad_dark_image.json"
         else:
-            raise NotImplementedError()
+            if 'ISPC' in image.ext_hdr: 
+                if image.ext_hdr['ISPC']:
+                    recipe_filename = "l2a_to_l2b_pc.json"
+                else:
+                    raise NotImplementedError()
+            else:
+                raise NotImplementedError()
     else:
         raise NotImplementedError()
 
diff --git a/tests/e2e_tests/astrom_e2e.py b/tests/e2e_tests/astrom_e2e.py
index ff491ebd..facf1856 100644
--- a/tests/e2e_tests/astrom_e2e.py
+++ b/tests/e2e_tests/astrom_e2e.py
@@ -178,7 +178,7 @@ def test_astrom_e2e(tvacdata_path, e2eoutput_path):
     assert dec == pytest.approx(target[1], abs=8.333e-7)
 
 if __name__ == "__main__":
-    tvacdata_dir = "/Users/macuser/Roman/corgi_contributions/Callibration_Notebooks/TVAC"
+    tvacdata_dir = "/Users/kevinludwick/Library/CloudStorage/Box-Box/CGI_TVAC_Data/Working_Folder/" #"/Users/macuser/Roman/corgi_contributions/Callibration_Notebooks/TVAC"
     outputdir = thisfile_dir
 
     ap = argparse.ArgumentParser(description="run the l1->l2b->boresight end-to-end test")
diff --git a/tests/e2e_tests/flatfield_e2e.py b/tests/e2e_tests/flatfield_e2e.py
index dc27edd4..3c9885c2 100644
--- a/tests/e2e_tests/flatfield_e2e.py
+++ b/tests/e2e_tests/flatfield_e2e.py
@@ -120,11 +120,21 @@ def test_flat_creation_neptune(tvacdata_path, e2eoutput_path):
     this_caldb.create_entry(nonlinear_cal)
 
     # KGain
+    # remove other KGain calibrations that may exist in case they don't have the added header keywords
+    for i in range(len(this_caldb._db['Type'])):
+        if this_caldb._db['Type'][i] == 'KGain':
+            this_caldb._db = this_caldb._db.drop(i)
+        elif this_caldb._db['Type'][i] == 'FlatField':
+            this_caldb._db = this_caldb._db.drop(i)
     kgain_val = 8.7
+    # add in keywords not provided by create_default_headers() (since L1 headers are simulated from that function)
+    ext_hdr['RN'] = 100
+    ext_hdr['RN_ERR'] = 0
     kgain = data.KGain(np.array([[kgain_val]]), pri_hdr=pri_hdr, ext_hdr=ext_hdr, 
                     input_dataset=mock_input_dataset)
     kgain.save(filedir=flat_outputdir, filename="mock_kgain.fits")
-    this_caldb.create_entry(kgain)
+    this_caldb.create_entry(kgain, to_disk=False)
+    this_caldb.save()
 
     # NoiseMap
     with fits.open(fpn_path) as hdulist:
@@ -290,11 +300,21 @@ def test_flat_creation_uranus(tvacdata_path, e2eoutput_path):
     this_caldb.create_entry(nonlinear_cal)
 
     # KGain
+    # remove other KGain calibrations that may exist in case they don't have the added header keywords
+    for i in range(len(this_caldb._db['Type'])):
+        if this_caldb._db['Type'][i] == 'KGain':
+            this_caldb._db = this_caldb._db.drop(i)
+        elif this_caldb._db['Type'][i] == 'FlatField':
+            this_caldb._db = this_caldb._db.drop(i)
     kgain_val = 8.7
+    # add in keywords not provided by create_default_headers() (since L1 headers are simulated from that function)
+    ext_hdr['RN'] = 100
+    ext_hdr['RN_ERR'] = 0
     kgain = data.KGain(np.array([[kgain_val]]), pri_hdr=pri_hdr, ext_hdr=ext_hdr, 
                     input_dataset=mock_input_dataset)
     kgain.save(filedir=flat_outputdir, filename="mock_kgain.fits")
-    this_caldb.create_entry(kgain)
+    this_caldb.create_entry(kgain, to_disk=False)
+    this_caldb.save()
 
     # NoiseMap
     with fits.open(fpn_path) as hdulist:
diff --git a/tests/e2e_tests/noisemap_cal_e2e.py b/tests/e2e_tests/noisemap_cal_e2e.py
index 916c782c..e75682aa 100644
--- a/tests/e2e_tests/noisemap_cal_e2e.py
+++ b/tests/e2e_tests/noisemap_cal_e2e.py
@@ -129,6 +129,15 @@ def test_noisemap_calibration_from_l1(tvacdata_path, e2eoutput_path):
     mock_input_dataset = data.Dataset(mock_cal_filelist)
 
     this_caldb = caldb.CalDB() # connection to cal DB
+    # remove other KGain calibrations that may exist in case they don't have the added header keywords
+    for i in range(len(this_caldb._db['Type'])):
+        if this_caldb._db['Type'][i] == 'KGain':
+            this_caldb._db = this_caldb._db.drop(i)
+        elif this_caldb._db['Type'][i] == 'Dark':
+            this_caldb._db = this_caldb._db.drop(i)
+        elif this_caldb._db['Type'][i] == 'NonLinearityCalibration':
+            this_caldb._db = this_caldb._db.drop(i)
+    this_caldb.save()
 
     pri_hdr, ext_hdr = mocks.create_default_headers()
     ext_hdr["DRPCTIME"] = time.Time.now().isot
@@ -145,6 +154,9 @@ def test_noisemap_calibration_from_l1(tvacdata_path, e2eoutput_path):
     kgain_val = 8.7 # From TVAC-20 noise characterization measurements
     kgain = data.KGain(np.array([[kgain_val]]), pri_hdr=pri_hdr, ext_hdr=ext_hdr, 
                     input_dataset=mock_input_dataset)
+    # add in keywords that didn't make it into mock_kgain.fits, using values used in mocks.create_photon_countable_frames()
+    kgain.ext_hdr['RN'] = 100
+    kgain.ext_hdr['RN_ERR'] = 0
     kgain.save(filedir=noisemap_outputdir, filename="mock_kgain.fits")
     this_caldb.create_entry(kgain)
 
@@ -167,16 +179,15 @@ def test_noisemap_calibration_from_l1(tvacdata_path, e2eoutput_path):
     ##### Check against II&T ("TVAC") data
     corgidrp_noisemap_fname = os.path.join(noisemap_outputdir,output_filename)
     # iit_noisemap_fname = os.path.join(iit_noisemap_datadir,"iit_test_noisemaps.fits")
-
     corgidrp_noisemap = data.autoload(corgidrp_noisemap_fname)
-    this_caldb.remove_entry(corgidrp_noisemap)
-
+    
     assert(np.nanmax(np.abs(corgidrp_noisemap.data[0]- F_map)) < 1e-11)
     assert(np.nanmax(np.abs(corgidrp_noisemap.data[1]- C_map)) < 1e-11)
     assert(np.nanmax(np.abs(corgidrp_noisemap.data[2]- D_map)) < 1e-11)
     assert(np.abs(corgidrp_noisemap.ext_hdr['B_O']- bias_offset) < 1e-11)
     pass
 
+    this_caldb.remove_entry(corgidrp_noisemap)
     # for noise_ext in ["FPN_map","CIC_map","DC_map"]:
         # corgi_dat = detector.imaging_slice('SCI', corgidrp_noisemap.__dict__[noise_ext])
         # iit_dat = detector.imaging_slice('SCI', iit_noisemap.__dict__[noise_ext])
@@ -325,11 +336,20 @@ def test_noisemap_calibration_from_l2a(tvacdata_path, e2eoutput_path):
     mock_input_dataset = data.Dataset(mock_cal_filelist)
 
     this_caldb = caldb.CalDB() # connection to cal DB
-
+    # remove other KGain calibrations that may exist in case they don't have the added header keywords
+    for i in range(len(this_caldb._db['Type'])):
+        if this_caldb._db['Type'][i] == 'KGain':
+            this_caldb._db = this_caldb._db.drop(i)
+        elif this_caldb._db['Type'][i] == 'Dark':
+            this_caldb._db = this_caldb._db.drop(i)
+    this_caldb.save()
     # KGain calibration
     kgain_val = 8.7 # From TVAC-20 noise characterization measurements
     kgain = data.KGain(np.array([[kgain_val]]), pri_hdr=pri_hdr, ext_hdr=ext_hdr, 
                     input_dataset=mock_input_dataset)
+    # add in keywords that didn't make it into mock_kgain.fits, using values used in mocks.create_photon_countable_frames()
+    kgain.ext_hdr['RN'] = 100
+    kgain.ext_hdr['RN_ERR'] = 0
     kgain.save(filedir=noisemap_outputdir, filename="mock_kgain.fits")
     this_caldb.create_entry(kgain)
 
diff --git a/tests/e2e_tests/photon_count_e2e.py b/tests/e2e_tests/photon_count_e2e.py
new file mode 100644
index 00000000..f507f695
--- /dev/null
+++ b/tests/e2e_tests/photon_count_e2e.py
@@ -0,0 +1,180 @@
+# Photon Counting E2E Test Code
+
+import argparse
+import os
+import pytest
+import numpy as np
+import astropy.time as time
+from astropy.io import fits
+
+import corgidrp
+import corgidrp.data as data
+import corgidrp.mocks as mocks
+import corgidrp.walker as walker
+import corgidrp.caldb as caldb
+import corgidrp.detector as detector
+
+@pytest.mark.e2e
+def test_expected_results_e2e(tvacdata_path, e2eoutput_path):
+    processed_cal_path = os.path.join(tvacdata_path, "TV-36_Coronagraphic_Data", "Cals")
+    flat_path = os.path.join(processed_cal_path, "flat.fits")
+    bp_path = os.path.join(processed_cal_path, "bad_pix.fits")
+
+    np.random.seed(1234)
+    ill_dataset, dark_dataset, ill_mean, dark_mean = mocks.create_photon_countable_frames(Nbrights=160, Ndarks=161, cosmic_rate=1, flux=0.5)
+    output_dir = os.path.join(e2eoutput_path, 'pc_sim_test_data')
+    output_ill_dir = os.path.join(output_dir, 'ill_frames')
+    output_dark_dir = os.path.join(output_dir, 'dark_frames')
+    if not os.path.exists(output_dir):
+        os.mkdir(output_dir)
+    # empty out directory of any previous files
+    for f in os.listdir(output_dir):
+        if os.path.isdir(os.path.join(output_dir,f)):
+            continue
+        os.remove(os.path.join(output_dir,f))
+    if not os.path.exists(output_ill_dir):
+        os.mkdir(output_ill_dir)
+    if not os.path.exists(output_dark_dir):
+        os.mkdir(output_dark_dir)
+    # empty out directory of any previous files
+    for f in os.listdir(output_ill_dir):
+        os.remove(os.path.join(output_ill_dir,f))
+    for f in os.listdir(output_dark_dir):
+        os.remove(os.path.join(output_dark_dir,f))
+    ill_dataset.save(output_ill_dir, ['pc_frame_ill_{0}.fits'.format(i) for i in range(len(ill_dataset))])
+    dark_dataset.save(output_dark_dir, ['pc_frame_dark_{0}.fits'.format(i) for i in range(len(dark_dataset))])
+    l1_data_ill_filelist = []
+    l1_data_dark_filelist = []
+    for f in os.listdir(output_ill_dir):
+        l1_data_ill_filelist.append(os.path.join(output_ill_dir, f))
+    for f in os.listdir(output_dark_dir):
+        l1_data_dark_filelist.append(os.path.join(output_dark_dir, f))
+
+    this_caldb = caldb.CalDB() # connection to cal DB
+    # remove other KGain calibrations that may exist in case they don't have the added header keywords
+    for i in range(len(this_caldb._db['Type'])):
+        if this_caldb._db['Type'][i] == 'KGain':
+            this_caldb._db = this_caldb._db.drop(i)
+        elif this_caldb._db['Type'][i] == 'Dark':
+            this_caldb._db = this_caldb._db.drop(i)
+    this_caldb.save()
+
+    # KGain
+    kgain_val = 7 # default value used in mocks.create_photon_countable_frames()
+    pri_hdr, ext_hdr = mocks.create_default_headers()
+    ext_hdr["DRPCTIME"] = time.Time.now().isot
+    ext_hdr['DRPVERSN'] =  corgidrp.__version__
+    mock_input_dataset = data.Dataset(l1_data_ill_filelist)
+    kgain = data.KGain(np.array([[kgain_val]]), pri_hdr=pri_hdr, ext_hdr=ext_hdr, 
+                    input_dataset=mock_input_dataset)
+    # add in keywords that didn't make it into mock_kgain.fits, using values used in mocks.create_photon_countable_frames()
+    kgain.ext_hdr['RN'] = 100
+    kgain.ext_hdr['RN_ERR'] = 0
+    kgain.save(filedir=output_dir, filename="mock_kgain.fits")
+    this_caldb.create_entry(kgain)
+
+    # NoiseMap
+    noise_map_dat = np.zeros((3, detector.detector_areas['SCI']['frame_rows'], detector.detector_areas['SCI']['frame_cols']))
+    noise_map_noise = np.zeros([1,] + list(noise_map_dat.shape))
+    noise_map_dq = np.zeros(noise_map_dat.shape, dtype=int)
+    err_hdr = fits.Header()
+    err_hdr['BUNIT'] = 'detected electrons'
+    ext_hdr['B_O'] = 0
+    ext_hdr['B_O_ERR'] = 0
+    noise_map = data.DetectorNoiseMaps(noise_map_dat, pri_hdr=pri_hdr, ext_hdr=ext_hdr,
+                                    input_dataset=mock_input_dataset, err=noise_map_noise,
+                                    dq = noise_map_dq, err_hdr=err_hdr)
+    noise_map.save(filedir=output_dir, filename="mock_detnoisemaps.fits")
+    this_caldb.create_entry(noise_map)
+
+    corgidrp_dir = os.path.split(corgidrp.__path__[0])[0]
+    tests_dir = os.path.join(corgidrp_dir, 'tests')
+    # Nonlinearity calibration
+    ### Create a dummy non-linearity file ####
+    #Create a mock dataset because it is a required input when creating a NonLinearityCalibration
+    dummy_dataset = mocks.create_prescan_files()
+    # Make a non-linearity correction calibration file
+    input_non_linearity_filename = "nonlin_table_TVAC.txt"
+    input_non_linearity_path = os.path.join(tests_dir, "test_data", input_non_linearity_filename)
+    test_non_linearity_filename = input_non_linearity_filename.split(".")[0] + ".fits"
+    nonlin_fits_filepath = os.path.join(tests_dir, "test_data", test_non_linearity_filename)
+    tvac_nonlin_data = np.genfromtxt(input_non_linearity_path, delimiter=",")
+
+    pri_hdr, ext_hdr = mocks.create_default_headers()
+    new_nonlinearity = data.NonLinearityCalibration(tvac_nonlin_data,pri_hdr=pri_hdr,ext_hdr=ext_hdr,input_dataset = dummy_dataset)
+    new_nonlinearity.filename = nonlin_fits_filepath
+    new_nonlinearity.pri_hdr = pri_hdr
+    new_nonlinearity.ext_hdr = ext_hdr
+    this_caldb.create_entry(new_nonlinearity)
+
+    ## Flat field
+    with fits.open(flat_path) as hdulist:
+        flat_dat = hdulist[0].data
+    flat = data.FlatField(flat_dat, pri_hdr=pri_hdr, ext_hdr=ext_hdr, input_dataset=mock_input_dataset)
+    flat.save(filedir=output_dir, filename="mock_flat.fits")
+    this_caldb.create_entry(flat)
+
+    # bad pixel map
+    with fits.open(bp_path) as hdulist:
+        bp_dat = hdulist[0].data
+    bp_map = data.BadPixelMap(bp_dat, pri_hdr=pri_hdr, ext_hdr=ext_hdr, input_dataset=mock_input_dataset)
+    bp_map.save(filedir=output_dir, filename="mock_bpmap.fits")
+    this_caldb.create_entry(bp_map)
+
+    # make PC dark
+    # below I leave out the template specification to check that the walker recipe guesser works as expected
+    walker.walk_corgidrp(l1_data_dark_filelist, '', output_dir)#, template="l1_to_l2b_pc_dark.json")
+    for f in os.listdir(output_dir):
+        if f.endswith('_pc_dark.fits'):
+            pc_dark_filename = f
+    # calDB was just updated with the PC Dark that was created with the walker above
+    pc_dark_file = os.path.join(output_dir, pc_dark_filename)
+
+    # make PC illuminated, subtracting the PC dark
+    # below I leave out the template specification to check that the walker recipe guesser works as expected
+    walker.walk_corgidrp(l1_data_ill_filelist, '', output_dir)#, template="l1_to_l2b_pc.json")
+    # get photon-counted frame
+    for f in os.listdir(output_dir):
+        if f.endswith('_pc.fits'):
+            pc_filename = f
+    pc_file = os.path.join(output_dir, pc_filename)
+    pc_frame = fits.getdata(pc_file)
+    pc_frame_err = fits.getdata(pc_file, 'ERR')
+    pc_dark_frame = fits.getdata(pc_dark_file)
+    pc_dark_frame_err = fits.getdata(pc_dark_file, 'ERR')
+
+    # more frames gets a better agreement; agreement to 1% for ~160 darks and illuminated
+    assert np.isclose(np.nanmean(pc_frame), ill_mean - dark_mean, rtol=0.02) 
+    assert np.isclose(np.nanmean(pc_dark_frame), dark_mean, rtol=0.01) 
+    assert pc_frame_err.min() >= 0
+    assert pc_dark_frame_err.min() >= 0
+
+    # load in CalDB again to reflect the PC Dark that was implicitly added in (but not found in this_caldb, which was loaded before the Dark was created)
+    post_caldb = caldb.CalDB()
+    post_caldb.remove_entry(kgain)
+    post_caldb.remove_entry(noise_map)
+    post_caldb.remove_entry(new_nonlinearity)
+    post_caldb.remove_entry(flat)
+    post_caldb.remove_entry(bp_map)
+    pc_dark = data.Dark(pc_dark_file)
+    post_caldb.remove_entry(pc_dark)
+
+
+if __name__ == "__main__":
+    # Use arguments to run the test. Users can then write their own scripts
+    # that call this script with the correct arguments and they do not need
+    # to edit the file. The arguments use the variables in this file as their
+    # defaults allowing the user to edit the file if that is their preferred
+    # workflow.
+    thisfile_dir = os.path.dirname(__file__)
+    outputdir = thisfile_dir
+    tvacdata_dir =  "/Users/kevinludwick/Library/CloudStorage/Box-Box/CGI_TVAC_Data/Working_Folder/"
+
+    ap = argparse.ArgumentParser(description="run the l1->l2a end-to-end test")
+    ap.add_argument("-tvac", "--tvacdata_dir", default=tvacdata_dir,
+                    help="Path to CGI_TVAC_Data Folder [%(default)s]")
+    ap.add_argument("-o", "--outputdir", default=outputdir,
+                    help="directory to write results to [%(default)s]")
+    args = ap.parse_args()
+    outputdir = args.outputdir
+    test_expected_results_e2e(tvacdata_dir, outputdir)
diff --git a/tests/test_badpixelmap.py b/tests/test_badpixelmap.py
index 022cd05d..645b8c07 100644
--- a/tests/test_badpixelmap.py
+++ b/tests/test_badpixelmap.py
@@ -7,7 +7,7 @@
 from corgidrp.bad_pixel_calibration import create_bad_pixel_map
 from corgidrp.darks import build_trad_dark
 
-
+np.random.seed(456)
 def test_badpixelmap(): 
     '''
 
@@ -70,7 +70,7 @@ def test_badpixelmap():
             flat_frame.data[i_col, i_row] = 0.3
 
     ###### make the badpixel map (input the flat_dataset just as a dummy):
-    badpixelmap = create_bad_pixel_map(flat_dataset, dark_frame,flat_frame)
+    badpixelmap = create_bad_pixel_map(flat_dataset, dark_frame,flat_frame, dthresh=6)
     # Use np.unpackbits to unpack the bits - big endien
     badpixelmap_bits = np.unpackbits(badpixelmap.data[:, :, np.newaxis], axis=2)
 
diff --git a/tests/test_build_synthesized_dark.py b/tests/test_build_synthesized_dark.py
index 51834677..d1550f55 100644
--- a/tests/test_build_synthesized_dark.py
+++ b/tests/test_build_synthesized_dark.py
@@ -26,6 +26,7 @@
 C = noise_maps.CIC_map
 D = noise_maps.DC_map
 # just need some a dataset for input
+np.random.seed(456)
 dataset = create_dark_calib_files()
 # values used in create_dark_calib_files()
 g = 1
diff --git a/tests/test_caldb.py b/tests/test_caldb.py
index e8a99455..ba6556f7 100644
--- a/tests/test_caldb.py
+++ b/tests/test_caldb.py
@@ -1,4 +1,5 @@
 import os
+import numpy as np
 import pytest
 import corgidrp
 import corgidrp.caldb as caldb
@@ -16,6 +17,7 @@
 testcaldb_filepath = os.path.join(calibdir, "test_caldb.csv")
 
 # make some fake test data to use
+np.random.seed(456)
 dark_dataset = mocks.create_dark_calib_files()
 master_dark = data.Dark(dark_dataset[0].data, dark_dataset[0].pri_hdr, dark_dataset[0].ext_hdr, dark_dataset)
 # save master dark to disk to be loaded later
@@ -163,4 +165,8 @@ def test_caldb_scan():
     os.remove(testcaldb_filepath)
 
 if __name__ == "__main__":
-    test_get_calib()
\ No newline at end of file
+    test_get_calib()
+    test_caldb_create_default()
+    test_caldb_custom_filepath()
+    test_caldb_insert_and_remove()
+    test_caldb_scan()
\ No newline at end of file
diff --git a/tests/test_calibrate_darks_lsq.py b/tests/test_calibrate_darks_lsq.py
index 1c22310a..8b6d8fa3 100644
--- a/tests/test_calibrate_darks_lsq.py
+++ b/tests/test_calibrate_darks_lsq.py
@@ -11,6 +11,8 @@
 from corgidrp.mocks import detector_areas_test as dat
 from corgidrp.data import DetectorNoiseMaps, DetectorParams, Dataset
 
+# make test reproducible
+np.random.seed(4567)
 # use default parameters
 detector_params = DetectorParams({})
 # specified parameters simulated in simulated data from
@@ -207,6 +209,7 @@ def test_mean_num():
 
 
 if __name__ == '__main__':
+    test_mean_num()
     test_expected_results_sub()
     test_sub_stack_len()
     test_g_arr_unique()
@@ -214,7 +217,6 @@ def test_mean_num():
     test_t_arr_unique()
     test_t_gtr_0()
     test_k_gtr_0()
-    test_mean_num()
 
 
 
diff --git a/tests/test_combine.py b/tests/test_combine.py
index 306b7131..a04761d8 100644
--- a/tests/test_combine.py
+++ b/tests/test_combine.py
@@ -305,4 +305,6 @@ def test_invalid_collapse():
         combined_dataset = combine.combine_subexposures(dataset, collapse="invalid_option")
 
 if __name__ == "__main__":
-    test_mean_combine_subexposures()
\ No newline at end of file
+    test_mean_combine_subexposures()
+    test_mean_combine_subexposures_with_bad()
+    test_median_combine_subexposures()
\ No newline at end of file
diff --git a/tests/test_dark_sub.py b/tests/test_dark_sub.py
index 0d09546c..380c7d84 100644
--- a/tests/test_dark_sub.py
+++ b/tests/test_dark_sub.py
@@ -10,6 +10,8 @@
 import corgidrp.l2a_to_l2b as l2a_to_l2b
 from corgidrp.darks import build_trad_dark
 
+np.random.seed(456)
+
 old_err_tracking = corgidrp.track_individual_errors
 # use default parameters
 detector_params = data.DetectorParams({})
@@ -81,6 +83,7 @@ def test_dark_sub():
     # load in the dark
     dark_filepath = os.path.join(calibdir, dark_filename)
     new_dark = data.Dark(dark_filepath)
+    assert new_dark.ext_hdr['PC_STAT'] == 'analog master dark'
 
     # check the dark can be pickled (for CTC operations)
     pickled = pickle.dumps(new_dark)
@@ -91,6 +94,11 @@ def test_dark_sub():
     darkest_dataset = l2a_to_l2b.dark_subtraction(dark_dataset, new_dark, outputdir=calibdir)
     assert(dark_filename in str(darkest_dataset[0].ext_hdr["HISTORY"]))
 
+    # PC dark cannot be used for this step function
+    new_dark.ext_hdr['PC_STAT'] = 'photon-counted master dark'
+    with pytest.raises(Exception):
+        l2a_to_l2b.dark_subtraction(dark_dataset, new_dark, outputdir=calibdir)
+
     # check the level of the dataset is now approximately 0, leaving off telemetry row
     assert np.mean(darkest_dataset.all_data[:,:-1,:]) == pytest.approx(0, abs=1e-2)
 
diff --git a/tests/test_data_levels.py b/tests/test_data_levels.py
index 588070c9..f0de7b03 100644
--- a/tests/test_data_levels.py
+++ b/tests/test_data_levels.py
@@ -2,12 +2,15 @@
 Test the data level specification
 """
 import pytest
+import numpy as np
 import corgidrp.mocks as mocks
 import corgidrp.l1_to_l2a as l1_to_l2a
 import corgidrp.l2a_to_l2b as l2a_to_l2b
 import corgidrp.l2b_to_l3 as l2b_to_l3
 import corgidrp.l3_to_l4 as l3_to_l4
 
+np.random.seed(456)
+
 def test_l1_to_l4():
     """
     Tests a successful upgrade of L1 to L4 in bookkeeping only
diff --git a/tests/test_dataset.py b/tests/test_dataset.py
index 6124dc2c..8316dd16 100644
--- a/tests/test_dataset.py
+++ b/tests/test_dataset.py
@@ -8,6 +8,8 @@
 from corgidrp.data import Image, Dataset
 from corgidrp.mocks import create_default_headers, create_dark_calib_files
 
+np.random.seed(123)
+
 data = np.ones([1024,1024]) * 2
 err = np.zeros([1024,1024])
 err1 = np.ones([1024,1024])
diff --git a/tests/test_err_dq.py b/tests/test_err_dq.py
index 31c29828..7d2dd4a4 100644
--- a/tests/test_err_dq.py
+++ b/tests/test_err_dq.py
@@ -10,6 +10,7 @@
 import corgidrp.caldb as caldb
 from corgidrp.darks import build_trad_dark
 
+np.random.seed(123)
 
 data = np.ones([1024,1024]) * 2
 err = np.zeros([1024,1024])
diff --git a/tests/test_flat_div.py b/tests/test_flat_div.py
index 48b6c80e..9b94c9dc 100644
--- a/tests/test_flat_div.py
+++ b/tests/test_flat_div.py
@@ -10,6 +10,7 @@
 
 old_err_tracking = corgidrp.track_individual_errors
 
+np.random.seed(9292)
 def test_flat_div():
     """
     Generate mock input data and pass into flat division function
@@ -61,7 +62,7 @@ def test_flat_div():
     # perform checks after the flat divison
     assert(flat_filename in str(flatdivided_dataset[0].ext_hdr["HISTORY"]))
     # check the level of the dataset is now approximately 100
-    assert np.mean(flatdivided_dataset.all_data) == pytest.approx(150, abs=1e-2)
+    assert np.mean(flatdivided_dataset.all_data) == pytest.approx(150, abs=2e-2)
     # check the propagated errors
     assert flatdivided_dataset[0].err_hdr["Layer_2"] == "FlatField_error"
     print("mean of all simulated data",np.mean(simflat_dataset.all_data))
diff --git a/tests/test_frame_selection.py b/tests/test_frame_selection.py
index b6f41b15..89e49edd 100644
--- a/tests/test_frame_selection.py
+++ b/tests/test_frame_selection.py
@@ -1,7 +1,9 @@
 import pytest
+import numpy as np
 import corgidrp.mocks as mocks
 from corgidrp.l2a_to_l2b import frame_select
 
+np.random.seed(123)
 
 def test_no_selection():
     """
diff --git a/tests/test_kgain.py b/tests/test_kgain.py
index 990ada5a..9da40a3d 100644
--- a/tests/test_kgain.py
+++ b/tests/test_kgain.py
@@ -17,6 +17,8 @@ def test_kgain():
     dat = np.ones([1024,1024]) * 2
     err = np.ones([1,1024,1024]) * 0.5
     ptc = np.ones([2,1024])
+    exthd["RN"] = 100
+    exthd["RN_ERR"] = 2
     ptc_hdr = fits.Header()
     image1 = data.Image(dat,pri_hdr = prhd, ext_hdr = exthd, err = err)
     image2 = image1.copy()
@@ -84,6 +86,9 @@ def test_kgain():
     assert gain_dataset[0].ext_hdr["BUNIT"] == "detected EM electrons"
     assert gain_dataset[0].err_hdr["BUNIT"] == "detected EM electrons"
     assert gain_dataset[0].ext_hdr["KGAIN"] == k_gain
+    assert gain_dataset[0].ext_hdr["KGAIN_ER"] == kgain.error[0]
+    assert gain_dataset[0].ext_hdr["RN"] > 0
+    assert gain_dataset[0].ext_hdr["RN_ERR"] > 0
     assert("converted" in str(gain_dataset[0].ext_hdr["HISTORY"]))
 
 if __name__ == "__main__":
diff --git a/tests/test_kgain_cal.py b/tests/test_kgain_cal.py
index 4783724b..23ebfa8c 100644
--- a/tests/test_kgain_cal.py
+++ b/tests/test_kgain_cal.py
@@ -17,8 +17,10 @@
 from corgidrp import check
 from corgidrp.data import Image, Dataset
 from corgidrp.mocks import (create_default_headers, make_fluxmap_image, nonlin_coefs)
-from corgidrp.calibrate_kgain import (calibrate_kgain, CalKgainException, kgain_params)
+from corgidrp.calibrate_kgain import (calibrate_kgain, CalKgainException, kgain_params_default)
 
+
+np.random.seed(8585)
 ######################## function definitions ###############################
 
 def count_contiguous_repeats(arr):
@@ -135,11 +137,18 @@ def test_expected_results_sub():
     """Outputs are as expected, for imported frames."""
     kgain = calibrate_kgain(dataset_kg, n_cal, n_mean, min_val, max_val, binwidth)
         
-    signal_bins_N = kgain_params['signal_bins_N']
+    signal_bins_N = kgain_params_default['signal_bins_N']
     # kgain - should be close to the assumed value
     assert np.isclose(round(kgain.value,1), kgain_in, atol=0.5)
     assert np.all(np.equal(kgain.ptc.shape, (signal_bins_N,2)))
 
+    # test bad input for kgain_params
+    kgain_params_bad = kgain_params_default.copy()
+    kgain_params_bad['colroi2'] = 'foo'
+    with pytest.raises(TypeError):
+        calibrate_kgain(dataset_kg, n_cal, n_mean, min_val, max_val, binwidth,
+                        kgain_params=kgain_params_bad)
+
 def test_psi():
     """These three below must be positive scalar integers."""
     check_list = test_check.psilist
diff --git a/tests/test_mean_combine.py b/tests/test_mean_combine.py
index da1932c1..7b8b0cbd 100644
--- a/tests/test_mean_combine.py
+++ b/tests/test_mean_combine.py
@@ -34,7 +34,7 @@ def test_mean_im():
 
     check_combined_im = np.mean(check_ims, axis=0)
     check_combined_im_err = np.sqrt(np.sum(np.array(check_ims)**2,
-                                            axis=0))/len(check_ims)
+                                            axis=0)/len(check_ims))
     # For the pixel that is masked throughout
     check_combined_im[all_masks_i] = 0
     check_combined_im_err[all_masks_i] = 0
@@ -42,7 +42,7 @@ def test_mean_im():
     # For pixel that is only masked once
     check_combined_im[single_mask_i] = np.mean(unmasked_vals)
     check_combined_im_err[single_mask_i] = \
-        np.sqrt(np.sum(unmasked_vals**2))/unmasked_vals.size
+        np.sqrt(np.sum(unmasked_vals**2)/unmasked_vals.size)
 
     combined_im, _, _, _ = mean_combine(check_ims, check_masks)
     combined_im_err, _, _, _ = mean_combine(check_ims,
diff --git a/tests/test_nonlin_cal.py b/tests/test_nonlin_cal.py
index 9a5122c8..a1461a16 100644
--- a/tests/test_nonlin_cal.py
+++ b/tests/test_nonlin_cal.py
@@ -15,7 +15,7 @@
 import test_check
 from corgidrp import check
 from corgidrp.data import Image, Dataset
-from corgidrp.calibrate_nonlin import (calibrate_nonlin, CalNonlinException)
+from corgidrp.calibrate_nonlin import (calibrate_nonlin, CalNonlinException, nonlin_params_default)
 from corgidrp.mocks import (create_default_headers, make_fluxmap_image, nonlin_coefs)
 
 ############################# prepare simulated frames #######################
@@ -248,8 +248,19 @@ def test_psi():
     for mmerr in check_list:
         with pytest.raises(TypeError):
             calibrate_nonlin(dataset_nl, n_cal, n_mean, mmerr, min_write, max_write)
+
+def test_nonlin_params():
+    '''Check that the check function for nonlin_params works as expected.'''
+    # test bad input for nonlin_params
+    nonlin_params_bad = nonlin_params_default.copy()
+    nonlin_params_bad['colroi2'] = 'foo'
+    with pytest.raises(TypeError):
+        calibrate_nonlin(dataset_nl, n_cal, n_mean, norm_val, min_write, max_write,
+                        nonlin_params=nonlin_params_bad)
  
 if __name__ == '__main__':
+    print('Running test_nonlin_params')
+    test_nonlin_params()
     print('Running test_expected_results_nom_sub')
     test_expected_results_nom_sub()
     print('Running test_expected_results_time_sub')
diff --git a/tests/test_photon_counting.py b/tests/test_photon_counting.py
new file mode 100644
index 00000000..ac63e819
--- /dev/null
+++ b/tests/test_photon_counting.py
@@ -0,0 +1,164 @@
+import pytest
+import os
+import astropy.time as time
+import corgidrp.mocks as mocks
+import corgidrp
+from astropy.io import fits
+import numpy as np
+import corgidrp.walker as walker
+import corgidrp.caldb as caldb
+import corgidrp.data as data
+import corgidrp.detector as detector
+from corgidrp.photon_counting import get_pc_mean, photon_count, PhotonCountException
+
+def test_negative():
+    """Values at or below the 
+    threshold are photon-counted as 0, including negative values."""
+    fr = np.array([-3,2,0])
+    pc_fr = photon_count(fr, thresh=2)
+    assert (pc_fr == np.array([0,0,0])).all()
+    
+
+def test_pc():
+    '''
+    Tests that a pixel that is masked heavily (reducing the number of usable frames for that pixel) has a bigger err than the average pixel.
+
+    Tests that if masking is all the way through for a certain pixel, that pixel will 
+    be flagged in the DQ. 
+    
+    Make sure there is failure when the number of iterations niter<1.
+
+    Tests that the use of a mask to mask out bright pixels when photon-counting, to ensure failure for bright pixels outside region of interest when no mask used
+    and success when mask used to mask those bright pixels.
+
+    Tests safemode, which makes the function issue warnings instead of crashing (useful for HOWFSC's iterative digging of dark holes).
+
+    Checks various inputs: threshold<0 gives exception, thresh>=em_gain gives exception, providing dataset of VISTYPE='DARK' while also inputting a photon-counted master dark raises exception, exception raised when 'CMDGAIN' not the same for all frames.
+
+    Checks various metadata changes: the expected filename and history edit for the output is done (for case of dark subtraction and the case of no dark subtraction), the master dark is indicated as such via a header keyword 'PC_STAT'.
+
+    Tests that output average value of the photon-counted frames agrees with what was simulated (for the case of dark subtraction and no dark subtraction).
+
+    Tests that if 'ISPC' is in image header and if it is False, get_pc_mean() raises an exception.
+
+    Tests that an exception occurs if the user input "inputmode" is inconsistent with the input dataset type.
+
+    Tests that an exception occurs if the photon-counted master dark does not have the header 'PCTHRESH'.
+
+    Tests that an exception occurs if the photon-counted master dark has a 'PCTHRESH' value than the one to be used on the illuminated dataset.
+    '''
+    # exposure time too long to get reasonable photon-counted result (after photometric correction)
+    np.random.seed(555)
+    dataset_err, dark_dataset_err, ill_mean, dark_mean = mocks.create_photon_countable_frames(Nbrights=160, Ndarks=160, cosmic_rate=0, full_frame=False, smear=False) 
+    # instead of running through walker, just do the pre-processing steps simply
+    # using EM gain=5000 and kgain=7 and bias=20000 and read noise = 100 and QE=0.9 (quantum efficiency), from mocks.create_photon_countable_frames()
+    for f in dataset_err.frames:
+        f.data = f.data.astype(float)*7 - 20000.
+    for f in dark_dataset_err.frames:
+        f.data = f.data.astype(float)*7 - 20000.
+    dataset_err.all_data = dataset_err.all_data.astype(float)*7 - 20000.
+    dark_dataset_err.all_data = dark_dataset_err.all_data.astype(float)*7 - 20000.
+    # masked for half of the bright and half of the dark frames
+    dataset_err.all_dq[80:,3,3] = 1
+    dark_dataset_err.all_dq[80:,3,3] = 1
+    dataset_err.all_dq[:, 2, 2] = 1 #masked all the way through for (2,2)
+    for f in dataset_err.frames: 
+        f.ext_hdr['RN'] = 100
+        f.ext_hdr['KGAIN'] = 7
+    for f in dark_dataset_err.frames: 
+        f.ext_hdr['RN'] = 100
+        f.ext_hdr['KGAIN'] = 7
+    # process the frames to make PC dark
+    pc_dark = get_pc_mean(dark_dataset_err, inputmode='darks')
+    assert pc_dark.ext_hdr['PC_STAT'] == 'photon-counted master dark'
+    # now process illuminated frames and subtract the PC dark
+    pc_dataset_err = get_pc_mean(dataset_err, pc_master_dark=pc_dark)
+    assert pc_dataset_err.frames[0].filename[-7:] == 'pc.fits'
+    assert pc_dataset_err.frames[0].filepath[-7:] == 'pc.fits'
+    history = ''
+    for line in pc_dataset_err.frames[0].ext_hdr["HISTORY"]:
+        history += line
+    assert "Dark-subtracted: yes" in history
+    assert np.isclose(pc_dataset_err.all_data.mean(), ill_mean - dark_mean, rtol=0.01) 
+    # the DQ for that pixel should be 1
+    assert pc_dataset_err[0].dq[2,2] == 1
+    # err for (3,3) is above the 95th percentile of error: 
+    assert pc_dataset_err[0].err[0][3,3]>np.nanpercentile(pc_dataset_err[0].err,95)
+
+    # also when niter<1, exception
+    with pytest.raises(PhotonCountException):
+        get_pc_mean(dataset_err, niter=0)
+    # when thresh<0, exception
+    with pytest.raises(PhotonCountException):
+        get_pc_mean(dataset_err, T_factor=-1, niter=2)
+    # when thresh>=em_gain, exception
+    with pytest.raises(Exception):
+        get_pc_mean(dataset_err, T_factor=50, niter=2)
+    # can't provide master dark while "VISTYPE" of input dataset is 'DARK'
+    with pytest.raises(PhotonCountException):
+        get_pc_mean(dark_dataset_err, pc_master_dark=pc_dark, inputmode='darks')
+    # must have same 'CMDGAIN' and other header values throughout input dataset
+    dark_dataset_err.frames[0].ext_hdr['CMDGAIN'] = 4999
+    with pytest.raises(PhotonCountException):
+        get_pc_mean(dark_dataset_err, inputmode='dark')
+    # change back:
+    dark_dataset_err.frames[0].ext_hdr['CMDGAIN'] = 5000
+
+    # test to make sure PC dark's threshold matches the one used for illuminated frames 
+    with pytest.raises(PhotonCountException):
+        pc_dark.ext_hdr['PCTHRESH'] = 499
+        get_pc_mean(dataset_err, pc_master_dark=pc_dark, inputmode='illuminated')
+    # PC dark should have header 'PCTHRESH'
+    with pytest.raises(PhotonCountException):
+        pc_dark.ext_hdr.pop("PCTHRESH")
+        get_pc_mean(dataset_err, pc_master_dark=pc_dark, inputmode='illuminated')
+    # set it back:
+    pc_dark.ext_hdr['PCTHRESH'] = 500
+    # to trigger pc_ecount_max error
+    for f in dataset_err.frames[:-2]: #all but 2 of the frames
+        f.data[22:40,23:49] = np.abs(f.data.astype(float)[22:40,23:49]*3000)
+    dataset_err.all_data[:-2,22:40,23:49] = np.abs(dataset_err.all_data.astype(float)[:-2,22:40,23:49]*3000)
+    with pytest.raises(Exception):
+        get_pc_mean(dataset_err, pc_ecount_max=1)
+    # now use mask to exclude 22,23 from region of interest:
+    mask = np.zeros_like(dataset_err.all_data[0])
+    mask[22:40,23:49] = 1 #exclude
+    thisfile_dir = os.path.dirname(__file__)
+    mask_filepath = os.path.join(thisfile_dir, 'test_data', 'pc_mask.fits')
+    hdr = fits.Header()
+    prim = fits.PrimaryHDU(header=hdr)
+    img = fits.ImageHDU(mask)
+    hdul = fits.HDUList([prim,img])
+    hdul.writeto(mask_filepath, overwrite=True)
+    pc_dataset_err = get_pc_mean(dataset_err, pc_ecount_max=1, mask_filepath=mask_filepath)
+    copy_pc = pc_dataset_err.all_data.copy()
+    # mask out all the irrelevant pixels:
+    copy_pc[0,22:40,23:49] = np.nan
+    # didn't dark-subtract this time:
+    assert pc_dataset_err.frames[0].filename.endswith('pc_no_ds.fits')
+    assert pc_dataset_err.frames[0].filepath.endswith('pc_no_ds.fits')
+    history = ''
+    for line in pc_dataset_err.frames[0].ext_hdr["HISTORY"]:
+        history += line
+    assert "Dark-subtracted: no" in history
+    assert np.isclose(np.nanmean(copy_pc), ill_mean, rtol=0.01) 
+    # test safemode=False: should get warning instead of exception
+    with pytest.warns(UserWarning):
+        get_pc_mean(dataset_err, safemode=False)
+    # test ISPC header value
+    with pytest.raises(PhotonCountException):
+        dataset_err.frames[0].ext_hdr['ISPC'] = False
+        get_pc_mean(dataset_err, pc_master_dark=pc_dark)
+    # set to True now
+    dataset_err.frames[0].ext_hdr['ISPC'] = True
+    # test inputmode's compatibility with dataset type
+    with pytest.raises(PhotonCountException):
+        get_pc_mean(dataset_err, pc_master_dark=pc_dark, inputmode='darks')
+    with pytest.raises(PhotonCountException):
+        get_pc_mean(dark_dataset_err, inputmode='illuminated')
+
+if __name__ == '__main__':
+    test_pc()
+    test_negative()
+    
+    
\ No newline at end of file
diff --git a/tests/test_walker.py b/tests/test_walker.py
index 53ba5045..daa19d06 100644
--- a/tests/test_walker.py
+++ b/tests/test_walker.py
@@ -14,6 +14,7 @@
 import corgidrp.walker as walker
 import corgidrp.detector as detector
 
+np.random.seed(456)
 
 def test_autoreducing():
     """