Merge branch 'main' into use_effective_focal_length_for_srcdep

.github/release-drafter.yml

@@ -1,5 +1,5 @@
-name-template: '$NEXT_PATCH_VERSION'
-tag-template: '$NEXT_PATCH_VERSION'
+name-template: 'v$NEXT_PATCH_VERSION'
+tag-template: 'v$NEXT_PATCH_VERSION'
 
 template: |
   ## What’s Changed

environment.yml

@@ -30,8 +30,8 @@ dependencies:
   - pandas>=2
   - pymongo
   - seaborn
-  - ctapipe_io_lst=0.21.1
+  - ctapipe_io_lst=0.22
   - pytest
+  - pyirf=0.8
   - pip:
       - ctaplot~=0.6.2
-      - pyirf~=0.8.0

lstchain/calib/camera/calibration_calculator.py

@@ -48,6 +48,11 @@ class CalibrationCalculator(Component):
         help='Excess noise factor squared: 1+ Var(gain)/Mean(Gain)**2'
     ).tag(config=True)
 
+    relative_qe_dispersion = traits.Float(
+        0.07,
+        help='Relative (effective) quantum efficiency dispersion of PMs over the camera'
+    ).tag(config=True)
+
     pedestal_product = traits.create_class_enum_trait(
         PedestalCalculator,
         default_value='PedestalIntegrator'
@@ -210,10 +215,12 @@ 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 uncertanty (neglect the 7% spread due to detection QE)
+        # 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)
+
         npe_outliers = (
-            np.logical_or(npe_deviation < self.npe_median_cut_outliers[0] * np.sqrt(npe_median)[:,np.newaxis],
-                          npe_deviation > self.npe_median_cut_outliers[1] * np.sqrt(npe_median)[:,np.newaxis]))
+            np.logical_or(npe_deviation < self.npe_median_cut_outliers[0] * tot_std[:,np.newaxis],
+                          npe_deviation > self.npe_median_cut_outliers[1] * tot_std[:,np.newaxis]))
 
         # calibration unusable pixels are an OR of all masks
         calib_data.unusable_pixels = np.logical_or(unusable_pixels, npe_outliers)

lstchain/calib/camera/flatfield.py

@@ -183,6 +183,7 @@ def calculate_relative_gain(self, event):
         self.collect_sample(charge, pixel_mask, arrival_time)
 
         sample_age = (self.trigger_time - self.time_start).to_value(u.s)
+
         # check if to create a calibration event
         if (self.num_events_seen > 0 and
                 (sample_age > self.sample_duration or

lstchain/calib/camera/pixel_threshold_estimation.py

@@ -24,10 +24,10 @@ def get_bias_and_std(dl1_file):
     """
     with tables.open_file(dl1_file) as f:
         ped = f.root[dl1_params_tel_mon_ped_key]
-        ped_charge_mean = np.array(ped.cols.charge_mean)
-        ped_charge_std = np.array(ped.cols.charge_std)
+        ped_charge_mean = ped.col('charge_mean')
+        ped_charge_std = ped.col('charge_std')
         calib = f.root[dl1_params_tel_mon_cal_key]
-        dc_to_pe = np.array(calib.cols.dc_to_pe[ORIGINAL_CALIBRATION_ID])
+        dc_to_pe = calib.col('dc_to_pe')[ORIGINAL_CALIBRATION_ID]
         ped_charge_mean_pe = ped_charge_mean * dc_to_pe
         ped_charge_std_pe = ped_charge_std * dc_to_pe
 

lstchain/data/onsite_camera_calibration_param.json

@@ -7,6 +7,9 @@
     "log_level": "INFO"
   },
   "LSTEventSource": {
+    "MultiFiles": {
+      "all_subruns": true
+     },
     "allowed_tels": [
       1
     ],

lstchain/reco/dl1_to_dl2.py

@@ -14,7 +14,7 @@
 import joblib
 import numpy as np
 import pandas as pd
-from astropy.coordinates import SkyCoord
+from astropy.coordinates import SkyCoord, Angle
 from astropy.time import Time
 from pathlib import Path
 from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor
@@ -417,6 +417,16 @@ def build_models(filegammas, fileprotons,
                                                   + config['disp_classification_features'],
                                     )
 
+    # Normalize all azimuth angles to the range [0, 360) degrees
+    df_gamma.az_tel = Angle(df_gamma.az_tel, u.rad).wrap_at(360 * u.deg).rad
+    df_proton.az_tel = Angle(df_proton.az_tel, u.rad).wrap_at(360 * u.deg).rad
+
+    # Dealing with `sin_az_tel` missing data because of the former version of lstchain
+    if 'sin_az_tel' not in df_gamma.columns:
+        df_gamma['sin_az_tel'] = np.sin(df_gamma.az_tel)
+    if 'sin_az_tel' not in df_proton.columns:
+        df_proton['sin_az_tel'] = np.sin(df_proton.az_tel)
+
     # Training MC gammas in reduced viewcone
     src_r_min = config['train_gamma_src_r_deg'][0]
     src_r_max = config['train_gamma_src_r_deg'][1]

lstchain/scripts/lstchain_create_run_summary.py

@@ -54,6 +54,13 @@
     default="/fefs/aswg/data/real/monitoring/RunSummary",
 )
 
+parser.add_argument(
+    "--overwrite",
+    action="store_true",
+    help="Overwrite existing Run Summary file",
+    default=False,
+)
+
 dtypes = {
     "ucts_timestamp": np.int64,
     "run_start": np.int64,
@@ -306,7 +313,10 @@ def main():
     run_summary.add_column(n_subruns, name="n_subruns", index=1)
     run_summary.add_column(run_types, name="run_type", index=2)
     run_summary.write(
-        args.output_dir / f"RunSummary_{args.date}.ecsv", format="ascii.ecsv", delimiter=",",
+        args.output_dir / f"RunSummary_{args.date}.ecsv",
+        format="ascii.ecsv",
+        delimiter=",",
+        overwrite=args.overwrite,
     )
 
 

lstchain/scripts/lstchain_dl1_to_dl2.py

@@ -12,6 +12,8 @@
 
 import numpy as np
 import pandas as pd
+import astropy.units as u
+from astropy.coordinates import Angle
 from ctapipe.instrument import SubarrayDescription
 from ctapipe_io_lst import OPTICS
 from tables import open_file
@@ -104,6 +106,14 @@ def apply_to_file(filename, models_dict, output_dir, config):
         print("subarray table is not readable because of the version inompatibility.")
         print("Use the effective focal lentgh for the standard LST optics")
         effective_focal_length = OPTICS.effective_focal_length
+
+    # Normalize all azimuth angles to the range [0, 360) degrees 
+    data.az_tel = Angle(data.az_tel, u.rad).wrap_at(360 * u.deg).rad
+
+    # Dealing with `sin_az_tel` missing data because of the former version of lstchain
+    if 'sin_az_tel' not in data.columns:
+        data['sin_az_tel'] = np.sin(data.az_tel)
+
 
     # Apply the models to the data
 

lstchain/scripts/onsite/onsite_create_calibration_file.py

@@ -236,7 +236,7 @@ def main():
         f"--PedestalCalculator.sample_size={stat_events}",
         f"--config={config_file}",
         f"--log-file={log_file}",
-        "--log-file-level=DEBUG",
+        "--log-file-level=INFO",
         *remaining_args,
     ]
 

lstchain/tools/lstchain_create_calibration_file.py

@@ -103,10 +103,11 @@ def __init__(self, **kwargs):
         self.processor = None
         self.writer = None
         self.simulation = False
+        self.n_calib = 0
 
     def setup(self):
 
-        self.log.debug("Opening file")
+        self.log.info("Opening file")
         self.eventsource = EventSource.from_config(parent=self)
 
         self.processor = CalibrationCalculator.from_name(
@@ -130,7 +131,7 @@ def setup(self):
 
         group_name = 'tel_' + str(tel_id)
 
-        self.log.debug(f"Open output file {self.output_file}")
+        self.log.info(f"Open output file {self.output_file}")
 
         self.writer = HDF5TableWriter(
             filename=self.output_file, group_name=group_name, overwrite=True
@@ -202,16 +203,16 @@ def start(self):
 
             # write flatfield results when enough statistics (also for pedestals)
             if (new_ff and new_ped):
-                self.log.debug(f"Write calibration at event n. {count+1}, event id {event.index.event_id} ")
+                self.log.info(f"Write calibration at event n. {count+1}, event id {event.index.event_id} ")
 
-                self.log.debug(f"Ready flatfield data at event n. {count_ff} "
-                                f"stat = {ff_data.n_events} events")
+                self.log.info(f"Ready flatfield data at event n. {count_ff} "
+                              f"stat = {ff_data.n_events} events")
 
                 # write on file
                 self.writer.write('flatfield', ff_data)
 
-                self.log.debug(f"Ready pedestal data at event n. {count_ped}"
-                                f"stat = {ped_data.n_events} events")
+                self.log.info(f"Ready pedestal data at event n. {count_ped} "
+                              f"stat = {ped_data.n_events} events")
 
                 # write only pedestal data used for calibration
                 self.writer.write('pedestal', ped_data)
@@ -223,15 +224,25 @@ def start(self):
                 self.processor.calculate_calibration_coefficients(event)
 
                 # write calib and pixel status
-                self.log.debug("Write pixel_status data")
+                self.log.info("Write pixel_status data")
                 self.writer.write('pixel_status', status_data)
 
-                self.log.debug("Write calibration data")
+                self.log.info("Write calibration data")
                 self.writer.write('calibration', calib_data)
+                self.n_calib += 1
+
                 if self.one_event:
                     break
 
     def finish(self):
+
+        self.log.info(f"Written {self.n_calib} calibration events")
+        if self.n_calib == 0:
+             self.log.critical("!!! No calibration events in the output file !!! : ")
+             self.log.critical(f"flatfield collected statistics = {self.processor.flatfield.num_events_seen} events")
+             self.log.critical(f"pedestal collected statistics = {self.processor.pedestal.num_events_seen} events")
+             self.exit(1)
+
         Provenance().add_output_file(
             self.output_file,
             role='mon.tel.calibration'

setup.py

@@ -42,7 +42,7 @@ def find_scripts(script_dir, prefix):
         'astropy~=5.0',
         'bokeh~=2.0',
         'ctapipe~=0.19.2',
-        'ctapipe_io_lst~=0.21.1',
+        'ctapipe_io_lst~=0.22.0',
         'ctaplot~=0.6.2',
         'eventio>=1.9.1,<2.0.0a0',  # at least 1.1.1, but not 2
         'gammapy~=1.1',