Merge tag 'v0.10.19' into merge-v011

Author: rlopezcoto

.github/workflows/ci.yml

@@ -39,8 +39,8 @@ jobs:
 
     strategy:
       matrix:
-        python-version: ["3.10", "3.11"]
-        ctapipe-version: ["v0.19.2"]
+        python-version: ["3.11", "3.12", "3.13"]
+        ctapipe-version: ["v0.25.1"]
 
     steps:
       - uses: actions/checkout@v4

environment.yml

@@ -2,25 +2,26 @@ name: lst-dev
 channels:
   - conda-forge
 dependencies:
-  - python=3.11
+  - python=3.13
   - numpy
-  - astropy=5
+  - astropy>=6.1,<8
   - pip
-  - ctapipe=0.19.2
+  - ctapipe >=0.25.0
+  - ctapipe_io_lst >=0.26.0,<0.28.0a0
   - ctaplot~=0.6.4
-  - gammapy=1.1
+  - gammapy=1.3
   - h5py
   - ipython
   - jupyter
-  - matplotlib=3.7
+  - matplotlib=3.10
   - notebook
   - iminuit>=2
-  - joblib~=1.2.0
+  - joblib~=1.4.0
   - toml
-  - protozfits=2.6
+  - protozfits=2.7
   - pyparsing
-  - scipy~=1.11.4
-  - scikit-learn=1.2
+  - scipy
+  - scikit-learn=1.6
   - sphinx
   - sphinx-automodapi
   - sphinx_rtd_theme
@@ -32,6 +33,5 @@ dependencies:
   - pandas>=2
   - pymongo
   - seaborn
-  - ctapipe_io_lst >=0.25.1,<0.26a0
   - pytest
-  - pyirf~=0.10.0
+  - pyirf~=0.12.0

lstchain/calib/camera/flatfield.py

@@ -135,15 +135,15 @@ def _extract_charge(self, event):
         peak_pos = 0
         if self.extractor:
             broken_pixels = event.mon.tel[self.tel_id].pixel_status.hardware_failing_pixels
-            dl1 = self.extractor(waveforms, self.tel_id, no_gain_selection, broken_pixels=broken_pixels)
+            dl1 = self.extractor(waveforms, self.tel_id, selected_gain_channel=None, broken_pixels=broken_pixels)
             charge = dl1.image
             peak_pos = dl1.peak_time
 
         # shift the time if time shift is already defined
         # (e.g. drs4 waveform time shifts for LST)
         time_shift = event.calibration.tel[self.tel_id].dl1.time_shift
         if time_shift is not None:
-                peak_pos -= time_shift
+            peak_pos -= time_shift
 
         return charge, peak_pos
 

lstchain/calib/camera/pedestals.py

@@ -137,7 +137,10 @@ def _extract_charge(self, event):
         peak_pos = 0
         if self.extractor:
             broken_pixels = event.mon.tel[self.tel_id].pixel_status.hardware_failing_pixels
-            dl1 = self.extractor(waveforms, self.tel_id, no_gain_selection, broken_pixels=broken_pixels)
+         
+            dl1 = self.extractor(waveforms, self.tel_id, 
+                                 selected_gain_channel=None, 
+                                 broken_pixels=broken_pixels)
             charge = dl1.image
             peak_pos = dl1.peak_time
 

lstchain/data/lstchain_standard_config.json

@@ -287,13 +287,14 @@
   "write_interleaved_events":{
     "DataWriter": {
       "overwrite":  true,
-      "write_images": false,
-      "write_parameters": false,
-      "write_waveforms": true,
+      "write_dl1_images": false,
+      "write_dl1_parameters": false,
+      "write_r0_waveforms": false,
+      "write_r1_waveforms": true,
       "transform_waveform": true,
       "waveform_dtype": "uint16",
       "waveform_offset": 400,
-      "waveform_scale": 80  
+      "waveform_scale": 60  
     }
   }
 }

lstchain/datachecks/containers.py

@@ -329,7 +329,8 @@ def fill_pixel_wise_info(self, table, mask, histogram_binnings,
         pointing = SkyCoord(az=self.mean_az_tel,
                             alt=self.mean_alt_tel,
                             frame=horizon_frame)
-        bright_stars = get_bright_stars(pointing=pointing, radius=3*u.deg,
+        bright_stars = get_bright_stars(time=obstime,
+                                        pointing=pointing, radius=3*u.deg,
                                         magnitude_cut=8)
         # Account for average relative spot shift (outwards) due to coma
         # aberration:

lstchain/datachecks/dl1_checker.py

@@ -387,8 +387,8 @@ def plot_datacheck(datacheck_filename, out_path=None, batch=False,
     engineering_geom = geom.transform_to(EngineeringCameraFrame())
 
     # For future bokeh-based display, turned off for now:
-    # page1 = Panel()
-    # page2 = Panel()
+    # page1 = TabPanel()
+    # page2 = TabPanel()
 
     with PdfPages(pdf_filename) as pdf:
         # first deal with the DL1 datacheck file, created from DL1 event data:

lstchain/image/modifier.py

@@ -711,7 +711,7 @@ def initialise_waveforms(self, waveform, dt, tel_id):
 
         """
         log.info(f"Pre-generating nsb waveforms for nsb tuning and telescope id {tel_id}.")
-        n_pixels, n_samples = waveform.shape
+        _, n_pixels, n_samples = waveform.shape
         baseline_correction = -(self.added_nsb[tel_id] * dt).to_value("")
         nsb_waveforms = np.full((self.multiplicity * n_pixels, 2, n_samples), baseline_correction)
         duration = (self.extra_samples + n_samples) * dt
@@ -768,12 +768,12 @@ def _tune_nsb_on_waveform_precomputed(self, waveform, tel_id, is_high_gain, suba
             sampling_rate = readout.sampling_rate
             dt = (1.0 / sampling_rate).to(u.s)
             self.initialise_waveforms(waveform, dt, tel_id)
-        n_pixels, _ = waveform.shape
+        _, n_pixels, _ = waveform.shape
         # The nsb_waveform array is randomised along the first axis,
         # then the n_pixels first elements with correct gain are used for the injection
         self.rng.shuffle(self.nsb_waveforms[tel_id])
-        waveform += ((is_high_gain == 0)[:, np.newaxis] * self.nsb_waveforms[tel_id][:n_pixels, 0] +
-                     (is_high_gain == 1)[:, np.newaxis] * self.nsb_waveforms[tel_id][:n_pixels, 1])
+        waveform[0] += ((is_high_gain == 0)[:, np.newaxis] * self.nsb_waveforms[tel_id][:n_pixels, 0] +
+                        (is_high_gain == 1)[:, np.newaxis] * self.nsb_waveforms[tel_id][:n_pixels, 1])
 
     def _tune_nsb_on_waveform_direct(self, waveform, tel_id, is_high_gain, subarray):
         """
@@ -790,7 +790,7 @@ def _tune_nsb_on_waveform_direct(self, waveform, tel_id, is_high_gain, subarray)
         subarray: `ctapipe.instrument.subarray.SubarrayDescription`
 
         """
-        n_pixels, n_samples = waveform.shape
+        _, n_pixels, n_samples = waveform.shape
         readout = subarray.tel[tel_id].camera.readout
         sampling_rate = readout.sampling_rate
         dt = (1.0 / sampling_rate).to(u.s)

lstchain/image/muon/muon_analysis.py

@@ -248,9 +248,7 @@ def analyze_muon_event(subarray, tel_id, event_id, image, good_ring_config, plot
     pix_ring_2 = image[dist_mask_2]
 
     muonparameters = MuonParametersContainer()
-    muonparameters.containment = ring_containment(
-        muonringparam.radius,
-        muonringparam.center_fov_lon, muonringparam.center_fov_lat, cam_rad)
+    muonparameters.containment = ring_containment(muonringparam, cam_rad)
 
     radial_distribution = radial_light_distribution(
         muonringparam.center_fov_lon,
@@ -291,12 +289,12 @@ def analyze_muon_event(subarray, tel_id, event_id, image, good_ring_config, plot
 
         # We do the calculation of the ring completeness (i.e. fraction of whole circle) using the pixels
         # within the "width" fitted using MuonIntensityFitter
+
         muonparameters.completeness = ring_completeness(
-            x[dist_ringwidth_mask], y[dist_ringwidth_mask],
-            image[dist_ringwidth_mask],
-            muonringparam.radius,
-            muonringparam.center_fov_lon,
-            muonringparam.center_fov_lat,
+            pixel_fov_lon=x[dist_ringwidth_mask], 
+            pixel_fov_lat=y[dist_ringwidth_mask],
+            weights=image[dist_ringwidth_mask],
+            ring=muonringparam,
             threshold=params['ring_completeness_threshold'],
             bins=30)
 

lstchain/image/tests/test_modifier.py

@@ -85,10 +85,7 @@ def test_tune_nsb_on_waveform():
     from lstchain.image.modifier import WaveformNsbTuner
     from lstchain.data.normalised_pulse_template import NormalizedPulseTemplate
 
-    waveform = np.array(
-        [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
-         [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]]
-    )
+    waveform = np.zeros((2, 2, 22))
     added_nsb_rate = {1: 0.1 * u.GHz}
     subarray = LSTEventSource.create_subarray(1)
     amplitude_HG = np.zeros(40)
@@ -108,21 +105,19 @@ def test_tune_nsb_on_waveform():
     charge_spe_cumulative_pdf = interp1d(spe_integral, spe[0], kind='cubic',
                                          bounds_error=False, fill_value=0.,
                                          assume_sorted=True)
-    nsb_tunner = WaveformNsbTuner(added_nsb_rate,
-                                  pulse_templates,
-                                  charge_spe_cumulative_pdf,
-                                  pre_computed_multiplicity=0)
-    nsb_tunner.tune_nsb_on_waveform(waveform, 1, gain, subarray)
+    nsb_tuner = WaveformNsbTuner(added_nsb_rate,
+                                 pulse_templates,
+                                 charge_spe_cumulative_pdf,
+                                 pre_computed_multiplicity=0)
+    nsb_tuner.tune_nsb_on_waveform(waveform, 1, gain, subarray)
     assert np.any(waveform != 0)
     assert np.isclose(np.mean(waveform), 0.0, atol=0.2)
-    waveform = np.array(
-        [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
-         [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]]
-    )
-    nsb_tunner = WaveformNsbTuner(added_nsb_rate,
-                                  pulse_templates,
-                                  charge_spe_cumulative_pdf,
-                                  pre_computed_multiplicity=10)
-    nsb_tunner.tune_nsb_on_waveform(waveform, 1, gain, subarray)
+    
+    waveform = np.zeros((2, 2, 22))
+    nsb_tuner = WaveformNsbTuner(added_nsb_rate,
+                                 pulse_templates,
+                                 charge_spe_cumulative_pdf,
+                                 pre_computed_multiplicity=10)
+    nsb_tuner.tune_nsb_on_waveform(waveform, 1, gain, subarray)
     assert np.any(waveform != 0)
     assert np.isclose(np.mean(waveform), 0.0, atol=0.2)

lstchain/mc/sensitivity.py

@@ -1,7 +1,7 @@
 import astropy.units as u
 import numpy as np
 import pandas as pd
-from astropy.coordinates.angle_utilities import angular_separation
+from astropy.coordinates import angular_separation
 from gammapy.stats import WStatCountsStatistic
 
 from lstchain.io import read_simu_info_merged_hdf5

lstchain/reco/r0_to_dl1.py

@@ -510,7 +510,7 @@ def r0_to_dl1(
                 event.simulation.prefix = 'mc'
 
             dl1_container.reset()
-
+            
             # write sub tables
             if is_simu:
                 write_subarray_tables(writer, event, metadata)
@@ -536,6 +536,13 @@ def r0_to_dl1(
                                            event.mon.tel[tel_id],
                                            new_ped=True, new_ff=True)
 
+                # Default of absolute_factor is None - needs to be set here
+                event.calibration.tel[tel_id].dl1.absolute_factor = np.ones((2, PIXEL_INDEX.size))
+                # PATCH: ctapipe expects relative_factor to be always (ngains, npixels)
+                if event.calibration.tel[tel_id].dl1.relative_factor is not None:
+                    if event.calibration.tel[tel_id].dl1.relative_factor.ndim == 1:
+                        event.calibration.tel[tel_id].dl1.relative_factor = np.array(2*[event.calibration.tel[tel_id].dl1.relative_factor])
+
                 # flat-field or pedestal:
                 if ((event.trigger.event_type == EventType.FLATFIELD or
                      event.trigger.event_type == EventType.SKY_PEDESTAL) and
@@ -568,14 +575,10 @@ def r0_to_dl1(
 
                     r1 = event.r1.tel[tel_id]
                     r1.selected_gain_channel = source.r0_r1_calibrator.gain_selector(event.r0.tel[tel_id].waveform)
-                    r1.waveform = r1.waveform[r1.selected_gain_channel, PIXEL_INDEX]
-
-                    event.calibration.tel[tel_id].dl1.time_shift = \
-                        event.calibration.tel[tel_id].dl1.time_shift[r1.selected_gain_channel, PIXEL_INDEX]
-
-                    event.calibration.tel[tel_id].dl1.relative_factor = \
-                        event.calibration.tel[tel_id].dl1.relative_factor[r1.selected_gain_channel, PIXEL_INDEX]
+                    # select gain but keep waveform 3d
+                    r1.waveform = r1.waveform[np.newaxis, r1.selected_gain_channel, PIXEL_INDEX]
 
+    
             # Option to add nsb in waveforms
             if nsb_tuning:
                 # FIXME? assumes same correction ratio for all telescopes
@@ -702,15 +705,15 @@ def r0_to_dl1(
                     if tag_pix_thr(image):
 
                         # re-calibrate r1 to obtain new dl1, using a more adequate pulse integrator for muon rings
-                        numsamples = event.r1.tel[telescope_id].waveform.shape[1]  # not necessarily the same as in r0!
+                        numsamples = event.r1.tel[telescope_id].waveform.shape[2]  # not necessarily the same as in r0!
                         bad_pixels_hg = calibration_mon.unusable_pixels[0]
                         bad_pixels_lg = calibration_mon.unusable_pixels[1]
 
                         # Now set to 0 all samples in unreliable pixels. Important for global peak
                         # integrator in case of crazy pixels!  TBD: can this be done in a simpler
                         # way?
                         bad_pixels = bad_pixels_hg | bad_pixels_lg
-                        bad_waveform = np.transpose(np.array(numsamples * [bad_pixels]))
+                        bad_waveform = np.array([np.transpose(np.array(numsamples * [bad_pixels]))]) # (1, npixels, nsamples)
 
                         # print('hg bad pixels:',np.where(bad_pixels_hg))
                         # print('lg bad pixels:',np.where(bad_pixels_lg))
@@ -745,8 +748,8 @@ def r0_to_dl1(
                             mask_hg = bright_pixels & (selected_gain == 0)
                             mask_lg = bright_pixels & (selected_gain == 1)
 
-                            bright_pixels_waveforms_hg = event.r1.tel[telescope_id].waveform[mask_hg, :]
-                            bright_pixels_waveforms_lg = event.r1.tel[telescope_id].waveform[mask_lg, :]
+                            bright_pixels_waveforms_hg = event.r1.tel[telescope_id].waveform[0, mask_hg, :]
+                            bright_pixels_waveforms_lg = event.r1.tel[telescope_id].waveform[0, mask_lg, :]
                             stacked_waveforms_hg = np.sum(bright_pixels_waveforms_hg, axis=0)
                             stacked_waveforms_lg = np.sum(bright_pixels_waveforms_lg, axis=0)
 

lstchain/reco/reconstructor.py

@@ -154,14 +154,14 @@ def get_ped_from_true_signal_less(self, source, nsb_tunner):
             if event.trigger.event_type == EventType.SUBARRAY:
                 for tel_id in event.r1.tel.keys():
                     if i == 0:
-                        n_pix[tel_id] = event.r1.tel[tel_id].waveform.shape[0]
+                        n_pix[tel_id] = event.r1.tel[tel_id].waveform.shape[1]
                     mask = event.simulation.tel[tel_id].true_image == 0
                     wave = event.r1.tel[tel_id].waveform
                     if nsb_tunner is not None:
                         selected_gains = event.r1.tel[tel_id].selected_gain_channel
                         mask_high = (selected_gains == 0)
                         nsb_tunner.tune_nsb_on_waveform(wave, tel_id, mask_high, self.subarray)
-                    waveforms[tel_id].append(wave[mask].flatten())
+                    waveforms[tel_id].append(wave[0, mask].flatten())
         for tel_id, tel_waveforms in waveforms.items():
             self.error[tel_id] = np.full(n_pix[tel_id], np.nanstd(np.concatenate(tel_waveforms)))
             self.allowed_pixels = np.ones(n_pix[tel_id], dtype=bool)
@@ -210,14 +210,15 @@ def call_setup(self, event, telescope_id, dl1_container):
         waveform = event.r1.tel[telescope_id].waveform
 
         dl1_calib = event.calibration.tel[telescope_id].dl1
-        time_shift = dl1_calib.time_shift
         # TODO check if this is correct here or if it is applied to r1 waveform earlier
         if dl1_calib.pedestal_offset is not None:
             waveform = waveform - dl1_calib.pedestal_offset[:, np.newaxis]
 
-        n_pixels, n_samples = waveform.shape
+        n_samples = waveform.shape[2]
         times = np.arange(0, n_samples) * dt
         selected_gains = event.r1.tel[telescope_id].selected_gain_channel
+        time_shift = dl1_calib.time_shift
+        
         is_high_gain = (selected_gains == 0)
 
         # We assume that the time gradient is given in unit of 'geometry spatial unit'/ns
@@ -283,13 +284,15 @@ def call_setup(self, event, telescope_id, dl1_container):
         crosstalks = spatial_ones * self.crosstalk.tel[telescope_id]
 
         times = (np.arange(0, n_samples) * dt)[mask_time]
+
+        time_shift = np.choose(selected_gains, time_shift)
         time_shift = time_shift[mask_pixel]
 
         p_x = pix_x[mask_pixel]
         p_y = pix_y[mask_pixel]
         pix_area = geometry.pix_area[mask_pixel].to_value(unit ** 2)
 
-        data = waveform
+        data = waveform[0]
         error = self.error
 
         filter_pixels = np.nonzero(~mask_pixel)

lstchain/reco/tests/test_utils.py

@@ -30,8 +30,8 @@ def test_camera_to_altaz():
 def test_radec_to_camera():
     pointing_radec = SkyCoord.from_name("Crab")
     obstime = Time("2020-01-27T23:00", scale="utc")
-    pointing_alt = u.Quantity(1.3748, u.rad, copy=False)
-    pointing_az = u.Quantity(4.0975, u.rad, copy=False)
+    pointing_alt = u.Quantity(1.3748, u.rad)
+    pointing_az = u.Quantity(4.0975, u.rad)
     focal = 28 * u.m
     expected_source_pos_camera = np.array([0.0, 0.0]) * u.m
     pointing_pos_camera = utils.radec_to_camera(