Upgrade nectarchain to ctapipe 0.19 (#68)

* Adapt package dependencies to prepare for ctapipe-0.19

* Upgrade to ctapipe 0.19.2

* Pin ctapipe_io_nectarcam version to 0.1.3

* Install pinned version of ctapipe-io-nectarcam via PyPI, the time for the conda-forge package to be built.

* Fixed typo

* Put ctapipe_io_nectarcam dependency back into conda dependencies, instead of a pip dependency (was a temporary hack, until the last release of ctapipe_io_nectarcam was published as a conda-forge package).

* Pin ctapipe 0.19.1 (same as for ctapipe_io_nectarcam)

* - minimize the RAM needed by not initializing arrays before reading events in waveformsContainer(s)
 - Update to ctapipe0.19

* Bump ctapipe version into CI tests.

* COMDIRAC is already embedded in DIRAC v8, no need to install it separately

* upgrading DQM to CTAPIPE 0.19.
(Next commit should ake care of flagged bad pixels and how to deal with them)

* upgrading DQM to CTAPIPE 0.19
(Fixing bad pixel handling hsould come next)

* Updated notebook

* bug fic for mean_camera_display.py

* Change the way the python version is specified, otherwise CI tests fail in testing the code for different python versions.

* Python 3.10 interpreted as 3.1 in CI

* clean commit history :
-update SPE fit parameters initialization
following updates of the ctapipe extractors
-improve the limits of generated charge histogram plots

---------

Co-authored-by: Jean-Philippe Lenain <[email protected]>
Co-authored-by: Halim Ashkar <[email protected]>
Co-authored-by: guillaume.grolleron <[email protected]>

.github/workflows/ci.yml

@@ -31,8 +31,8 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.8]
-        ctapipe-version: [v0.12.0]
+        python-version: [3.8, 3.9, "3.10", 3.11]
+        ctapipe-version: [v0.19.1]
 
     defaults:
       run:

README.md

@@ -84,7 +84,7 @@ conda env config vars set X509_CERT_DIR=${CONDA_PREFIX}/etc/grid-security/certif
 # The following is needed for the environment variables, used for DIRAC configuration, to be available:
 mamba deactivate
 mamba activate nectarchain
-pip install CTADIRAC COMDIRAC
+pip install CTADIRAC
 dirac-configure
 ```
 

environment.yml

@@ -3,13 +3,10 @@ channels:
   - default
   - conda-forge
 dependencies:
-  - python=3.8 # nail the python version, so conda does not try upgrading / downgrading
-  - ctapipe=0.12
+  - python=3.11
+  - ctapipe=0.19.1
   - ctapipe-io-nectarcam
-  - numpy=1.22 # higher versions >=1.23 don't work due to astropy4 dependency in ctapipe0.12
   - jupyterlab
-  - protozfits=2.0
-  - pytables>=3.7
   - pytest-runner
   - pip
   - pandas

pyproject.toml

@@ -15,15 +15,17 @@ classifiers = [
     "Topic :: Scientific/Engineering :: Physics",
     "Programming Language :: Python :: 3 :: Only",
     "Programming Language :: Python :: 3.8",
+    "Programming Language :: Python :: 3.9",
+    "Programming Language :: Python :: 3.10",
+    "Programming Language :: Python :: 3.11"
 ]
-requires-python = "~=3.8"
+requires-python = ">=3.8"
 dependencies = [
-    "astropy~=4.2",
-    "ctapipe~=0.12",
+    "browser_cookie3",
+    "ctapipe~=0.19",
     "ctapipe-io-nectarcam",
-    "numpy~=1.22",
-    "protozfits~=2.0",
-    "tables>=3.7",
+    "mechanize",
+    "pandas",
 ]
 dynamic = ["version"]
 

src/nectarchain/calibration/NectarGain/SPEfit/NectarGainSPE.py

@@ -30,6 +30,9 @@
 
 class NectarGainSPE(ABC) :
     _Ncall = 4000000
+    _Windows_lenght = 40
+    _Order = 2
+
     def __init__(self) :
         #set parameters value for fit
         self.__parameters = Parameters()
@@ -104,8 +107,8 @@ def _get_mean_gaussian_fit(charge_in, histo_in ,extension = ""):
         charge = charge_in.data[~histo_in.mask]
         histo = histo_in.data[~histo_in.mask]
 
-        windows_lenght = 80
-        order = 2
+        windows_lenght = NectarGainSPE._Windows_lenght
+        order = NectarGainSPE._Order
         histo_smoothed = savgol_filter(histo, windows_lenght, order)
 
         peaks = find_peaks(histo_smoothed,10)
@@ -126,15 +129,16 @@ def _get_mean_gaussian_fit(charge_in, histo_in ,extension = ""):
 
         if log.getEffectiveLevel() == logging.DEBUG :
             log.debug('plotting figures with prefit parameters computation') 
-            fig,ax = plt.subplots(1,1,figsize = (8,8))
+            fig,ax = plt.subplots(1,1,figsize = (5,5))
             ax.errorbar(charge,histo,np.sqrt(histo),zorder=0,fmt=".",label = "data")
             ax.plot(charge,histo_smoothed,label = f'smoothed data with savgol filter (windows lenght : {windows_lenght}, order : {order})')
             ax.plot(charge,weight_gaussian(charge,coeff_mean[0],coeff_mean[1],coeff_mean[2]),label = 'gaussian fit of the SPE')
             ax.vlines(coeff_mean[1],0,peak_value,label = f'mean initial value = {coeff_mean[1] - coeff[1]:.0f}',color = "red")
             ax.add_patch(Rectangle((coeff_mean[1]-coeff_mean[2], 0), 2 * coeff_mean[2], peak_value_mean,fc=to_rgba('red', 0.5)))
+            ax.set_xlim([peak_pos - 500,None])
             ax.set_xlabel("Charge (ADC)", size=15)
             ax.set_ylabel("Events", size=15)
-            ax.legend(fontsize=15)
+            ax.legend(fontsize=7)
             os.makedirs(f"{os.environ.get('NECTARCHAIN_LOG')}/{os.getpid()}/figures/",exist_ok=True)
             fig.savefig(f"{os.environ.get('NECTARCHAIN_LOG')}/{os.getpid()}/figures/initialization_mean_pixel{extension}_{os.getpid()}.pdf")
             fig.clf()
@@ -149,8 +153,8 @@ def _get_pedestal_gaussian_fit(charge_in, histo_in ,extension = "") :
         charge = charge_in.data[~histo_in.mask]
         histo = histo_in.data[~histo_in.mask]
 
-        windows_lenght = 80
-        order = 2
+        windows_lenght = NectarGainSPE._Windows_lenght
+        order = NectarGainSPE._Order
         histo_smoothed = savgol_filter(histo, windows_lenght, order)
 
         peaks = find_peaks(histo_smoothed,10)
@@ -161,15 +165,16 @@ def _get_pedestal_gaussian_fit(charge_in, histo_in ,extension = "") :
 
         if log.getEffectiveLevel() == logging.DEBUG :
             log.debug('plotting figures with prefit parameters computation') 
-            fig,ax = plt.subplots(1,1,figsize = (8,8))
+            fig,ax = plt.subplots(1,1,figsize = (5,5))
             ax.errorbar(charge,histo,np.sqrt(histo),zorder=0,fmt=".",label = "data")
             ax.plot(charge,histo_smoothed,label = f'smoothed data with savgol filter (windows lenght : {windows_lenght}, order : {order})')
             ax.plot(charge,weight_gaussian(charge,coeff[0],coeff[1],coeff[2]),label = 'gaussian fit of the pedestal, left tail only')
+            ax.set_xlim([peak_pos - 500,None])
             ax.vlines(coeff[1],0,peak_value,label = f'pedestal initial value = {coeff[1]:.0f}',color = 'red')
             ax.add_patch(Rectangle((coeff[1]-coeff[2], 0), 2 * coeff[2], peak_value,fc=to_rgba('red', 0.5)))
             ax.set_xlabel("Charge (ADC)", size=15)
             ax.set_ylabel("Events", size=15)
-            ax.legend(fontsize=15)
+            ax.legend(fontsize=7)
             os.makedirs(f"{os.environ.get('NECTARCHAIN_LOG')}/{os.getpid()}/figures/",exist_ok=True)
             fig.savefig(f"{os.environ.get('NECTARCHAIN_LOG')}/{os.getpid()}/figures/initialization_pedestal_pixel{extension}_{os.getpid()}.pdf")
             fig.clf()

src/nectarchain/calibration/NectarGain/SPEfit/NectarGainSPE_singlerun.py

@@ -457,17 +457,18 @@ def _run_obs_static(cls,it : int, funct,parameters : Parameters, pixels_id : int
 
 
             if kwargs.get('figpath',0) != 0 :
-                fig,ax = plt.subplots(1,1,figsize=(8, 6))
+                fig,ax = plt.subplots(1,1,figsize=(8, 8))
                 ax.errorbar(charge,histo,np.sqrt(histo),zorder=0,fmt=".",label = "data")
                 ax.plot(charge,
                     np.trapz(histo,charge)*MPE2(charge,parameters['pp'].value, parameters['resolution'].value, parameters['mean'].value, parameters['n'].value, parameters['pedestal'].value, parameters['pedestalWidth'].value, parameters['luminosity'].value),
                     zorder=1,
                     linewidth=2,
-                    label = f"SPE model fit \n gain : {gain[0] - gain[1]:.2f} < {gain[0]:.2f} < {gain[0] + gain[2]:.2f} ADC/pe, pvalue = {Statistics.chi2_pvalue(ndof,fit.fval)},\n likelihood = {fit.fval:.2f}")
+                    label = f"SPE model fit \n gain : {gain[0] - gain[1]:.2f} < {gain[0]:.2f} < {gain[0] + gain[2]:.2f} ADC/pe,\n pvalue = {Statistics.chi2_pvalue(ndof,fit.fval)},\n likelihood = {fit.fval:.2f}")
                 ax.set_xlabel("Charge (ADC)", size=15)
                 ax.set_ylabel("Events", size=15)
                 ax.set_title(f"SPE fit pixel {it} with pixel_id : {pixels_id}")
-                ax.legend(fontsize=15)
+                ax.set_xlim([parameters['pedestal'].value - 6 * parameters['pedestalWidth'].value, None])
+                ax.legend(fontsize=18)
                 os.makedirs(kwargs.get('figpath'),exist_ok = True)
                 fig.savefig(f"{kwargs.get('figpath')}/fit_SPE_pixel{pixels_id}.pdf")
                 fig.clf()
@@ -503,17 +504,18 @@ def _run_obs(self,pixel,prescan = False,**kwargs) :
 
 
             if kwargs.get('figpath',0) != 0 :
-                fig,ax = plt.subplots(1,1,figsize=(8, 6))
+                fig,ax = plt.subplots(1,1,figsize=(8, 8))
                 ax.errorbar(self.charge[pixel],self.histo[pixel],np.sqrt(self.histo[pixel]),zorder=0,fmt=".",label = "data")
                 ax.plot(self.charge[pixel],
                     np.trapz(self.histo[pixel],self.charge[pixel])*MPE2(self.charge[pixel],self.__pp.value,self.__resolution.value,self.__mean.value,self.__n.value,self.pedestal.value,self.__pedestalWidth.value,self.__luminosity.value),
                     zorder=1,
                     linewidth=2,
-                    label = f"SPE model fit \n gain : {self.__gain[pixel,0] - self.__gain[pixel,1]:.2f} < {self.__gain[pixel,0]:.2f} < {self.__gain[pixel,0] + self.__gain[pixel,2]:.2f} ADC/pe, pvalue = {Statistics.chi2_pvalue(ndof,fit.fval)},\n likelihood = {fit.fval:.2f}")
+                    label = f"SPE model fit \n gain : {self.__gain[pixel,0] - self.__gain[pixel,1]:.2f} < {self.__gain[pixel,0]:.2f} < {self.__gain[pixel,0] + self.__gain[pixel,2]:.2f} ADC/pe,\n pvalue = {Statistics.chi2_pvalue(ndof,fit.fval)},\n likelihood = {fit.fval:.2f}")
                 ax.set_xlabel("Charge (ADC)", size=15)
                 ax.set_ylabel("Events", size=15)
                 ax.set_title(f"SPE fit pixel : {pixel} (pixel id : {self.pixels_id[pixel]})")
-                ax.legend(fontsize=15)
+                ax.set_xlim([self.pedestal.value - 6 * self.pedestalWidth.value, None])
+                ax.legend(fontsize=18)
                 os.makedirs(kwargs.get('figpath'),exist_ok = True)
                 fig.savefig(f"{kwargs.get('figpath')}/fit_SPE_pixel{self.pixels_id[pixel]}.pdf")
                 fig.clf()
@@ -571,7 +573,7 @@ def _update_parameters_prefit_static(cls,it : int, parameters : Parameters, char
         log.debug(f"pedestalWidth updated : {pedestalWidth}")
         try : 
             coeff,var_matrix =  NectarGainSPE._get_mean_gaussian_fit(charge,histo,f'{it}_nominal')
-            if coeff[1] - pedestal.value < 0 : raise Exception("mean gaussian fit not good")
+            if (coeff[1] - pedestal.value < 0) or ((coeff[1] - coeff[2]) - pedestal.max < 0) : raise Exception("mean gaussian fit not good")
             mean = parameters['mean']
             mean.value = coeff[1] - pedestal.value
             mean.min = (coeff[1] - coeff[2]) - pedestal.max
@@ -785,7 +787,7 @@ def run(self,pixel : int = None,multiproc = False, **kwargs):
         super().run(pixel,multiproc,**kwargs)
 
     def _run_obs(self,pixel,prescan = False,**kwargs) : 
-        if self.__nectarGainSPEresult[pixel]['is_valid'].value : 
+        if self.__nectarGainSPEresult[pixel]['is_valid'] : 
             kwargs['pixel_id'] = self.pixels_id[pixel]
             super()._run_obs(pixel,prescan,**kwargs)
         else :

src/nectarchain/calibration/container/charge.py

@@ -34,6 +34,7 @@
 from numba import guvectorize, float64, int64, bool_
 
 from .waveforms import WaveformsContainer,WaveformsContainers
+from .utils import CtaPipeExtractor
 
 
 
@@ -131,6 +132,7 @@ def __init__(self,charge_hg,charge_lg,peak_hg,peak_lg,run_number,pixels_id,neven
         self.event_id = np.zeros((self.nevents),dtype = np.uint16)
         self.trig_pattern_all = np.zeros((self.nevents,self.CAMERA.n_pixels,4),dtype = bool)
 
+
     @classmethod
     def from_waveforms(cls,waveformContainer : WaveformsContainer,method : str = "FullWaveformSum",**kwargs) : 
         """ create a new ChargeContainer from a WaveformsContainer
@@ -299,11 +301,12 @@ def compute_charge(waveformContainer : WaveformsContainer,channel : int,method :
 
         log.debug(f"Extracting charges with method {method} and extractor_kwargs {extractor_kwargs}")
         ImageExtractor = eval(method)(waveformContainer.subarray,**extractor_kwargs)
+
         if channel == constants.HIGH_GAIN:
-            out = np.array([ImageExtractor(waveformContainer.wfs_hg[i],waveformContainer.TEL_ID,channel) for i in range(len(waveformContainer.wfs_hg))]).transpose(1,0,2)
+            out = np.array([CtaPipeExtractor.get_image_peak_time(ImageExtractor(waveformContainer.wfs_hg[i],waveformContainer.TEL_ID,channel,waveformContainer.broken_pixels_hg)) for i in range(len(waveformContainer.wfs_hg))]).transpose(1,0,2)
             return out[0],out[1]
         elif channel == constants.LOW_GAIN:
-            out = np.array([ImageExtractor(waveformContainer.wfs_lg[i],waveformContainer.TEL_ID,channel) for i in range(len(waveformContainer.wfs_lg))]).transpose(1,0,2)
+            out = np.array([CtaPipeExtractor.get_image_peak_time(ImageExtractor(waveformContainer.wfs_lg[i],waveformContainer.TEL_ID,channel,waveformContainer.broken_pixels_lg)) for i in range(len(waveformContainer.wfs_lg))]).transpose(1,0,2)
             return out[0],out[1]
         else :
             raise ArgumentError(f"channel must be {constants.LOW_GAIN} or {constants.HIGH_GAIN}")

src/nectarchain/calibration/container/utils.py

@@ -10,6 +10,8 @@
 from pathlib import Path
 from typing import List,Tuple
 
+from ctapipe.containers import DL1CameraContainer 
+
 
 import logging
 logging.basicConfig(format='%(asctime)s %(name)s %(levelname)s %(message)s')
@@ -180,3 +182,6 @@ def chainEventSource(list : list,max_events : int = None) : #useless with ctapip
         else :
             return ChainGenerator.chain(list[0],ChainGenerator.chainEventSource(list[1:]))
 
+class CtaPipeExtractor():
+    def get_image_peak_time(cameraContainer : DL1CameraContainer) : 
+        return cameraContainer.image, cameraContainer.peak_time

src/nectarchain/calibration/container/waveforms.py

@@ -40,7 +40,7 @@ def __new__(cls,*args,**kwargs) :
         obj = object.__new__(cls)
         return obj
 
-    def __init__(self,run_number : int,max_events : int = None,nevents : int = -1,run_file = None):
+    def __init__(self,run_number : int,max_events : int = None,nevents : int = -1,run_file = None, init_arrays : bool = False):
         """construtor
 
         Args:
@@ -52,7 +52,9 @@ def __init__(self,run_number : int,max_events : int = None,nevents : int = -1,ru
         """
 
         self.__run_number = run_number
-        #gerer ici le fait de traiter plusieurs fichiers ou simplement 1 par 1
+        self.__run_file = run_file
+        self.__max_events = max_events
+
         self.__reader = WaveformsContainer.load_run(run_number,max_events,run_file = run_file)
 
         #from reader members
@@ -68,12 +70,20 @@ def __init__(self,run_number : int,max_events : int = None,nevents : int = -1,ru
         #run properties
         if nevents != -1 :
             self.__nevents = nevents if max_events is None else min(max_events,nevents) #self.check_events()
+            self.__reader = None
         else :
             self.__nevents = self.check_events()
-            #reload file (bc check_events has drained reader generator)
-            self.__reader = WaveformsContainer.load_run(run_number,max_events,run_file = run_file)
+
         log.info(f"N_events : {self.nevents}")
 
+        if init_arrays : 
+            self.__init_arrays()
+
+    def __init_arrays(self,**kwargs) : 
+        log.debug('creation of the EventSource reader')
+        self.__reader = WaveformsContainer.load_run(self.__run_number,self.__max_events,run_file = self.__run_file)
+
+        log.debug("create wfs, ucts, event properties and triger pattern arrays")
         #define zeros members which will be filled therafter
         self.wfs_hg = np.zeros((self.nevents,self.npixels,self.nsamples),dtype = np.uint16)
         self.wfs_lg = np.zeros((self.nevents,self.npixels,self.nsamples),dtype = np.uint16)
@@ -86,6 +96,14 @@ def __init__(self,run_number : int,max_events : int = None,nevents : int = -1,ru
         #self.trig_pattern = np.zeros((self.nevents,self.npixels),dtype = bool)
         #self.multiplicity = np.zeros((self.nevents,self.npixels),dtype = np.uint16)
 
+        self.__broken_pixels_hg = np.zeros((self.npixels),dtype = bool)
+        self.__broken_pixels_lg = np.zeros((self.npixels),dtype = bool)
+
+
+    def __compute_broken_pixels(self,**kwargs) : 
+        log.warning("computation of broken pixels is not yet implemented")
+        self.__broken_pixels_hg = np.zeros((self.npixels),dtype = bool)
+        self.__broken_pixels_lg = np.zeros((self.npixels),dtype = bool)
 
     @staticmethod
     def load_run(run_number : int,max_events : int = None, run_file = None) : 
@@ -98,8 +116,8 @@ def load_run(run_number : int,max_events : int = None, run_file = None) :
         Returns:
             List[ctapipe_io_nectarcam.NectarCAMEventSource]: List of EventSource for each run files
         """
-        generic_filename,filenames = DataManagement.findrun(run_number)
         if run_file is None : 
+            generic_filename,_ = DataManagement.findrun(run_number)
             log.info(f"{str(generic_filename)} will be loaded")
             eventsource = NectarCAMEventSource(input_url=generic_filename,max_events=max_events)
         else :  
@@ -135,6 +153,9 @@ def load_wfs(self,compute_trigger_patern = False):
         Args:
             compute_trigger_patern (bool, optional): To recompute on our side the trigger patern. Defaults to False.
         """
+        if not(hasattr(self, "wfs_hg")) : 
+            self.__init_arrays()
+
         wfs_hg_tmp=np.zeros((self.npixels,self.nsamples),dtype = np.uint16)
         wfs_lg_tmp=np.zeros((self.npixels,self.nsamples),dtype = np.uint16)
 
@@ -159,7 +180,7 @@ def load_wfs(self,compute_trigger_patern = False):
             self.wfs_hg[i] = wfs_hg_tmp
             self.wfs_lg[i] = wfs_lg_tmp
 
-
+        self.__compute_broken_pixels()
 
         #if compute_trigger_patern and np.max(self.trig_pattern) == 0:
         #    self.compute_trigger_patern()
@@ -262,6 +283,8 @@ def load(path : str) :
             table_trigger = hdul[4].data
             cls.trig_pattern_all = table_trigger["trig_pattern_all"]
 
+        cls.__compute_broken_pixels()
+
         return cls
 
 
@@ -345,6 +368,13 @@ def select_waveforms_lg(self,pixel_id : np.ndarray) :
         res = res.transpose(res.shape[1],res.shape[0],res.shape[2])
         return res
 
+
+    @property 
+    def _run_file(self) : return self.__run_file 
+
+    @property
+    def _max_events(self) : return self.__max_events
+
     @property
     def reader(self) : return self.__reader
 
@@ -369,7 +399,12 @@ def nevents(self) : return self.__nevents
     @property
     def run_number(self) : return self.__run_number
 
+    @property
+    def broken_pixels_hg(self) : return self.__broken_pixels_hg
 
+    @property
+    def broken_pixels_lg(self) : return self.__broken_pixels_lg
+
     #physical properties
     @property
     def multiplicity(self) :  return np.uint16(np.count_nonzero(self.trig_pattern,axis = 1))
@@ -395,7 +430,7 @@ def __new__(cls,*args,**kwargs) :
         obj = object.__new__(cls)
         return obj
 
-    def __init__(self,run_number : int,max_events : int = None) :
+    def __init__(self,run_number : int,max_events : int = None, init_arrays : bool = False) :
         """initialize the waveformsContainer list inside the main object
 
         Args:
@@ -407,7 +442,7 @@ def __init__(self,run_number : int,max_events : int = None) :
         self.waveformsContainer = []
         self.__nWaveformsContainer = 0
         for i,file in enumerate(filenames) : 
-            self.waveformsContainer.append(WaveformsContainer(run_number,max_events=max_events,run_file=file))
+            self.waveformsContainer.append(WaveformsContainer(run_number,max_events=max_events,run_file=file, init_arrays= init_arrays))
             self.__nWaveformsContainer += 1
             if not(max_events is None) : max_events -= self.waveformsContainer[i].nevents
             log.info(f'WaveformsContainer number {i} is created')