Merge pull request #1149 from cta-observatory/fix_flatfield_calculator

Fix arrays not being masked correctly in flatfield calculator

lstchain/calib/camera/calibration_calculator.py

@@ -172,7 +172,7 @@ def calculate_calibration_coefficients(self, event):
 
         # find unusable pixel from pedestal and flat-field data
         unusable_pixels = np.logical_or(status_data.pedestal_failing_pixels,
-                                                   status_data.flatfield_failing_pixels)
+                                        status_data.flatfield_failing_pixels)
 
         signal = ff_data.charge_median - ped_data.charge_median
 
@@ -237,6 +237,10 @@ def calculate_calibration_coefficients(self, event):
         fill_array = np.ones((constants.N_GAINS, constants.N_PIXELS)) * median_pedestal_per_sample
         calib_data.pedestal_per_sample = np.ma.filled(pedestal_per_sample_masked, fill_array)
 
+        # set to zero time corrections of unusable pixels
+        time_correction_masked =  np.ma.array(calib_data.time_correction, mask=calib_data.unusable_pixels)
+        calib_data.time_correction = time_correction_masked.filled(0)
+
         # in the case FF intensity is not sufficiently high, better to scale low gain calibration from high gain results
         if self.use_scaled_low_gain:
             calib_data.unusable_pixels[constants.LOW_GAIN] = calib_data.unusable_pixels[constants.HIGH_GAIN]

lstchain/calib/camera/flatfield.py

@@ -300,11 +300,11 @@ def calculate_time_results(
             'sample_time': (time_start + (trigger_time - time_start) / 2).unix * u.s,
             'sample_time_min': time_start.unix*u.s,
             'sample_time_max': trigger_time.unix*u.s,
-            'time_mean': np.ma.getdata(pixel_mean)*u.ns,
-            'time_median': np.ma.getdata(pixel_median)*u.ns,
-            'time_std': np.ma.getdata(pixel_std)*u.ns,
-            'relative_time_median': np.ma.getdata(relative_median)*u.ns,
-            'time_median_outliers': np.ma.getdata(time_median_outliers),
+            'time_mean': pixel_mean.filled(np.nan)*u.ns,
+            'time_median': pixel_median.filled(np.nan)*u.ns,
+            'time_std': pixel_std.filled(np.nan)*u.ns,
+            'relative_time_median': relative_median.filled(np.nan)*u.ns,
+            'time_median_outliers': time_median_outliers.filled(True),
 
         }
 
@@ -330,13 +330,19 @@ def calculate_relative_gain_results(
             axis=0,
         )
 
+        # mask pixels without defined statistical values (mainly due to hardware problems)
+        pixel_mean = np.ma.array(pixel_mean, mask=np.isnan(pixel_mean))
+        pixel_median = np.ma.array(pixel_median, mask=np.isnan(pixel_median))
+        pixel_std = np.ma.array(pixel_std, mask=np.isnan(pixel_std))
+
         unused_values = np.abs(masked_trace_integral - pixel_mean) > (max_sigma * pixel_std)
+
         # only warn for values discard in the sigma clipping, not those from before
         outliers = unused_values & (~masked_trace_integral.mask)
         check_outlier_mask(outliers, self.log, "flatfield")
 
-        # ignore outliers identified by sigma clipping also for following operations
-        masked_trace_integral.mask = unused_values
+        # add outliers identified by sigma clipping for following operations
+        masked_trace_integral.mask |= unused_values
 
         # median of the median over the camera
         median_of_pixel_median = np.ma.median(pixel_median, axis=1)
@@ -356,26 +362,20 @@ def calculate_relative_gain_results(
         charge_median_outliers = (
             np.logical_or(charge_deviation < self.charge_median_cut_outliers[0] * median_of_pixel_median[:,np.newaxis],
                           charge_deviation > self.charge_median_cut_outliers[1] * median_of_pixel_median[:,np.newaxis]))
-
-
+
         # outliers from standard deviation
         deviation = pixel_std - median_of_pixel_std[:, np.newaxis]
         charge_std_outliers = (
             np.logical_or(deviation < self.charge_std_cut_outliers[0] * std_of_pixel_std[:, np.newaxis],
                           deviation > self.charge_std_cut_outliers[1] * std_of_pixel_std[:, np.newaxis]))
 
-        # mask pixels with NaN mean, due to missing statistics
-        pixels_without_stat = np.where(np.isnan(pixel_mean)==True)
-        charge_median_outliers[pixels_without_stat] = True
-        charge_std_outliers[pixels_without_stat] = True
-
         return {
-            'relative_gain_median': np.ma.getdata(np.ma.median(relative_gain_event, axis=0)),
-            'relative_gain_mean': np.ma.getdata(np.ma.mean(relative_gain_event, axis=0)),
-            'relative_gain_std': np.ma.getdata(np.ma.std(relative_gain_event, axis=0)),
-            'charge_median': np.ma.getdata(pixel_median),
-            'charge_mean': np.ma.getdata(pixel_mean),
-            'charge_std': np.ma.getdata(pixel_std),
-            'charge_std_outliers': np.ma.getdata(charge_std_outliers),
-            'charge_median_outliers': np.ma.getdata(charge_median_outliers),
+            'relative_gain_median': np.ma.median(relative_gain_event, axis=0).filled(np.nan),
+            'relative_gain_mean': np.ma.mean(relative_gain_event, axis=0).filled(np.nan),
+            'relative_gain_std': np.ma.std(relative_gain_event, axis=0).filled(np.nan),
+            'charge_median': pixel_median.filled(np.nan),
+            'charge_mean': pixel_mean.filled(np.nan),
+            'charge_std': pixel_std.filled(np.nan),
+            'charge_std_outliers': charge_std_outliers.filled(True),
+            'charge_median_outliers': charge_median_outliers.filled(True)
         }

lstchain/calib/camera/pedestals.py

@@ -268,13 +268,19 @@ def calculate_pedestal_results(self, trace_integral, masked_pixels_of_sample):
             axis=0,
         )
 
+        # mask pixels without defined statistical values (mainly due to hardware problems)
+        pixel_mean = np.ma.array(pixel_mean, mask=np.isnan(pixel_mean))
+        pixel_median = np.ma.array(pixel_median, mask=np.isnan(pixel_median))
+        pixel_std = np.ma.array(pixel_std, mask=np.isnan(pixel_std))
+
         unused_values = np.abs(masked_trace_integral - pixel_mean) > (max_sigma * pixel_std)
+
         # only warn for values discard in the sigma clipping, not those from before
         outliers = unused_values & (~masked_trace_integral.mask)
         check_outlier_mask(outliers, self.log, "pedestal")
 
-        # ignore outliers identified by sigma clipping also for following operations
-        masked_trace_integral.mask = unused_values
+        # add outliers identified by sigma clipping for following operations
+        masked_trace_integral.mask |= unused_values
 
         # median over the camera
         median_of_pixel_median = np.ma.median(pixel_median, axis=1)
@@ -301,18 +307,13 @@ def calculate_pedestal_results(self, trace_integral, masked_pixels_of_sample):
             deviation < self.charge_median_cut_outliers[0] * std_of_pixel_median[:,np.newaxis],
             deviation > self.charge_median_cut_outliers[1] * std_of_pixel_median[:,np.newaxis],
         )
-
-        # mask pixels with NaN mean, due to missing statistics
-        pixels_without_stat = np.where(np.isnan(pixel_mean)==True)
-        charge_median_outliers[pixels_without_stat] = True
-        charge_std_outliers[pixels_without_stat] = True
-
+
         return {
-            'charge_median': np.ma.getdata(pixel_median),
-            'charge_mean': np.ma.getdata(pixel_mean),
-            'charge_std': np.ma.getdata(pixel_std),
-            'charge_std_outliers': np.ma.getdata(charge_std_outliers),
-            'charge_median_outliers': np.ma.getdata(charge_median_outliers)
+            'charge_median': pixel_median.filled(np.nan),
+            'charge_mean': pixel_mean.filled(np.nan),
+            'charge_std': pixel_std.filled(np.nan),
+            'charge_std_outliers': charge_std_outliers.filled(True),
+            'charge_median_outliers': charge_median_outliers.filled(True)
         }
 
 

lstchain/visualization/plot_calib.py

@@ -1,4 +1,5 @@
 import numpy as np
+from astropy import units as u
 from ctapipe.coordinates import EngineeringCameraFrame
 from ctapipe.visualization import CameraDisplay
 from ctapipe_io_lst import load_camera_geometry
@@ -139,7 +140,7 @@ def plot_calibration_results(ped_data, ff_data, calib_data, run=0, plot_file=Non
                 plt.tight_layout()
                 disp = CameraDisplay(camera)
                 disp.highlight_pixels(mask[chan], linewidth=2)
-                disp.image = image[chan]
+                disp.image = image[chan].to_value(u.ns) 
                 disp.cmap = plt.cm.coolwarm
                 # disp.axes.text(lposx, 0, f'{channel[chan]} time', rotation=90)
                 plt.title(f"{channel[chan]} time")