[AGIPD][CORRECT] New PC data sanitization

Pulse Capacitor sanitization developed by Jola.

TESTS:

• Tested changes by @jsztuk for better relative gain correction, I understand that there was a good feedback as well from beam-line scientists for such modifications

• Extra Xfel-calibrate CL test

xfel-calibrate AGIPD CORRECT \
--slurm-mem 750 \
--karabo-da -1 \
--receiver-id {}CH0 \
--karabo-id-control MID_IRU_AGIPD1M1 \
--karabo-da-control AGIPD1MCTRL00 \
--h5path-ctrl /CONTROL/{}/MDL/FPGA_COMP_TEST \
--overwrite --sequences-per-node 1 \
--in-folder /gpfs/exfel/exp/MID/202030/p900137/raw \
--out-folder /gpfs/exfel/data/scratch/ahmedk/test/AGIPD_MID-pc \
--karabo-id MID_DET_AGIPD1M-1 \
--run 449 \
--force-hg-if-below \
--force-mg-if-below \
--low-medium-gap \
--zero-nans \
--sequences 0,1,2 \
--rel-gain \
--common-mode \
--blc-stripes

Reviewers: @danilevc @kamile,

cal_tools/cal_tools/agipdlib.py

@@ -58,10 +58,10 @@ def get_acq_rate(fast_paths: Tuple[str, str, int],
         with h5py.File(slow_data_file, "r") as fin:
             if slow_data_path in fin:
                 # The acquisition rate value is stored in a 1D array of type
-                # float. Use the 3rd value, arbitrarily chosen.
-                # It is desired to loose precision here because the usage is
-                # about bucketing the rate for managing meta-data.
-                return round(float(fin[slow_data_path][3]), 1)
+                # float. Use the 3rd value, arbitrarily chosen. It's okay to
+                # loose precision here because the usage is about defining the
+                # rate for meta-data.
+                return round(fin[slow_data_path][3], 1)
 
     # Compute acquisition rate from fast data
     fast_data_file, karabo_id, module = fast_paths
@@ -204,7 +204,7 @@ class AgipdCorrections:
         self.min_pulse = self.pulses_lst[0]
         self.max_pulse = self.pulses_lst[-1]
         self.max_cells = max_cells
-
+        
         # Correction parameters
         self.baseline_corr_noise_threshold = -1000
         self.snow_resolution = SnowResolution.INTERPOLATE
@@ -273,7 +273,7 @@ class AgipdCorrections:
             f = h5py.File(file_name, 'r')
             group = f[data_path]
 
-            _, first_index, last_index, __, valid_indices = \
+            _, first_index, last_index, _, valid_indices = \
                 self.get_valid_image_idx(idx_base, f)
             firange = self.gen_valid_range(first_index, last_index,
                                            self.max_cells, agipd_base, f,
@@ -392,8 +392,7 @@ class AgipdCorrections:
 
                 # Cell common mode
                 cell_cm_sum, cell_cm_count = \
-                    calgs.sum_and_count_in_range_cell(asic_data, dark_min,
-                                                      dark_max)
+                    calgs.sum_and_count_in_range_cell(asic_data, dark_min, dark_max)
                 cell_cm = cell_cm_sum / cell_cm_count
 
                 cell_cm[cell_cm_count < fraction * 32 * 256] = 0
@@ -401,8 +400,7 @@ class AgipdCorrections:
 
                 # Asics common mode
                 asic_cm_sum, asic_cm_count = \
-                    calgs.sum_and_count_in_range_asic(asic_data, dark_min,
-                                                      dark_max)
+                    calgs.sum_and_count_in_range_asic(asic_data, dark_min, dark_max)
                 asic_cm = asic_cm_sum / asic_cm_count
 
                 asic_cm[asic_cm_count < fraction * 64 * 64] = 0
@@ -557,22 +555,23 @@ 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
-        # slopeFF = slopeFFpix/avarege(slopeFFpix)
+        # 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"):
-            data /= self.xray_cor[module_idx]
+            data /= self.xray_cor[module_idx][cellid, ...]
 
         if self.corr_bools.get('melt_snow'):
             melt_snowy_pixels(raw_data, data, gain, rgain,
@@ -848,29 +847,27 @@ class AgipdCorrections:
         exists of the current AGIPD instances.
 
         Relative gain is derived both from pulse capacitor as well as low
-        intensity flat field data, information from flat field data is 
-        needed to 'calibrate' pulse capacitor data, if there is no 
-        available FF data, relative gain for High Gain stage is set to 1:
+        intensity  flat field data, information from FF data is needed to 
+        'calibrate' pulse capacitor data, if there is no available FF data,
+        relative gain for High Gain stage is set to 1:
 
-        * Relative gain for High gain stage - from the FF data we get 
-          the relative slopes of a given pixel and memory cells with 
-          respect to all memory cells and all pixels in the module,
-          Please note: Current slopesFF avaialble in calibibration
-          constants are created per pixel only, not per memory cell:
+        * Relative gain for High gain stage - from the FF data we get the relative slopes 
+          of a given pixel and memory cells with respect to all memory cells and all pixels in the module,
+          Please note: Current slopesFF avaialble in calib. 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:
+        * 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:
 
              rfpc_high_medium = m_h/m_m          
 
-          where m_h and m_m is the medium gain slope of given memory cells 
-          and pixel and m_h is the high gain slope as above
+          where m_h and m_m is the medium gain slope of given memory cells and 
+          pixel and m_h is the high gain slope as above
            rel_gain_medium = rel_high_gain * rfpc_high_medium
 
         With this data the relative gain for the three gain stages evaluates
@@ -949,38 +946,49 @@ 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)
+            
+            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]
+                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]
+    
+                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]
+                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]
+            
+            # 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
@@ -1113,7 +1121,7 @@ class AgipdCorrections:
                                                              dtype='f4')
 
             self.mask[module_idx] = sharedmem.empty(constant_shape, dtype='i4')
-            self.xray_cor[module_idx] = sharedmem.empty(constant_shape[2:],
+            self.xray_cor[module_idx] = sharedmem.empty(constant_shape[0:],
                                                         dtype='f4')
 
     def allocate_images(self, shape, n_cores_files):