[AGIPD][CORRECT] New PC data sanitization

cal_tools/cal_tools/agipdlib.py

@@ -556,18 +556,20 @@ class AgipdCorrections:
             data *= rel_cor
             del rel_cor
 
-        # Set negative values for medium gain to 0
-        if self.corr_bools.get('blc_set_min'):
-            data[(data < 0) & (gain == 1)] = 0
-
         # Adjust medium gain baseline to match highest high gain value
         if self.corr_bools.get("adjust_mg_baseline"):
             mgbc = self.md_additional_offset[module_idx][cellid, ...]
             data[gain == 1] += mgbc[gain == 1]
             del mgbc
 
-        # Do xray correction if requested 
-        # The slopes we have in our constants are already relative
+        # Set negative values for medium gain to 0
+        # TODO: Probably it would be better to add it to badpixel maps,
+        # not just set to 0
+        if self.corr_bools.get('blc_set_min'):
+            data[(data < 0) & (gain == 1)] = 0
+            
+        # Do xray correction if requested
+        # The slopes we have in our constants are already relative 
         # slopeFF = slopeFFpix/avarege(slopeFFpix)
         # To apply them we have to / not *
         if self.corr_bools.get("xray_corr"):
@@ -857,11 +859,10 @@ class AgipdCorrections:
           Please note: Current slopesFF avaialble in calibibration
           constants are created per pixel only, not per memory cell:
 
+
              rel_high_gain = 1 if only PC data is available
              rel_high_gain = rel_slopesFF if FF data is also available
 
-          
-
         * Relative gain for Medium gain stage: we derive the factor
           between high and medium gain using slope information from
           fits to the linear part of high and medium gain:
@@ -948,38 +949,59 @@ class AgipdCorrections:
             pc_med_m = slopesPC[..., :self.max_cells, 3]
             pc_med_l = slopesPC[..., :self.max_cells, 4]
 
-            # calculate ratio high to medium gain 
-            pc_high_ave = np.nanmean(pc_high_m, axis=(0,1))
-            pc_med_ave = np.nanmean(pc_med_m, axis=(0,1))
-            # ration between HG and MG per pixel per mem cell 
-            # used for rel gain calculation
+            # calculate median for slopes  
+            pc_high_med = np.nanmedian(pc_high_m, axis=(0,1))
+            pc_med_med = np.nanmedian(pc_med_m, axis=(0,1))
+            # calculate median for intercepts:
+            pc_high_l_med = np.nanmedian(pc_high_l, axis=(0,1))
+            pc_med_l_med = np.nanmedian(pc_med_l, axis=(0,1))
+            
+            # sanitize PC data 
+            # (it should be done already on the level of constants)
+            # In the following loop, 
+            # replace `nan`s across memory cells with
+            # the median value calculated previously.
+            # Then, values outside of the valid range (0.8 and 1.2)
+            # are fixed to the median value.
+            # This is applied for high and medium gain stages
+            for i in range(self.max_cells):
+                pc_high_m[np.isnan(pc_high_m[..., i])] = pc_high_med[i]
+                pc_med_m[np.isnan(pc_med_m[..., i])] = pc_med_med[i]
+
+                pc_high_l[np.isnan(pc_high_l[..., i])] = pc_high_l_med[i]
+                pc_med_l[np.isnan(pc_med_l[..., i])] = pc_med_l_med[i]
+
+                pc_high_m[(pc_high_m[..., i] < 0.8 * pc_high_med[i]) |
+                          (pc_high_m[..., i] > 1.2 * pc_high_med[i])] = pc_high_med[i]  # noqa
+                pc_med_m[(pc_med_m[..., i] < 0.8 * pc_med_med[i]) |
+                         (pc_med_m[..., i] > 1.2 * pc_med_med[i])] = pc_med_med[i]  # noqa
+
+                pc_high_l[(pc_high_l[..., i] < 0.8 * pc_high_l_med[i]) |
+                          (pc_high_l[..., i] > 1.2 * pc_high_l_med[i])] = pc_high_l_med[i]  # noqa
+                pc_med_l[(pc_med_l[..., i] < 0.8 * pc_med_l_med[i]) |
+                         (pc_med_l[..., i] > 1.2 * pc_med_l_med[i])] = pc_med_l_med[i]  # noqa
+
+            # ration between HG and MG per pixel per mem cell used
+            # for rel gain calculation
             frac_high_med_pix = pc_high_m / pc_med_m
             # avarage ration between HG and MG as a function of 
             # mem cell (needed for bls_stripes)
             # TODO: Per pixel would be more optimal correction
-            frac_high_med = pc_high_ave / pc_med_ave
+            frac_high_med = pc_high_med / pc_med_med
             # calculate additional medium-gain offset
             md_additional_offset = pc_high_l - pc_med_l * pc_high_m / pc_med_m
 
-            # Calculate relative gain. If FF constants are available,
-            # use them for high gain
-            # if not rel_gain is calculated using PC data only
-            # if self.corr_bools.get("xray_corr"):
-            #     rel_gain[..., :self.max_cells, 0] /= xray_corr 
-
-            # PC data should be 'calibrated with X-ray data, 
-            # if it is not done, it is better to use 1 instead of bias 
-            # the results with PC arteffacts. 
+            # PC data should be 'calibrated with X-ray data,
+            # if it is not done, it is better to use 1 instead of bias
+            # the results with PC arteffacts.
             # rel_gain[..., 0] = 1./(pc_high_m / pc_high_ave)
-
-            # High-gain (rel_gain[..., 0]) stays the same
             rel_gain[..., 1] = rel_gain[..., 0] * frac_high_med_pix
             rel_gain[..., 2] = rel_gain[..., 1] * 4.48
 
             self.md_additional_offset[module_idx][...] = md_additional_offset.transpose()[...]  # noqa
             self.rel_gain[module_idx][...] = rel_gain[...].transpose()
             self.frac_high_med[module_idx][...] = frac_high_med
-            
+
         self.mask[module_idx][...] = bpixels.transpose()[...]
 
         return