Merge pull request #1147 from cta-observatory/ctaB_online_script

Add Cat-B calibration onsite scripts

lstchain/calib/camera/calibration_calculator.py

@@ -76,7 +76,8 @@ class CalibrationCalculator(Component):
 
     hg_lg_ratio = traits.Float(
         1.,
-        help='HG/LG ratio applied if use_scaled_low_gain is True. In case of calibrated data the ratio should be 1.'
+        help='HG/LG ratio applied if use_scaled_low_gain is True. The ratio is ~1 for calibrated data, ~17.4 for uncalibrated data.' 
+
     ).tag(config=True)
 
     classes = (
@@ -214,8 +215,8 @@ def calculate_calibration_coefficients(self, event):
         # define unusables on number of estimated pe
         npe_deviation =  calib_data.n_pe - npe_median[:,np.newaxis]
 
-        # cut on the base of pe statistical uncertainty (adding a 7% spread due to different detection QE among PMs) 
-        tot_std = np.sqrt(npe_median + (self.relative_qe_dispersion * npe_median)**2)
+        # cut on the base of the pe statistical uncertainty over the camera
+        tot_std = self.expected_npe_std(npe_median, ff_data.n_events)
 
         npe_outliers = (
             np.logical_or(npe_deviation < self.npe_median_cut_outliers[0] * tot_std[:,np.newaxis],
@@ -296,3 +297,41 @@ def output_interleaved_results(self, event):
                 new_ff = True
 
         return new_ped, new_ff
+
+    def expected_npe_std(self, npe_median, n_events):
+
+        """
+        The expected standard deviation of the estimated npe over the camera is given in principle by
+
+        std_pe_mean=std_npe/sqrt((n_events)+ (relative_qe_dispersion*npe)**2)
+
+        where the relative_qe_dispersion is mainly due to different detection QE among PMs.
+
+        However, due to the systematics correction associated to the B term, a linear and quadratic 
+        noise component  must be added, these components depend on the sample statistics per pixel (n_events).
+        The parameters in this function (linear_noise_params and quadratic_noise_params) have been obtained with 
+        a fit of the std of filter scan taken in date 2023/05/10 and considering 
+        n_events = [1000,2500,5000,7500,10000,30000]
+
+        """
+        basic_variance = npe_median/n_events + (self.relative_qe_dispersion * npe_median)**2
+
+        # [par0_HG, par0_LG],[par1_HG, par1_LG]
+        linear_noise_params=np.array(([1.79717813 ,1.72458305e+00],[0.0231544, -1.62036639e-03]))
+
+        # [par0_HG, par0_LG],[par1_HG, par1_LG]
+        quadratic_noise_params=np.array(([4.99670969e-04, 0.00142218],[ 2.49034290e-05,0.0001207]))
+
+        # function to estimate the added noise components as function of the sample statistcs
+        def noise_term(n_events, par):
+            return par[0]/(np.sqrt(n_events))+par[1]
+
+        linear_term = noise_term(n_events, linear_noise_params)
+        quadratic_term = noise_term(n_events, quadratic_noise_params)
+
+        added_variance = (linear_term * npe_median)**2 + (quadratic_term * npe_median**2)**2 
+
+        std = np.sqrt(basic_variance + added_variance)
+
+        return std
+

lstchain/calib/camera/pedestals.py

@@ -106,6 +106,7 @@ def __init__(self, subarray, **kwargs):
         self.extractor = ImageExtractor.from_name(
             self.charge_product, parent=self, subarray=subarray
         )
+
 
     def _extract_charge(self, event):
         """
@@ -117,7 +118,7 @@ def _extract_charge(self, event):
         event : general event container
 
         """
-
+        
         # copy the waveform be cause we do not want to change it for the moment
         waveforms = np.copy(event.r1.tel[self.tel_id].waveform)
 
@@ -155,17 +156,18 @@ def calculate_pedestals(self, event):
         Returns: True if the mon.tel[tel_id].pedestal is updated, False otherwise
 
         """
+
         # initialize the np array at each cycle
         waveform = event.r1.tel[self.tel_id].waveform
-
+    
         # re-initialize counter
         if self.num_events_seen == self.sample_size:
             self.num_events_seen = 0
 
         pixel_mask = event.mon.tel[self.tel_id].pixel_status.hardware_failing_pixels
 
         self.trigger_time = event.trigger.time
-
+        
         if self.num_events_seen == 0:
             self.time_start = self.trigger_time
             self.setup_sample_buffers(waveform, self.sample_size)
@@ -231,7 +233,7 @@ def store_results(self, event):
 
     def setup_sample_buffers(self, waveform, sample_size):
         """Initialize sample buffers"""
-
+       
         n_channels = waveform.shape[0]
         n_pix = waveform.shape[1]
         shape = (sample_size, n_channels, n_pix)
@@ -257,7 +259,7 @@ def calculate_pedestal_results(self, trace_integral, masked_pixels_of_sample):
             trace_integral,
             mask=masked_pixels_of_sample
         )
-
+        
         # mean and std over the sample per pixel
         max_sigma = self.sigma_clipping_max_sigma
         pixel_mean, pixel_median, pixel_std = sigma_clipped_stats(

lstchain/conftest.py

@@ -106,11 +106,13 @@ def observed_dl1_files(temp_dir_observed_files, run_summary_path):
     datacheck_file1 = temp_dir_observed_files / "datacheck_dl1_LST-1.Run02008.0000.h5"
     dvr_file1 = temp_dir_observed_files / "DVR_settings_LST-1.Run02008.h5"
     pixmasks_file1 = temp_dir_observed_files / "Pixel_selection_LST-1.Run02008.0000.h5"
-
+    interleaved_file1 = temp_dir_observed_files / "interleaved/interleaved_LST-1.Run02008.0000.h5"
+
     # Second set of files
     dl1_output_path2 = temp_dir_observed_files / "dl1_LST-1.Run02008.0100.h5"
     muons_file2 = temp_dir_observed_files / "muons_LST-1.Run02008.0100.fits"
     datacheck_file2 = temp_dir_observed_files / "datacheck_dl1_LST-1.Run02008.0100.h5"
+    interleaved_file2 = temp_dir_observed_files / "interleaved/interleaved_LST-1.Run02008.0100.h5"
 
     run_program(
         "lstchain_data_r0_to_dl1",
@@ -198,12 +200,42 @@ def observed_dl1_files(temp_dir_observed_files, run_summary_path):
         'datacheck1': datacheck_file1,
         'dvr_file1': dvr_file1,
         'pixmasks_file1': pixmasks_file1,
+        'interleaved_file1': interleaved_file1,
         'dl1_file2': dl1_output_path2,
         'muons2': muons_file2,
-        'datacheck2': datacheck_file2
+        'datacheck2': datacheck_file2,
+        'interleaved_file2': interleaved_file2
     }
 
 
+
+@pytest.mark.private_data
+@pytest.fixture(scope="session")
+def interleaved_r1_file(temp_dir_observed_files, run_summary_path):
+    test_pedcal_run = test_data / 'real/R0/20200218/LST-1.1.Run02006.0000_first50.fits.fz'
+
+    run_program(
+        "lstchain_data_r0_to_dl1",
+        "-f",
+        test_pedcal_run,
+        "-o",
+        temp_dir_observed_files,
+        "--pedestal-file",
+        test_drs4_pedestal_path,
+        "--calibration-file",
+        test_calib_path,
+        "--time-calibration-file",
+        test_time_calib_path,
+        "--pointing-file",
+        test_drive_report,
+        '--run-summary-path',
+        run_summary_path,
+        "--default-trigger-type=tib"
+    )
+
+    return temp_dir_observed_files / "interleaved/interleaved_LST-1.Run02006.0000.h5"
+
+
 @pytest.fixture(scope="session")
 def simulated_dl2_file(temp_dir_simulated_files, simulated_dl1_file, rf_models):
     """

lstchain/io/io.py

@@ -82,6 +82,10 @@
 dl1_params_tel_mon_ped_key = "/dl1/event/telescope/monitoring/pedestal"
 dl1_params_tel_mon_cal_key = "/dl1/event/telescope/monitoring/calibration"
 dl1_params_tel_mon_flat_key = "/dl1/event/telescope/monitoring/flatfield"
+dl1_mon_tel_CatB_ped_key = "/dl1/monitoring/telescope/catB/pedestal"
+dl1_mon_tel_CatB_cal_key = "/dl1/monitoring/telescope/catB/calibration"
+dl1_mon_tel_CatB_flat_key = "/dl1/monitoring/telescope/catB/flatfield"
+
 dl1_params_lstcam_key = "/dl1/event/telescope/parameters/LST_LSTCam"
 dl1_images_lstcam_key = "/dl1/event/telescope/image/LST_LSTCam"
 dl2_params_lstcam_key = "/dl2/event/telescope/parameters/LST_LSTCam"
@@ -881,6 +885,7 @@ def write_calibration_data(writer, mon_index, mon_event, new_ped=False, new_ff=F
     mon_event.flatfield.prefix = ''
     mon_event.calibration.prefix = ''
     mon_index.prefix = ''
+    monitoring_table='telescope/monitoring'
 
     # update index
     if new_ped:
@@ -893,20 +898,20 @@ def write_calibration_data(writer, mon_index, mon_event, new_ped=False, new_ff=F
     if new_ped:
         # write ped container
         writer.write(
-            table_name="telescope/monitoring/pedestal",
+            table_name=f"{monitoring_table}/pedestal",
             containers=[mon_index, mon_event.pedestal],
         )
 
     if new_ff:
         # write calibration container
         writer.write(
-            table_name="telescope/monitoring/flatfield",
+            table_name=f"{monitoring_table}/flatfield",
             containers=[mon_index, mon_event.flatfield],
         )
 
         # write ff container
         writer.write(
-            table_name="telescope/monitoring/calibration",
+            table_name=f"{monitoring_table}/calibration",
             containers=[mon_index, mon_event.calibration],
         )