diff --git a/.gitignore b/.gitignore index bb7534aa..ce9c8e3d 100644 --- a/.gitignore +++ b/.gitignore @@ -92,6 +92,7 @@ src/nectarchain/user_scripts/**/test **.log **.pdf **.csv +**.png #VScode .vscode/ diff --git a/notebooks/tool_implementation/tuto_SPE.py b/notebooks/tool_implementation/tuto_SPE.py new file mode 100644 index 00000000..f35c9c1f --- /dev/null +++ b/notebooks/tool_implementation/tuto_SPE.py @@ -0,0 +1,75 @@ +# --- +# jupyter: +# jupytext: +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.14.6 +# kernelspec: +# display_name: nectarchain +# language: python +# name: python3 +# --- + +# %% +import logging +import os +import pathlib + +logging.basicConfig( + format="%(asctime)s %(name)s %(levelname)s %(message)s", level=logging.INFO +) +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + +from nectarchain.makers.calibration import ( + FlatFieldSPEHHVStdNectarCAMCalibrationTool, + FlatFieldSPENominalStdNectarCAMCalibrationTool, +) + +# %% +run_number = 3936 + +# %% + +# !ls -lh $NECTARCAMDATA/runs/* + +# %% +tool = FlatFieldSPENominalStdNectarCAMCalibrationTool( + progress_bar=True, + method="LocalPeakWindowSum", + extractor_kwargs={"window_width": 12, "window_shift": 4}, + multiproc=True, + nproc=10, + run_number=run_number, + max_events=10000, + log_level=20, + reload_events=True, + # events_per_slice = 200, + asked_pixels_id=[52, 48], + output_path=pathlib.Path(os.environ.get("NECTARCAMDATA", "/tmp")) + / "tutorials/" + / f"SPEfit_{run_number}.h5", +) + +# %% +tool + +# %% +tool.initialize() + +# %% +tool.setup() + +# %% +tool.start() + +# %% +output = tool.finish(return_output_component=True, display=True, figpath=os.getcwd()) +output + +# %% +output[0].resolution + +# %% diff --git a/notebooks/tool_implementation/tuto_tools.py b/notebooks/tool_implementation/tuto_tools.py new file mode 100644 index 00000000..595e4a11 --- /dev/null +++ b/notebooks/tool_implementation/tuto_tools.py @@ -0,0 +1,389 @@ +# --- +# jupyter: +# jupytext: +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.14.6 +# kernelspec: +# display_name: nectarchain +# language: python +# name: python3 +# --- + +# %% [markdown] +# # Tutorial to use ctapipe tools and component + +# %% [markdown] +# ## Context + +# %% [markdown] +# Tool and Component are 2 modules of ctapipe, Tool is a high level module to analyse raw data (fits.fz files). This module use Component to perform computation on the raw data. Basically, we can create a class (MyTool) which inherits of Tool, where we can define 2 Component (Comp_A and Comp_B). Thus, with an instance of MyTool, we can loop over event within raw data, and for each event apply sucessively Comp_A, Comp_B. +# +# A ctapipe tutorial is accessible here : https://ctapipe.readthedocs.io/en/stable/auto_examples/core/command_line_tools.html#sphx-glr-auto-examples-core-command-line-tools-py +# +# You can find documentation of ctapipe Tool and Component : +# +# https://ctapipe.readthedocs.io/en/stable/api-reference/tools/index.html +# +# +# https://ctapipe.readthedocs.io/en/stable/api/ctapipe.core.Component.html +# +# Within nectarchain, we implemented within the nectarchain.makers module both a top level Tool and Component from which all the nectarchain Component and Tool should inherit. +# +# In this tutorial, we explain quickly how we can use Tool and Component to develop the nectarchain software, as an example there is the implementation of a PedestalTool which extract pedestal + +# %% [markdown] +# ### Imports + +# %% +import numpy as np +import pathlib +import os + +import matplotlib.pyplot as plt + +from ctapipe_io_nectarcam.containers import NectarCAMDataContainer +from ctapipe_io_nectarcam import constants + +from ctapipe.core.traits import ComponentNameList +from ctapipe.containers import Field +from ctapipe.core.traits import Integer +from ctapipe.core import Component +from ctapipe.io import HDF5TableReader + + +from nectarchain.makers import EventsLoopNectarCAMCalibrationTool +from nectarchain.makers.component import NectarCAMComponent, ArrayDataComponent +from nectarchain.data.container import ( + NectarCAMContainer, + ArrayDataContainer, + TriggerMapContainer, +) +from nectarchain.utils import ComponentUtils + +# %% +tools = EventsLoopNectarCAMCalibrationTool() +tools + +# %% +tools.classes + + +# %% [markdown] +# The only thing to add to to fill the componentList field, which contains the names of the component to be apply on events. +# +# Then we will define a very simple component to compute the pedestal of each events. + +# %% [markdown] +# ### Definition of container to store extracted data on disk + +# %% +class MyContainer(NectarCAMContainer): + run_number = Field( + type=np.uint16, + description="run number associated to the waveforms", + ) + npixels = Field( + type=np.uint16, + description="number of effective pixels", + ) + pixels_id = Field(type=np.ndarray, dtype=np.uint16, ndim=1, description="pixel ids") + ucts_timestamp = Field( + type=np.ndarray, dtype=np.uint64, ndim=1, description="events ucts timestamp" + ) + event_type = Field( + type=np.ndarray, dtype=np.uint8, ndim=1, description="trigger event type" + ) + event_id = Field(type=np.ndarray, dtype=np.uint32, ndim=1, description="event ids") + + pedestal_hg = Field( + type=np.ndarray, dtype=np.uint16, ndim=2, description="The high gain pedestal" + ) + pedestal_lg = Field( + type=np.ndarray, dtype=np.uint16, ndim=2, description="The low gain pedestal" + ) + + +# %% [markdown] +# ### Definition of our Component + +# %% +class MyComp(NectarCAMComponent): + window_shift = Integer( + default_value=4, + help="the time in ns before the peak to extract charge", + ).tag(config=True) + + window_width = Integer( + default_value=12, + help="the duration of the extraction window in ns", + ).tag(config=True) + + def __init__(self, subarray, config=None, parent=None, *args, **kwargs): + super().__init__( + subarray=subarray, config=config, parent=parent, *args, **kwargs + ) + ## If you want you can add here members of MyComp, they will contain interesting quantity during the event loop process + + self.__ucts_timestamp = [] + self.__event_type = [] + self.__event_id = [] + + self.__pedestal_hg = [] + self.__pedestal_lg = [] + + ##This method need to be defined ! + def __call__(self, event: NectarCAMDataContainer, *args, **kwargs): + self.__event_id.append(np.uint32(event.index.event_id)) + self.__event_type.append(event.trigger.event_type.value) + self.__ucts_timestamp.append(event.nectarcam.tel[0].evt.ucts_timestamp) + + wfs = [] + wfs.append(event.r0.tel[0].waveform[constants.HIGH_GAIN][self.pixels_id]) + wfs.append(event.r0.tel[0].waveform[constants.LOW_GAIN][self.pixels_id]) + + #####THE JOB IS HERE###### + for i, pedestal in enumerate([self.__pedestal_hg, self.__pedestal_lg]): + index_peak = np.argmax(wfs[i]) + signal_start = index_peak - self.window_shiftclasses + signal_stop = index_peak + self.window_width - self.window_shift + if signal_start < 0: + signal_stop = self.window_width + signal_start = 0 + if signal_stop > constants.N_SAMPLES: + signal_stop = constants.N_SAMPLES + signal_start = constants.N_SAMPLES - self.window_width + pedestal.append( + ( + np.sum(wfs[i][:, 0:signal_start], axis=1) + + np.sum(wfs[i][:, signal_stop:], axis=1) + ) + / (constants.N_SAMPLES - self.window_width) + ) + + ##This method need to be defined ! + def finish(self): + output = MyContainer( + run_number=MyContainer.fields["run_number"].type(self._run_number), + npixels=MyContainer.fields["npixels"].type(self._npixels), + pixels_id=MyContainer.fields["pixels_id"].dtype.type(self._pixels_id), + ucts_timestamp=MyContainer.fields["ucts_timestamp"].dtype.type( + self.__ucts_timestamp + ), + event_type=MyContainer.fields["event_type"].dtype.type(self.__event_type), + event_id=MyContainer.fields["event_id"].dtype.type(self.__event_id), + pedestal_hg=MyContainer.fields["pedestal_hg"].dtype.type( + self.__pedestal_hg + ), + pedestal_lg=MyContainer.fields["pedestal_lg"].dtype.type( + self.__pedestal_lg + ), + ) + return output + + +# %% [markdown] +# ### Definition of our Tool + +# %% [markdown] +# Now we can define out Tool, we have just to add our component "MyComp" in the ComponentList : + +# %% +def get_valid_component(): + return NectarCAMComponent.non_abstract_subclasses() + + +class MyTool(EventsLoopNectarCAMCalibrationTool): + name = "PedestalTutoNectarCAM" + + componentsList = ComponentNameList( + NectarCAMComponent, + default_value=["MyComp"], + help="List of Component names to be apply, the order will be respected", + ).tag(config=True) + + ####### THIS PART IS NEEDED NOW ##### + + # This part uis needed because the component defined outside of the nectarchain + # module is not accessible from the module (with eval(componentName)). + # This issus should be fixed soon. + def __new__(cls, *args, **kwargs): + _cls = super(EventsLoopNectarCAMCalibrationTool, cls).__new__( + cls, *args, **kwargs + ) + + for componentName in _cls.componentsList: + configurable_traits = ComponentUtils.get_configurable_traits( + eval(componentName) + ) + _cls.add_traits(**configurable_traits) + _cls.aliases.update( + {key: f"{componentName}.{key}" for key in configurable_traits.keys()} + ) + return _cls + + def _get_provided_component_kwargs(self, componentName: str): + component_kwargs = ComponentUtils.get_configurable_traits(eval(componentName)) + output_component_kwargs = {} + for key in component_kwargs.keys(): + if hasattr(self, key): + output_component_kwargs[key] = getattr(self, key) + return output_component_kwargs + + ##################################################################### + + # def __init__(self,*args,**kwargs) : + # super().__init__(*args,**kwargs) + + def _init_output_path(self): + if self.max_events is None: + filename = f"{self.name}_run{self.run_number}.h5" + else: + filename = f"{self.name}_run{self.run_number}_maxevents{self.max_events}.h5" + self.output_path = pathlib.Path( + f"{os.environ.get('NECTARCAMDATA','/tmp')}/tutorials/{filename}" + ) + + +# %% +tool = MyTool( + progress_bar=True, run_number=4943, max_events=500, log_level=20, window_width=14 +) + +# %% +tool.componentsList + +# %% +tool + +# %% [markdown] +# First we have to initialize the tool : + +# %% +tool.initialize() + +# %% [markdown] +# Then to setup, it will in particular setup the Components : + +# %% +tool.setup() + +# %% [markdown] +# The following command will just start the tool and apply components looping over events + +# %% +tool.start() + +# %% [markdown] +# Then, we finish the tool, behind thius command the component will be finilized and will create an output container whiich will be written on disk and can be returned + +# %% +output = tool.finish(return_output_component=True)[0] + +# %% +output + +# %% [markdown] +# The following file has been written : + +# %% +# !ls -lh $NECTARCAMDATA/tutorials + +# %% [markdown] +# The shape of pedestal is (n_events,n_pixels) + +# %% +output.pedestal_hg.shape + +# %% [markdown] +# To have a look to a random pixel pedestal evolution : + +# %% +fix, ax = plt.subplots(1, 1) +i = np.random.randint(0, len(output.pixels_id)) + +ax.plot( + (output.ucts_timestamp - output.ucts_timestamp[0]) / 1e6, + output.pedestal_hg[:, i], + linewidth=0.5, +) +ax.set_ylabel("Charge [ADC]") +ax.set_xlabel("time [ms]") + +# %% [markdown] +# If you want to load container thereafter : + +# %% +container_loaded = MyContainer._container_from_hdf5( + f"{os.environ.get('NECTARCAMDATA','/tmp')}/tutorials/PedestalTutoNectarCAM_run4943_maxevents500.h5", + MyContainer, +) +container_loaded.validate() +container_loaded + +# %% [markdown] +# ## Going further + +# %% [markdown] +# An argument that are implemented in EventsLoopNectarCAMCalibrationTool is the 'event_per_slice', this argument allows to split all the events within the raw data fits.fz file in slices. It allows to, for each slice, loop over events and write container on disk. This mechanism allows to save RAM. +# The resulting hdf5 file that is written on disk , can be easily loaded thereafter. There is only one hdf5 file for the whole run, which is a mapping between slices and containers filled by computed quantity from components. + +# %% +tool = MyTool( + progress_bar=True, + run_number=4943, + max_events=2000, + log_level=20, + events_per_slice=1000, + overwrite=True, +) +tool + +# %% +tool.initialize() +tool.setup() + +# %% +tool.start() + +# %% +output = tool.finish(return_output_component=True)[0] +output + + +# %% +# !h5ls -r $NECTARCAMDATA/tutorials/PedestalTutoNectarCAM_run4943_maxevents2000.h5 + +# %% +# container_loaded = ArrayDataContainer._container_from_hdf5(f"{os.environ.get('NECTARCAMDATA','/tmp')}/tutorials/PedestalTutoNectarCAM_run4943_maxevents2000.h5",MyContainer) +# container_loaded + +# %% +def read_hdf5_sliced(path): + container = MyContainer() + container_class = MyContainer + with HDF5TableReader(path) as reader: + for data in reader._h5file.root.__members__: + # print(data) + data_cont = eval(f"reader._h5file.root.{data}.__members__")[0] + # print(data_cont) + tableReader = reader.read( + table_name=f"/{data}/{data_cont}", containers=container_class + ) + # container.containers[data].containers[trigger] = next(tableReader) + container = next(tableReader) + yield container + + +# %% +container_loaded = read_hdf5_sliced( + f"{os.environ.get('NECTARCAMDATA','/tmp')}/tutorials/PedestalTutoNectarCAM_run4943_maxevents2000.h5" +) +for i, container in enumerate(container_loaded): + print( + f"Container {i} is filled by events from {container.event_id[0]} to {container.event_id[-1]}" + ) + +# %% diff --git a/notebooks/tool_implementation/tuto_waveforms_charges_extraction.py b/notebooks/tool_implementation/tuto_waveforms_charges_extraction.py new file mode 100644 index 00000000..b77b00ea --- /dev/null +++ b/notebooks/tool_implementation/tuto_waveforms_charges_extraction.py @@ -0,0 +1,69 @@ +# --- +# jupyter: +# jupytext: +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.14.6 +# kernelspec: +# display_name: nectarchain +# language: python +# name: python3 +# --- + +# %% +from nectarchain.makers import ( + ChargesNectarCAMCalibrationTool, + WaveformsNectarCAMCalibrationTool, +) +from nectarchain.makers.component import get_valid_component +from nectarchain.data.container import ( + ChargesContainers, + ChargesContainer, + WaveformsContainer, + WaveformsContainers, +) +from ctapipe.io import HDF5TableReader +from ctapipe.containers import EventType + +# %% +get_valid_component() + +# %% +run_number = 3942 + +# %% +tool = WaveformsNectarCAMCalibrationTool( + progress_bar=True, run_number=run_number, max_events=500, log_level=20 +) +tool + +# %% +tool = ChargesNectarCAMCalibrationTool( + progress_bar=True, + method="LocalPeakWindowSum", + extractor_kwargs={"window_width": 12, "window_shift": 4}, + run_number=run_number, + max_events=500, + log_level=20, +) +tool + +# %% +tool.initialize() + +# %% +tool.setup() + +# %% +tool.start() + +# %% +output = tool.finish(return_output_component=True)[0] +output + +# %% +output.containers + +# %% diff --git a/src/nectarchain/data/container/__init__.py b/src/nectarchain/data/container/__init__.py index e69de29b..464a8848 100644 --- a/src/nectarchain/data/container/__init__.py +++ b/src/nectarchain/data/container/__init__.py @@ -0,0 +1,11 @@ +from .chargesContainer import * +from .core import ( + ArrayDataContainer, + NectarCAMContainer, + TriggerMapContainer, + get_array_keys, + merge_map_ArrayDataContainer, +) +from .eventSource import * +from .gainContainer import * +from .waveformsContainer import * diff --git a/src/nectarchain/data/container/charges_container.py b/src/nectarchain/data/container/chargesContainer.py similarity index 91% rename from src/nectarchain/data/container/charges_container.py rename to src/nectarchain/data/container/chargesContainer.py index c00d1b5d..78929e13 100644 --- a/src/nectarchain/data/container/charges_container.py +++ b/src/nectarchain/data/container/chargesContainer.py @@ -1,13 +1,15 @@ import logging +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + import numpy as np from ctapipe.containers import Field, Map, partial from .core import ArrayDataContainer, TriggerMapContainer -logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") -log = logging.getLogger(__name__) -log.handlers = logging.getLogger("__main__").handlers +__all__ = ["ChargesContainer", "ChargesContainers"] class ChargesContainer(ArrayDataContainer): @@ -40,5 +42,5 @@ class ChargesContainer(ArrayDataContainer): class ChargesContainers(TriggerMapContainer): containers = Field( default_factory=partial(Map, ChargesContainer), - description="trigger mapping of ChargesContainer", + description="trigger or slices of data mapping of ChargesContainer", ) diff --git a/src/nectarchain/data/container/core.py b/src/nectarchain/data/container/core.py index 458ee396..48fe5b8e 100644 --- a/src/nectarchain/data/container/core.py +++ b/src/nectarchain/data/container/core.py @@ -1,20 +1,54 @@ import logging -import numpy as np -from ctapipe.containers import Container, Field, Map, partial - logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") log = logging.getLogger(__name__) log.handlers = logging.getLogger("__main__").handlers +import copy +import importlib +from pathlib import Path + +import numpy as np +from ctapipe.containers import Container, EventType, Field, Map, partial +from ctapipe.core.container import FieldValidationError +from ctapipe.io import HDF5TableReader +from tables.exceptions import NoSuchNodeError + +__all__ = [ + "ArrayDataContainer", + "TriggerMapContainer", + "get_array_keys", + "merge_map_ArrayDataContainer", +] -__all__ = ["ArrayDataContainer", "TriggerMapContainer"] + +def get_array_keys(container: Container): + keys = [] + for field in container.fields: + if field.type == np.ndarray: + keys.append(field.key) + return keys class NectarCAMContainer(Container): """base class for the NectarCAM containers. This contaner cannot berecursive, to be directly written with a HDF5TableWriter""" + @staticmethod + def _container_from_hdf5(path, container_class): + if isinstance(path, str): + path = Path(path) + + container = container_class() + with HDF5TableReader(path) as reader: + tableReader = reader.read( + table_name=f"/data/{container_class.__name__}", + containers=container_class, + ) + container = next(tableReader) + + yield container + class ArrayDataContainer(NectarCAMContainer): """ @@ -39,15 +73,15 @@ class ArrayDataContainer(NectarCAMContainer): """ run_number = Field( - type=int, + type=np.uint16, description="run number associated to the waveforms", ) nevents = Field( - type=int, + type=np.uint64, description="number of events", ) npixels = Field( - type=int, + type=np.uint16, description="number of effective pixels", ) pixels_id = Field(type=np.ndarray, dtype=np.uint16, ndim=1, description="pixel ids") @@ -84,9 +118,129 @@ class ArrayDataContainer(NectarCAMContainer): type=np.ndarray, dtype=np.uint16, ndim=1, description="events multiplicity" ) + @staticmethod + def _container_from_hdf5(path, container_class, slice_index=None): + if isinstance(path, str): + path = Path(path) + module = importlib.import_module(f"{container_class.__module__}") + container = eval(f"module.{container_class.__name__}s")() + + with HDF5TableReader(path) as reader: + if len(reader._h5file.root.__members__) > 1 and slice_index is None: + log.info( + f"reading {container_class.__name__}s containing {len(reader._h5file.root.__members__)} slices, will return a generator" + ) + for data in reader._h5file.root.__members__: + # container.containers[data] = eval(f"module.{container_class.__name__}s")() + for key, trigger in EventType.__members__.items(): + try: + waveforms_data = eval( + f"reader._h5file.root.{data}.__members__" + ) + _mask = [ + container_class.__name__ in _word + for _word in waveforms_data + ] + _waveforms_data = np.array(waveforms_data)[_mask] + if len(_waveforms_data) == 1: + tableReader = reader.read( + table_name=f"/{data}/{_waveforms_data[0]}/{trigger.name}", + containers=container_class, + ) + # container.containers[data].containers[trigger] = next(tableReader) + container.containers[trigger] = next(tableReader) + + else: + log.info( + f"there is {len(_waveforms_data)} entry corresponding to a {container_class} table save, unable to load" + ) + except NoSuchNodeError as err: + log.warning(err) + except Exception as err: + log.error(err, exc_info=True) + raise err + yield container + else: + if slice_index is None: + log.info( + f"reading {container_class.__name__}s containing a single slice, will return the {container_class.__name__}s instance" + ) + data = "data" + else: + log.info( + f"reading slice {slice_index} of {container_class.__name__}s, will return the {container_class.__name__}s instance" + ) + data = f"data_{slice_index}" + for key, trigger in EventType.__members__.items(): + try: + container_data = eval(f"reader._h5file.root.{data}.__members__") + _mask = [ + container_class.__name__ in _word + for _word in container_data + ] + _container_data = np.array(container_data)[_mask] + if len(_container_data) == 1: + tableReader = reader.read( + table_name=f"/{data}/{_container_data[0]}/{trigger.name}", + containers=container_class, + ) + container.containers[trigger] = next(tableReader) + else: + log.info( + f"there is {len(_container_data)} entry corresponding to a {container_class} table save, unable to load" + ) + except NoSuchNodeError as err: + log.warning(err) + except Exception as err: + log.error(err, exc_info=True) + raise err + yield container + return container + + @classmethod + def from_hdf5(cls, path, slice_index=None): + return cls._container_from_hdf5( + path, slice_index=slice_index, container_class=cls + ) + class TriggerMapContainer(Container): containers = Field( default_factory=partial(Map, Container), description="trigger mapping of Container", ) + + def is_empty(self): + return len(self.containers.keys()) == 0 + + def validate(self): + super().validate() + for i, container in enumerate(self.containers): + if i == 0: + container_type = type(container) + else: + if not (isinstance(container, container_type)): + raise FieldValidationError( + "all the containers mapped must have the same type to be merged " + ) + + +def merge_map_ArrayDataContainer(triggerMapContainer: TriggerMapContainer): + triggerMapContainer.validate() + keys = list(triggerMapContainer.containers.keys()) + output_container = copy.deepcopy(triggerMapContainer.containers[keys[0]]) + for key in keys[1:]: + for field in get_array_keys(triggerMapContainer.containers[key]): + output_container[field] = np.concatenate( + (output_container[field], triggerMapContainer.containers[key][field]), + axis=0, + ) + if "nevents" in output_container.fields: + output_container.nevents += triggerMapContainer.containers[key].nevents + return output_container + + +# class TriggerMapArrayDataContainer(TriggerMapContainer): +# containers = Field(default_factory=partial(Map, ArrayDataContainer), +# description = "trigger mapping of arrayDataContainer" +# ) diff --git a/src/nectarchain/data/container/eventSource.py b/src/nectarchain/data/container/eventSource.py new file mode 100644 index 00000000..0706b1b9 --- /dev/null +++ b/src/nectarchain/data/container/eventSource.py @@ -0,0 +1,260 @@ +import logging +import sys + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + + +import struct + +import numpy as np +from ctapipe.containers import EventType +from ctapipe_io_nectarcam import ( + NectarCAMDataContainer, + NectarCAMEventSource, + TriggerBits, + time_from_unix_tai_ns, +) +from ctapipe_io_nectarcam.anyarray_dtypes import ( + CDTS_AFTER_37201_DTYPE, + CDTS_BEFORE_37201_DTYPE, + TIB_DTYPE, +) +from ctapipe_io_nectarcam.constants import N_GAINS, N_PIXELS, N_SAMPLES +from ctapipe_io_nectarcam.containers import NectarCAMEventContainer + +__all__ = ["LightNectarCAMEventSource"] + +# DONE event.trigger.event_type +# DONE event.index.event_id +# DONE event.nectarcam.tel[__class__.TEL_ID.default_value].evt.ucts_timestamp +# DONE event.nectarcam.tel[__class__.TEL_ID.default_value].evt.ucts_busy_counter +# DONE event.nectarcam.tel[__class__.TEL_ID.default_value].evt.ucts_event_counter +# DONE event.nectarcam.tel[__class__.TEL_ID.default_value].evt.trigger_pattern +# DONE event.r0.tel[0].waveform + + +def fill_nectarcam_event_container_from_zfile(self, array_event, event): + tel_id = self.tel_id + event_container = NectarCAMEventContainer() + array_event.nectarcam.tel[tel_id].evt = event_container + # event_container.configuration_id = event.configuration_id + # event_container.event_id = event.event_id + # event_container.tel_event_id = event.tel_event_id + # event_container.pixel_status = event.pixel_status + # event_container.ped_id = event.ped_id + # event_container.module_status = event.nectarcam.module_status + event_container.extdevices_presence = event.nectarcam.extdevices_presence + # event_container.swat_data = event.nectarcam.swat_data + event_container.counters = event.nectarcam.counters + ## unpack TIB data + unpacked_tib = event.nectarcam.tib_data.view(TIB_DTYPE)[0] + # event_container.tib_event_counter = unpacked_tib[0] + # event_container.tib_pps_counter = unpacked_tib[1] + # event_container.tib_tenMHz_counter = unpacked_tib[2] + # event_container.tib_stereo_pattern = unpacked_tib[3] + event_container.tib_masked_trigger = unpacked_tib[4] + # unpack CDTS data + is_old_cdts = len(event.nectarcam.cdts_data) < 36 + if is_old_cdts: + unpacked_cdts = event.nectarcam.cdts_data.view(CDTS_BEFORE_37201_DTYPE)[0] + event_container.ucts_event_counter = unpacked_cdts[0] + # event_container.ucts_pps_counter = unpacked_cdts[1] + # event_container.ucts_clock_counter = unpacked_cdts[2] + event_container.ucts_timestamp = unpacked_cdts[3] + # event_container.ucts_camera_timestamp = unpacked_cdts[4] + event_container.ucts_trigger_type = unpacked_cdts[5] + # event_container.ucts_white_rabbit_status = unpacked_cdts[6] + else: + unpacked_cdts = event.nectarcam.cdts_data.view(CDTS_AFTER_37201_DTYPE)[0] + event_container.ucts_timestamp = unpacked_cdts[0] + # event_container.ucts_address = unpacked_cdts[1] # new + event_container.ucts_event_counter = unpacked_cdts[2] + event_container.ucts_busy_counter = unpacked_cdts[3] # new + # event_container.ucts_pps_counter = unpacked_cdts[4] + # event_container.ucts_clock_counter = unpacked_cdts[5] + event_container.ucts_trigger_type = unpacked_cdts[6] + # event_container.ucts_white_rabbit_status = unpacked_cdts[7] + # event_container.ucts_stereo_pattern = unpacked_cdts[8] # new + # event_container.ucts_num_in_bunch = unpacked_cdts[9] # new + # event_container.cdts_version = unpacked_cdts[10] # new + # Unpack FEB counters and trigger pattern + self.unpack_feb_data(event_container, event) + + +def unpack_feb_data(self, event_container, event): + """Unpack FEB counters and trigger pattern""" + # Deduce data format version + bytes_per_module = ( + len(event.nectarcam.counters) // self.camera_config.nectarcam.num_modules + ) + # Remain compatible with data before addition of trigger pattern + module_fmt = "IHHIBBBBBBBB" if bytes_per_module > 16 else "IHHIBBBB" + n_fields = len(module_fmt) + rec_fmt = "=" + module_fmt * self.camera_config.nectarcam.num_modules + # Unpack + unpacked_feb = struct.unpack(rec_fmt, event.nectarcam.counters) + # Initialize field containers + # n_camera_modules = N_PIXELS // 7 + # event_container.feb_abs_event_id = np.zeros(shape=(n_camera_modules,), dtype=np.uint32) + # event_container.feb_event_id = np.zeros(shape=(n_camera_modules,), dtype=np.uint16) + # event_container.feb_pps_cnt = np.zeros(shape=(n_camera_modules,), dtype=np.uint16) + # event_container.feb_ts1 = np.zeros(shape=(n_camera_modules,), dtype=np.uint32) + # event_container.feb_ts2_trig = np.zeros(shape=(n_camera_modules,), dtype=np.int16) + # event_container.feb_ts2_pps = np.zeros(shape=(n_camera_modules,), dtype=np.int16) + if bytes_per_module > 16: + n_patterns = 4 + event_container.trigger_pattern = np.zeros( + shape=(n_patterns, N_PIXELS), dtype=bool + ) + # Unpack absolute event ID + # event_container.feb_abs_event_id[ + # self.camera_config.nectarcam.expected_modules_id] = unpacked_feb[0::n_fields] + ## Unpack PPS counter + # event_container.feb_pps_cnt[ + # self.camera_config.nectarcam.expected_modules_id] = unpacked_feb[1::n_fields] + ## Unpack relative event ID + # event_container.feb_event_id[ + # self.camera_config.nectarcam.expected_modules_id] = unpacked_feb[2::n_fields] + ## Unpack TS1 counter + # event_container.feb_ts1[ + # self.camera_config.nectarcam.expected_modules_id] = unpacked_feb[3::n_fields] + ## Unpack TS2 counters + # ts2_decimal = lambda bits: bits - (1 << 8) if bits & 0x80 != 0 else bits + # ts2_decimal_vec = np.vectorize(ts2_decimal) + # event_container.feb_ts2_trig[ + # self.camera_config.nectarcam.expected_modules_id] = ts2_decimal_vec( + # unpacked_feb[4::n_fields]) + # event_container.feb_ts2_pps[ + # self.camera_config.nectarcam.expected_modules_id] = ts2_decimal_vec( + # unpacked_feb[5::n_fields]) + # Loop over modules + + for module_idx, module_id in enumerate( + self.camera_config.nectarcam.expected_modules_id + ): + offset = module_id * 7 + if bytes_per_module > 16: + field_id = 8 + # Decode trigger pattern + for pattern_id in range(n_patterns): + value = unpacked_feb[n_fields * module_idx + field_id + pattern_id] + module_pattern = [ + int(digit) for digit in reversed(bin(value)[2:].zfill(7)) + ] + event_container.trigger_pattern[ + pattern_id, offset : offset + 7 + ] = module_pattern + + # Unpack native charge + # if len(event.nectarcam.charges_gain1) > 0: + # event_container.native_charge = np.zeros(shape=(N_GAINS, N_PIXELS), + # dtype=np.uint16) + # rec_fmt = '=' + 'H' * self.camera_config.num_pixels + # for gain_id in range(N_GAINS): + # unpacked_charge = struct.unpack(rec_fmt, getattr(event.nectarcam, + # f'charges_gain{gain_id + 1}')) + # event_container.native_charge[ + # gain_id, self.camera_config.expected_pixels_id] = unpacked_charge + + +def fill_trigger_info(self, array_event): + tel_id = self.tel_id + nectarcam = array_event.nectarcam.tel[tel_id] + tib_available = nectarcam.evt.extdevices_presence & 1 + ucts_available = nectarcam.evt.extdevices_presence & 2 + # fill trigger time using UCTS timestamp + trigger = array_event.trigger + trigger_time = nectarcam.evt.ucts_timestamp + trigger_time = time_from_unix_tai_ns(trigger_time) + trigger.time = trigger_time + trigger.tels_with_trigger = [tel_id] + trigger.tel[tel_id].time = trigger.time + # decide which source to use, if both are available, + # the option decides, if not, fallback to the avilable source + # if no source available, warn and do not fill trigger info + if tib_available and ucts_available: + if self.default_trigger_type == "ucts": + trigger_bits = nectarcam.evt.ucts_trigger_type + else: + trigger_bits = nectarcam.evt.tib_masked_trigger + elif tib_available: + trigger_bits = nectarcam.evt.tib_masked_trigger + elif ucts_available: + trigger_bits = nectarcam.evt.ucts_trigger_type + else: + self.log.warning("No trigger info available.") + trigger.event_type = EventType.UNKNOWN + return + if ( + ucts_available + and nectarcam.evt.ucts_trigger_type == 42 # TODO check if it's correct + and self.default_trigger_type == "ucts" + ): + self.log.warning( + "Event with UCTS trigger_type 42 found." + " Probably means unreliable or shifted UCTS data." + ' Consider switching to TIB using `default_trigger_type="tib"`' + ) + # first bit mono trigger, second stereo. + # If *only* those two are set, we assume it's a physics event + # for all other we only check if the flag is present + if (trigger_bits & TriggerBits.PHYSICS) and not (trigger_bits & TriggerBits.OTHER): + trigger.event_type = EventType.SUBARRAY + elif trigger_bits & TriggerBits.CALIBRATION: + trigger.event_type = EventType.FLATFIELD + elif trigger_bits & TriggerBits.PEDESTAL: + trigger.event_type = EventType.SKY_PEDESTAL + elif trigger_bits & TriggerBits.SINGLE_PE: + trigger.event_type = EventType.SINGLE_PE + else: + self.log.warning( + f"Event {array_event.index.event_id} has unknown event type, trigger: {trigger_bits:08b}" + ) + trigger.event_type = EventType.UNKNOWN + + +class LightNectarCAMEventSource(NectarCAMEventSource): + def _generator(self): + # container for NectarCAM data + array_event = NectarCAMDataContainer() + array_event.meta["input_url"] = self.input_url + array_event.meta["max_events"] = self.max_events + array_event.meta["origin"] = "NectarCAM" + + # also add service container to the event section + array_event.nectarcam.tel[self.tel_id].svc = self.nectarcam_service + + # initialize general monitoring container + self.initialize_mon_container(array_event) + + # loop on events + for count, event in enumerate(self.multi_file): + array_event.count = count + array_event.index.event_id = event.event_id + array_event.index.obs_id = self.obs_ids[0] + + # fill R0/R1 data + self.fill_r0r1_container(array_event, event) + # fill specific NectarCAM event data + # fill specific NectarCAM event data + self.fill_nectarcam_event_container_from_zfile(array_event, event) + + if self.trigger_information: + self.fill_trigger_info(array_event) + + # fill general monitoring data + # self.fill_mon_container_from_zfile(array_event, event) + # + ## gain select and calibrate to pe + # if self.r0_r1_calibrator.calibration_path is not None: + # # skip flatfield and pedestal events if asked + # if ( + # self.calibrate_flatfields_and_pedestals + # or array_event.trigger.event_type not in {EventType.FLATFIELD, + # EventType.SKY_PEDESTAL} + # ): + # self.r0_r1_calibrator.calibrate(array_event) + + yield array_event diff --git a/src/nectarchain/data/container/gainContainer.py b/src/nectarchain/data/container/gainContainer.py new file mode 100644 index 00000000..7f2fb36f --- /dev/null +++ b/src/nectarchain/data/container/gainContainer.py @@ -0,0 +1,45 @@ +import logging + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + +import numpy as np +from ctapipe.containers import Field + +from .core import NectarCAMContainer + +__all__ = ["GainContainer", "SPEfitContainer"] + + +class GainContainer(NectarCAMContainer): + is_valid = Field(type=np.ndarray, dtype=bool, ndim=1, description="is_valid") + high_gain = Field( + type=np.ndarray, dtype=np.float64, ndim=2, description="high gain" + ) + low_gain = Field(type=np.ndarray, dtype=np.float64, ndim=2, description="low gain") + pixels_id = Field(type=np.ndarray, dtype=np.uint16, ndim=1, description="pixel ids") + + @classmethod + def from_hdf5(cls, path): + return super(__class__, cls)._container_from_hdf5(path, container_class=cls) + + +class SPEfitContainer(GainContainer): + likelihood = Field( + type=np.ndarray, dtype=np.float64, ndim=1, description="likelihood" + ) + p_value = Field(type=np.ndarray, dtype=np.float64, ndim=1, description="p_value") + pedestal = Field(type=np.ndarray, dtype=np.float64, ndim=2, description="pedestal") + pedestalWidth = Field( + type=np.ndarray, dtype=np.float64, ndim=2, description="pedestalWidth" + ) + resolution = Field( + type=np.ndarray, dtype=np.float64, ndim=2, description="resolution" + ) + luminosity = Field( + type=np.ndarray, dtype=np.float64, ndim=2, description="luminosity" + ) + mean = Field(type=np.ndarray, dtype=np.float64, ndim=2, description="mean") + n = Field(type=np.ndarray, dtype=np.float64, ndim=2, description="n") + pp = Field(type=np.ndarray, dtype=np.float64, ndim=2, description="pp") diff --git a/src/nectarchain/data/container/waveforms_container.py b/src/nectarchain/data/container/waveformsContainer.py similarity index 85% rename from src/nectarchain/data/container/waveforms_container.py rename to src/nectarchain/data/container/waveformsContainer.py index d479196e..45da1296 100644 --- a/src/nectarchain/data/container/waveforms_container.py +++ b/src/nectarchain/data/container/waveformsContainer.py @@ -1,14 +1,14 @@ import logging +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + import numpy as np from ctapipe.containers import Field, Map, partial from .core import ArrayDataContainer, TriggerMapContainer -logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") -log = logging.getLogger(__name__) -log.handlers = logging.getLogger("__main__").handlers - class WaveformsContainer(ArrayDataContainer): """ @@ -22,11 +22,10 @@ class WaveformsContainer(ArrayDataContainer): """ nsamples = Field( - type=int, + type=np.uint8, description="number of samples in the waveforms", ) - # subarray = Field(type=SubarrayDescription, - # description="The subarray description") + # subarray = Field(type=SubarrayDescription, description="The subarray description") wfs_hg = Field( type=np.ndarray, dtype=np.uint16, ndim=3, description="high gain waveforms" ) @@ -38,5 +37,5 @@ class WaveformsContainer(ArrayDataContainer): class WaveformsContainers(TriggerMapContainer): containers = Field( default_factory=partial(Map, WaveformsContainer), - description="trigger mapping of WaveformContainer", + description="trigger or slices of data mapping of WaveformContainer", ) diff --git a/src/nectarchain/data/management.py b/src/nectarchain/data/management.py index b524bc40..76736968 100644 --- a/src/nectarchain/data/management.py +++ b/src/nectarchain/data/management.py @@ -6,9 +6,12 @@ import glob import os +import pathlib from pathlib import Path from typing import List, Tuple +import numpy as np + __all__ = ["DataManagement"] @@ -162,3 +165,74 @@ def get_GRID_location( return lfns else: return url_data + + def find_waveforms(run_number, max_events=None): + return __class__.__find_computed_data( + run_number=run_number, max_events=max_events, data_type="waveforms" + ) + + def find_charges( + run_number, method="FullWaveformSum", str_extractor_kwargs="", max_events=None + ): + return __class__.__find_computed_data( + run_number=run_number, + max_events=max_events, + ext=f"_{method}_{str_extractor_kwargs}.h5", + data_type="charges", + ) + + def find_SPE_HHV(run_number, method="FullWaveformSum", str_extractor_kwargs=""): + full_file = glob.glob( + pathlib.Path( + f"{os.environ.get('NECTARCAMDATA','/tmp')}/SPEfit/FlatFieldSPEHHVStdNectarCAM_run{run_number}_{method}_{str_extractor_kwargs}.h5" + ).__str__() + ) + if len(full_file) != 1: + all_files = glob.glob( + pathlib.Path( + f"{os.environ.get('NECTARCAMDATA','/tmp')}/SPEfit/FlatFieldSPEHHVStdNectarCAM_run{run_number}_maxevents*_{method}_{str_extractor_kwargs}.h5" + ).__str__() + ) + max_events = 0 + for i, file in enumerate(all_files): + data = file.split("/")[-1].split(".h5")[0].split("_") + for _data in data: + if "maxevents" in _data: + _max_events = int(_data.split("maxevents")[-1]) + break + if _max_events >= max_events: + max_events = _max_events + index = i + return [all_files[index]] + else: + return full_file + + def __find_computed_data( + run_number, max_events=None, ext=".h5", data_type="waveforms" + ): + out = glob.glob( + pathlib.Path( + f"{os.environ.get('NECTARCAMDATA','/tmp')}/runs/{data_type}/*_run{run_number}{ext}" + ).__str__() + ) + if not (max_events is None): + all_files = glob.glob( + pathlib.Path( + f"{os.environ.get('NECTARCAMDATA','/tmp')}/runs/{data_type}/*_run{run_number}_maxevents*{ext}" + ).__str__() + ) + best_max_events = np.inf + best_index = None + for i, file in enumerate(all_files): + data = file.split("/")[-1].split(".h5")[0].split("_") + for _data in data: + if "maxevents" in _data: + _max_events = int(_data.split("maxevents")[-1]) + break + if _max_events >= max_events: + if _max_events < best_max_events: + best_max_events = _max_events + best_index = i + if not (best_index is None): + out = [all_files[best_index]] + return out diff --git a/src/nectarchain/dqm/start_dqm.py b/src/nectarchain/dqm/start_dqm.py index 8afe6e6f..b375ee59 100644 --- a/src/nectarchain/dqm/start_dqm.py +++ b/src/nectarchain/dqm/start_dqm.py @@ -51,7 +51,6 @@ if args.runnb is not None: # Grab runs automatically from DIRAC is the -r option is provided from nectarchain.data.container import utils - dm = utils.DataManagement() _, filelist = dm.findrun(args.runnb) args.input_files = [s.name for s in filelist] diff --git a/src/nectarchain/makers/__init__.py b/src/nectarchain/makers/__init__.py index 118008af..d690d0f6 100644 --- a/src/nectarchain/makers/__init__.py +++ b/src/nectarchain/makers/__init__.py @@ -1,3 +1,4 @@ -# from .charges_makers import * -# from .core import * -# from .waveforms_makers import * +# from .chargesMakers import * +from .chargesMakers import * +from .core import * +from .waveformsMakers import * diff --git a/src/nectarchain/makers/calibration/__init__.py b/src/nectarchain/makers/calibration/__init__.py index 81f96f56..0ded83f5 100644 --- a/src/nectarchain/makers/calibration/__init__.py +++ b/src/nectarchain/makers/calibration/__init__.py @@ -1,3 +1,3 @@ -# from .flatfieldMakers import * -# from .gain import * -# from .pedestalMakers import * +from .flatfieldMakers import * +from .gain import * +from .pedestalMakers import * diff --git a/src/nectarchain/makers/calibration/core.py b/src/nectarchain/makers/calibration/core.py index 7c4f7db1..c6c46e5f 100644 --- a/src/nectarchain/makers/calibration/core.py +++ b/src/nectarchain/makers/calibration/core.py @@ -1,78 +1,45 @@ import logging -from collections.abc import Iterable -from copy import copy -from datetime import date -from pathlib import Path - -import astropy.units as u -import numpy as np -from astropy.table import Column, QTable - -from ..core import BaseMaker logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") log = logging.getLogger(__name__) log.handlers = logging.getLogger("__main__").handlers -__all__ = ["CalibrationMaker"] +from ctapipe.core.traits import List +from ..core import EventsLoopNectarCAMCalibrationTool -class CalibrationMaker(BaseMaker): - """ - Mother class for all calibration makers that can be defined to compute - calibration coefficients from data. +__all__ = [""] - Attributes: - _reduced_name (str): A string representing the name of the calibration. - PIXELS_ID_COLUMN (str): A string representing the name of the column in the - result table that stores the pixels id. - NP_PIXELS (str): A string representing the key in the metadata that stores - the number of pixels. - Members: - _pixels_id (ndarray): A private property that stores the pixels id. - _results (QTable): A private property that stores the result table. - """ +class NectarCAMCalibrationTool(EventsLoopNectarCAMCalibrationTool): + name = "CalibrationTool" + # PIXELS_ID_COLUMN = "pixels_id" + # NP_PIXELS = "npixels" - _reduced_name = "Calibration" - PIXELS_ID_COLUMN = "pixels_id" - NP_PIXELS = "npixels" + pixels_id = List( + default_value=None, + help="the list of pixel id to apply the components", + allow_none=True, + ).tag(config=True) - def __new__(cls, *args, **kwargs): - """ - Constructor. - Returns: - CalibrationMaker: An instance of the CalibrationMaker class. - """ - return super(CalibrationMaker, cls).__new__(cls) - - def __init__(self, pixels_id, *args, **kwargs) -> None: - """ - Initialize the CalibrationMaker object. - - Args: - pixels_id (iterable, np.ndarray): The list of pixels id. - """ - super().__init__() - if not (isinstance(pixels_id, Iterable)): - raise TypeError("pixels_id must be iterable") - self.__pixels_id = np.array(pixels_id) +''' + def setup(self) -> None: + super().setup() self.__results = QTable() self.__results.add_column( Column( - self.__pixels_id, - __class__.PIXELS_ID_COLUMN, + self.pixels_id, + self.pixels_id.name, unit=u.dimensionless_unscaled, ) ) self.__results.meta[__class__.NP_PIXELS] = self.npixels - self.__results.meta["comments"] = ( - f"Produced with NectarChain, Credit : CTA NectarCam" - f' {date.today().strftime("%B %d, %Y")}' - ) + self.__results.meta[ + "comments" + ] = f'Produced with NectarChain, Credit : CTA NectarCam {date.today().strftime("%B %d, %Y")}' - def save(self, path, **kwargs): + def finish(self, path, **kwargs): """ Saves the results to a file in the specified path. @@ -92,45 +59,6 @@ def save(self, path, **kwargs): overwrite=kwargs.get("overwrite", False), ) - @property - def _pixels_id(self): - """ - Get the pixels id. - - Returns: - ndarray: The pixels id. - """ - return self.__pixels_id - - @_pixels_id.setter - def _pixels_id(self, value): - """ - Set the pixels id. - - Args: - value (ndarray): The pixels id. - """ - self.__pixels_id = value - - @property - def pixels_id(self): - """ - Get a copy of the pixels id. - - Returns: - ndarray: A copy of the pixels id. - """ - return copy(self.__pixels_id) - - @property - def npixels(self): - """ - Get the number of pixels. - - Returns: - int: The number of pixels. - """ - return len(self.__pixels_id) @property def _results(self): @@ -151,3 +79,4 @@ def results(self): QTable: A copy of the result table. """ return copy(self.__results) +''' diff --git a/src/nectarchain/makers/calibration/flatfield_makers.py b/src/nectarchain/makers/calibration/flatfieldMakers.py similarity index 51% rename from src/nectarchain/makers/calibration/flatfield_makers.py rename to src/nectarchain/makers/calibration/flatfieldMakers.py index e38ee805..01da9958 100644 --- a/src/nectarchain/makers/calibration/flatfield_makers.py +++ b/src/nectarchain/makers/calibration/flatfieldMakers.py @@ -1,20 +1,16 @@ import logging -from .core import CalibrationMaker - logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") log = logging.getLogger(__name__) log.handlers = logging.getLogger("__main__").handlers -__all__ = ["FlatfieldMaker"] +from .core import NectarCAMCalibrationTool +__all__ = ["FlatfieldNectarCAMCalibrationTool"] -class FlatfieldMaker(CalibrationMaker): - def __init__(self, *args, **kwargs) -> None: - super().__init__(*args, **kwargs) - def make(self): +class FlatfieldNectarCAMCalibrationTool(NectarCAMCalibrationTool): + def start(self): raise NotImplementedError( - "The computation of the flatfield calibration is not yet implemented, " - "feel free to contribute !:)" + "The computation of the flatfield calibration is not yet implemented, feel free to contribute !:)" ) diff --git a/src/nectarchain/makers/calibration/gain/FlatFieldSPEMakers.py b/src/nectarchain/makers/calibration/gain/FlatFieldSPEMakers.py new file mode 100644 index 00000000..8f6948e9 --- /dev/null +++ b/src/nectarchain/makers/calibration/gain/FlatFieldSPEMakers.py @@ -0,0 +1,197 @@ +import logging + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + +import os +import pathlib + +import numpy as np +from ctapipe.containers import Container +from ctapipe.core.traits import ComponentNameList + +from ....data.container import ChargesContainer, ChargesContainers +from ....data.container.core import merge_map_ArrayDataContainer +from ....data.management import DataManagement +from ...component import ArrayDataComponent, NectarCAMComponent +from ...extractor.utils import CtapipeExtractor +from .core import GainNectarCAMCalibrationTool + +__all__ = [ + "FlatFieldSPENominalNectarCAMCalibrationTool", + "FlatFieldSPENominalStdNectarCAMCalibrationTool", + "FlatFieldSPEHHVNectarCAMCalibrationTool", + "FlatFieldSPEHHVStdNectarCAMCalibrationTool", + "FlatFieldSPECombinedStdNectarCAMCalibrationTool", +] + + +class FlatFieldSPENominalNectarCAMCalibrationTool(GainNectarCAMCalibrationTool): + name = "FlatFieldSPEHHVNectarCAM" + componentsList = ComponentNameList( + NectarCAMComponent, + default_value=["FlatFieldSingleHHVSPENectarCAMComponent"], + help="List of Component names to be apply, the order will be respected", + ).tag(config=True) + + # events_per_slice = Integer( + # help="feature desactivated for this class", + # default_value=None, + # allow_none=True, + # read_only = True, + # ).tag(config=True) + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + str_extractor_kwargs = CtapipeExtractor.get_extractor_kwargs_str( + self.extractor_kwargs + ) + if not (self.reload_events): + files = DataManagement.find_charges( + run_number=self.run_number, + method=self.method, + str_extractor_kwargs=str_extractor_kwargs, + max_events=self.max_events, + ) + if len(files) == 1: + log.warning( + "You asked events_per_slice but you don't want to reload events and a charges file is on disk, then events_per_slice is set to None" + ) + self.events_per_slice = None + + def _init_output_path(self): + str_extractor_kwargs = CtapipeExtractor.get_extractor_kwargs_str( + self.extractor_kwargs + ) + if self.events_per_slice is None: + ext = ".h5" + else: + ext = f"_sliced{self.events_per_slice}.h5" + if self.max_events is None: + filename = f"{self.name}_run{self.run_number}_{self.method}_{str_extractor_kwargs}{ext}" + else: + filename = f"{self.name}_run{self.run_number}_maxevents{self.max_events}_{self.method}_{str_extractor_kwargs}{ext}" + + self.output_path = pathlib.Path( + f"{os.environ.get('NECTARCAMDATA','/tmp')}/SPEfit/{filename}" + ) + + def start( + self, + n_events=np.inf, + # trigger_type: list = None, + restart_from_begining: bool = False, + *args, + **kwargs, + ): + str_extractor_kwargs = CtapipeExtractor.get_extractor_kwargs_str( + self.extractor_kwargs + ) + files = DataManagement.find_charges( + run_number=self.run_number, + method=self.method, + str_extractor_kwargs=str_extractor_kwargs, + max_events=self.max_events, + ) + if self.reload_events or len(files) != 1: + if len(files) != 1: + self.log.info( + f"{len(files)} computed charges files found with max_events > {self.max_events} for run {self.run_number} with extraction method {self.method} and {str_extractor_kwargs},\n reload charges from event loop" + ) + super().start( + n_events=n_events, + restart_from_begining=restart_from_begining, + *args, + **kwargs, + ) + else: + self.log.info(f"reading computed charge from files {files[0]}") + chargesContainers = ChargesContainer.from_hdf5(files[0]) + if isinstance(chargesContainers, ChargesContainer): + self.components[0]._chargesContainers = chargesContainers + elif isinstance(chargesContainers, ChargesContainers): + self.log.debug("merging along TriggerType") + self.components[0]._chargesContainers = merge_map_ArrayDataContainer( + chargesContainers + ) + else: + self.log.debug("merging along slices") + chargesContaienrs_merdes_along_slices = ( + ArrayDataComponent.merge_along_slices( + containers_generator=chargesContainers + ) + ) + self.log.debug("merging along TriggerType") + self.components[0]._chargesContainers = merge_map_ArrayDataContainer( + chargesContaienrs_merdes_along_slices + ) + + def _write_container(self, container: Container, index_component: int = 0) -> None: + # if isinstance(container,SPEfitContainer) : + # self.writer.write(table_name = f"{self.method}_{CtapipeExtractor.get_extractor_kwargs_str(self.extractor_kwargs)}", + # containers = container, + # ) + # else : + super()._write_container(container=container, index_component=index_component) + + +class FlatFieldSPEHHVNectarCAMCalibrationTool( + FlatFieldSPENominalNectarCAMCalibrationTool +): + name = "FlatFieldSPEHHVNectarCAM" + componentsList = ComponentNameList( + NectarCAMComponent, + default_value=["FlatFieldSingleHHVSPENectarCAMComponent"], + help="List of Component names to be apply, the order will be respected", + ).tag(config=True) + + +class FlatFieldSPEHHVStdNectarCAMCalibrationTool( + FlatFieldSPENominalNectarCAMCalibrationTool +): + name = "FlatFieldSPEHHVStdNectarCAM" + componentsList = ComponentNameList( + NectarCAMComponent, + default_value=["FlatFieldSingleHHVSPEStdNectarCAMComponent"], + help="List of Component names to be apply, the order will be respected", + ).tag(config=True) + + +class FlatFieldSPENominalStdNectarCAMCalibrationTool( + FlatFieldSPENominalNectarCAMCalibrationTool +): + name = "FlatFieldSPENominalStdNectarCAM" + componentsList = ComponentNameList( + NectarCAMComponent, + default_value=["FlatFieldSingleNominalSPEStdNectarCAMComponent"], + help="List of Component names to be apply, the order will be respected", + ).tag(config=True) + + +class FlatFieldSPECombinedStdNectarCAMCalibrationTool( + FlatFieldSPENominalNectarCAMCalibrationTool +): + name = "FlatFieldCombinedStddNectarCAM" + componentsList = ComponentNameList( + NectarCAMComponent, + default_value=["FlatFieldCombinedSPEStdNectarCAMComponent"], + help="List of Component names to be apply, the order will be respected", + ).tag(config=True) + + def _init_output_path(self): + for word in self.SPE_result.stem.split("_"): + if "run" in word: + HHVrun = int(word.split("run")[-1]) + str_extractor_kwargs = CtapipeExtractor.get_extractor_kwargs_str( + self.extractor_kwargs + ) + if self.max_events is None: + filename = f"{self.name}_run{self.run_number}_HHV{HHVrun}_{self.method}_{str_extractor_kwargs}.h5" + else: + filename = f"{self.name}_run{self.run_number}_maxevents{self.max_events}_HHV{HHVrun}_{self.method}_{str_extractor_kwargs}.h5" + + self.output_path = pathlib.Path( + f"{os.environ.get('NECTARCAMDATA','/tmp')}/SPEfit/{filename}" + ) diff --git a/src/nectarchain/makers/calibration/gain/whitetarget_spe_makers.py b/src/nectarchain/makers/calibration/gain/WhiteTargetSPEMakers.py similarity index 76% rename from src/nectarchain/makers/calibration/gain/whitetarget_spe_makers.py rename to src/nectarchain/makers/calibration/gain/WhiteTargetSPEMakers.py index 139c68be..be726219 100644 --- a/src/nectarchain/makers/calibration/gain/whitetarget_spe_makers.py +++ b/src/nectarchain/makers/calibration/gain/WhiteTargetSPEMakers.py @@ -1,12 +1,12 @@ import logging -from .gain_makers import GainMaker - logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") log = logging.getLogger(__name__) log.handlers = logging.getLogger("__main__").handlers -__all__ = ["WhiteTargetSPEMaker"] +from .gainMakers import GainMaker + +__all__ = ["FlatfieldMaker"] class WhiteTargetSPEMaker(GainMaker): @@ -15,6 +15,5 @@ def __init__(self, *args, **kwargs) -> None: def make(self): raise NotImplementedError( - "The computation of the white target calibration is not yet implemented, " - "feel free to contribute !:)" + "The computation of the white target calibration is not yet implemented, feel free to contribute !:)" ) diff --git a/src/nectarchain/makers/calibration/gain/__init__.py b/src/nectarchain/makers/calibration/gain/__init__.py index f8a87d7d..518e6584 100644 --- a/src/nectarchain/makers/calibration/gain/__init__.py +++ b/src/nectarchain/makers/calibration/gain/__init__.py @@ -1,3 +1,5 @@ -# from .flatfield_spe_makers import * +from .FlatFieldSPEMakers import * +from .photostat_makers import * + # from .WhiteTargetSPEMakers import * -# from .photostatistic_makers import * +# from .PhotoStatisticMakers import * diff --git a/src/nectarchain/makers/calibration/gain/core.py b/src/nectarchain/makers/calibration/gain/core.py new file mode 100644 index 00000000..095fc32c --- /dev/null +++ b/src/nectarchain/makers/calibration/gain/core.py @@ -0,0 +1,17 @@ +import logging + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + +from ctapipe.core.traits import Bool + +from ..core import NectarCAMCalibrationTool + +__all__ = ["GainNectarCAMCalibrationTool"] + + +class GainNectarCAMCalibrationTool(NectarCAMCalibrationTool): + reload_events = Bool( + default_value=False, help="a flag to re compute the charge from raw data" + ) diff --git a/src/nectarchain/makers/calibration/gain/gain_makers.py b/src/nectarchain/makers/calibration/gain/gain_makers.py deleted file mode 100644 index 6743ff11..00000000 --- a/src/nectarchain/makers/calibration/gain/gain_makers.py +++ /dev/null @@ -1,138 +0,0 @@ -import logging -from copy import copy - -import astropy.units as u -import numpy as np -from astropy.table import Column - -from ..core import CalibrationMaker - -logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") -log = logging.getLogger(__name__) -log.handlers = logging.getLogger("__main__").handlers - -__all__ = ["GainMaker"] - - -class GainMaker(CalibrationMaker): - """ - A class for gain calibration calculations on data. - - Inherits from the `CalibrationMaker` class and adds functionality specific to - gain calibration. - - Members: - __high_gain (ndarray): Private field to store the high gain values. - __low_gain (ndarray): Private field to store the low gain values. - - Methods: - __init__(self, *args, **kwargs): Initializes the `GainMaker` object and sets - up the result table with columns for high gain, high gain error, low gain, - low gain error, and validity flag. - _high_gain.setter: Sets the high gain values. - high_gain(self): Returns a copy of the high gain values. - _low_gain.setter: Sets the low gain values. - low_gain(self): Returns a copy of the low gain values. - """ - - def __init__(self, *args, **kwargs): - """ - Initializes the `GainMaker` object and sets up the result table with columns - for high gain, high gain error, low gain, low gain error, and validity flag. - - Args: - *args: Variable length argument list. - **kwargs: Arbitrary keyword arguments. - """ - super().__init__(*args, **kwargs) - self.__high_gain = np.empty((self.npixels), dtype=np.float64) - self.__low_gain = np.empty((self.npixels), dtype=np.float64) - self._results.add_column( - Column( - data=self.__high_gain, name="high_gain", unit=u.dimensionless_unscaled - ) - ) - self._results.add_column( - Column( - np.empty((self.npixels, 2), dtype=np.float64), - "high_gain_error", - unit=u.dimensionless_unscaled, - ) - ) - self._results.add_column( - Column(data=self.__low_gain, name="low_gain", unit=u.dimensionless_unscaled) - ) - self._results.add_column( - Column( - np.empty((self.npixels, 2), dtype=np.float64), - "low_gain_error", - unit=u.dimensionless_unscaled, - ) - ) - self._results.add_column( - Column( - np.zeros((self.npixels), dtype=bool), - "is_valid", - unit=u.dimensionless_unscaled, - ) - ) - - @property - def _high_gain(self): - """ - Getter for the high gain values. - - Returns: - ndarray: A copy of the high gain values. - """ - return self.__high_gain - - @_high_gain.setter - def _high_gain(self, value): - """ - Setter for the high gain values. - - Args: - value (ndarray): The high gain values. - """ - self.__high_gain = value - - @property - def high_gain(self): - """ - Getter for the high gain values. - - Returns: - ndarray: A copy of the high gain values. - """ - return copy(self.__high_gain) - - @property - def _low_gain(self): - """ - Getter for the low gain values. - - Returns: - ndarray: A copy of the low gain values. - """ - return self.__low_gain - - @_low_gain.setter - def _low_gain(self, value): - """ - Setter for the low gain values. - - Args: - value (ndarray): The low gain values. - """ - self.__low_gain = value - - @property - def low_gain(self): - """ - Getter for the low gain values. - - Returns: - ndarray: A copy of the low gain values. - """ - return copy(self.__low_gain) diff --git a/src/nectarchain/makers/calibration/gain/photostat_makers.py b/src/nectarchain/makers/calibration/gain/photostat_makers.py new file mode 100644 index 00000000..3249170c --- /dev/null +++ b/src/nectarchain/makers/calibration/gain/photostat_makers.py @@ -0,0 +1,186 @@ +import logging + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + +import os +import pathlib + +import numpy as np +from ctapipe.containers import Container, EventType +from ctapipe.core.traits import ComponentNameList, Integer, Path + +from ....data.container import ChargesContainer +from ....data.container.core import NectarCAMContainer, merge_map_ArrayDataContainer +from ....data.management import DataManagement +from ...component import ArrayDataComponent, NectarCAMComponent +from ...extractor.utils import CtapipeExtractor +from .core import GainNectarCAMCalibrationTool + +__all__ = ["PhotoStatisticNectarCAMCalibrationTool"] + + +class PhotoStatisticNectarCAMCalibrationTool(GainNectarCAMCalibrationTool): + ###TO DO : IMPLEMENT a MOTHER PHOTOSTAT CLASS WITH ONLY 1 RUN WITH FF AND PEDESTAL INTERLEAVED. + + name = "PhotoStatisticNectarCAM" + componentsList = ComponentNameList( + NectarCAMComponent, + default_value=["PhotoStatisticNectarCAMComponent"], + help="List of Component names to be apply, the order will be respected", + ).tag(config=True) + run_number = Integer(help="the FF run number to be treated", default_value=-1).tag( + config=True + ) + Ped_run_number = Integer( + help="the FF run number to be treated", default_value=-1 + ).tag(config=True) + + run_file = Path( + help="desactivated for PhotoStatistic maker with FF and pedestal runs separated", + default_value=None, + allow_none=True, + read_only=True, + ).tag(config=False) + + events_per_slice = Integer( + help="desactivated for PhotoStatistic maker with FF and pedestal runs separated", + default_value=None, + allow_none=True, + read_only=True, + ).tag(config=False) + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + str_extractor_kwargs = CtapipeExtractor.get_extractor_kwargs_str( + self.extractor_kwargs + ) + if not (self.reload_events): + FF_files = DataManagement.find_charges( + run_number=self.run_number, + method=self.method, + str_extractor_kwargs=str_extractor_kwargs, + max_events=self.max_events, + ) + Ped_files = DataManagement.find_charges( + run_number=self.run_number, + max_events=self.max_events, + ) + + if len(FF_files) == 1 and len(Ped_files) == 1: + log.warning( + "You asked events_per_slice but you don't want to reload events and a charges file is on disk for FF and Pedestals, then events_per_slice is set to None" + ) + self.events_per_slice = None + + def _init_output_path(self): + str_extractor_kwargs = CtapipeExtractor.get_extractor_kwargs_str( + self.extractor_kwargs + ) + if self.max_events is None: + filename = f"{self.name}_FFrun{self.run_number}_{self.method}_{str_extractor_kwargs}_Pedrun{self.Ped_run_number}_FullWaveformSum.h5" + else: + filename = f"{self.name}_FFrun{self.run_number}_{self.method}_{str_extractor_kwargs}_Pedrun{self.Ped_run_number}_FullWaveformSum_maxevents{self.max_events}.h5" + self.output_path = pathlib.Path( + f"{os.environ.get('NECTARCAMDATA','/tmp')}/PhotoStat/{filename}" + ) + + def _load_eventsource(self, FF_run=True): + if FF_run: + self.log.debug("loading FF event source") + self.event_source = self.enter_context( + self.load_run(self.run_number, self.max_events) + ) + else: + self.log.debug("loading Ped event source") + self.event_source = self.enter_context( + self.load_run(self.Ped_run_number, self.max_events) + ) + + def start( + self, + n_events=np.inf, + *args, + **kwargs, + ): + str_extractor_kwargs = CtapipeExtractor.get_extractor_kwargs_str( + self.extractor_kwargs + ) + FF_files = DataManagement.find_charges( + run_number=self.run_number, + method=self.method, + str_extractor_kwargs=str_extractor_kwargs, + max_events=self.max_events, + ) + Ped_files = DataManagement.find_charges( + run_number=self.run_number, + max_events=self.max_events, + ) + if self.reload_events or len(FF_files) != 1 or len(Ped_files) != 1: + if len(FF_files) != 1 or len(Ped_files) != 1: + self.log.info( + f"{len(FF_files)} computed charges FF files found with max_events > {self.max_events} for run {self.run_number} with extraction method {self.method} and {str_extractor_kwargs},\n reload charges from event loop" + ) + self.log.info( + f"{len(Ped_files)} computed charges FF files found with max_events > {self.max_events} for run {self.Ped_run_number} with extraction method FullWaveformSum,\n reload charges from event loop" + ) + + super().start( + n_events=n_events, restart_from_begining=False, *args, **kwargs + ) + self._setup_eventsource(FF_run=False) + super().start( + n_events=n_events, restart_from_begining=False, *args, **kwargs + ) + else: + self.log.info(f"reading computed charge from FF file {FF_files[0]}") + chargesContainers = ChargesContainer.from_hdf5(FF_files[0]) + if isinstance(chargesContainers, NectarCAMContainer): + self.components[0]._FF_chargesContainers = chargesContainers + elif isinstance(list(chargesContainers.containers.keys())[0], EventType): + self.log.debug("merging along TriggerType") + self.components[0]._FF_chargesContainers = merge_map_ArrayDataContainer( + chargesContainers + ) + else: + self.log.debug("merging along slices") + chargesContaienrs_merdes_along_slices = ( + ArrayDataComponent.merge_along_slices(chargesContainers) + ) + self.log.debug("merging along TriggerType") + self.components[0]._FF_chargesContainers = merge_map_ArrayDataContainer( + chargesContaienrs_merdes_along_slices + ) + + self.log.info(f"reading computed charge from Ped file {Ped_files[0]}") + chargesContainers = ChargesContainer.from_hdf5(Ped_files[0]) + if isinstance(chargesContainers, NectarCAMContainer): + self.components[0]._Ped_chargesContainers = chargesContainers + elif isinstance(list(chargesContainers.containers.keys())[0], EventType): + self.log.debug("merging along TriggerType") + self.components[ + 0 + ]._Ped_chargesContainers = merge_map_ArrayDataContainer( + chargesContainers + ) + else: + self.log.debug("merging along slices") + chargesContaienrs_merdes_along_slices = ( + ArrayDataComponent.merge_along_slices(chargesContainers) + ) + self.log.debug("merging along TriggerType") + self.components[ + 0 + ]._Ped_chargesContainers = merge_map_ArrayDataContainer( + chargesContaienrs_merdes_along_slices + ) + + def _write_container(self, container: Container, index_component: int = 0) -> None: + # if isinstance(container,SPEfitContainer) : + # self.writer.write(table_name = f"{self.method}_{CtapipeExtractor.get_extractor_kwargs_str(self.extractor_kwargs)}", + # containers = container, + # ) + # else : + super()._write_container(container=container, index_component=index_component) diff --git a/src/nectarchain/makers/calibration/gain/photostatistic_makers.py b/src/nectarchain/makers/calibration/gain/photostatistic_makers.py deleted file mode 100644 index e472401c..00000000 --- a/src/nectarchain/makers/calibration/gain/photostatistic_makers.py +++ /dev/null @@ -1,620 +0,0 @@ -import copy -import logging -import os -from pathlib import Path - -import numpy as np -from astropy.table import QTable -from astropy.visualization import quantity_support -from ctapipe_io_nectarcam import constants -from matplotlib import pyplot as plt -from scipy.stats import linregress - -from ....data.container import ChargesContainer, ChargesContainerIO -from ...charges_makers import ChargesMaker -from .gain_makers import GainMaker - -logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") -log = logging.getLogger(__name__) -log.handlers = logging.getLogger("__main__").handlers - -__all__ = ["PhotoStatisticMaker"] - - -class PhotoStatisticMaker(GainMaker): - """ - The `PhotoStatisticMaker` class is a subclass of `GainMaker` and is used to - calculate photo statistics for a given set of charge data. It provides methods to - create an instance from charge containers or run numbers, as well as methods to - calculate various statistics such as gain and standard deviation. - - Example Usage: - # Create an instance of PhotoStatisticMaker using charge containers - FFcharge = ChargesContainer(...) - Pedcharge = ChargesContainer(...) - coefCharge_FF_Ped = 0.5 - SPE_result = "path/to/SPE_results" - photo_stat = PhotoStatisticMaker.create_from_chargeContainer(FFcharge, - Pedcharge, coefCharge_FF_Ped, SPE_result) - - # Calculate and retrieve the gain values - gain_hg = photo_stat.gainHG - gain_lg = photo_stat.gainLG - - # Plot the correlation between photo statistic gain and SPE gain - photo_stat_gain = np.array(...) - SPE_gain = np.array(...) - fig = PhotoStatisticMaker.plot_correlation(photo_stat_gain, SPE_gain) - - Methods: - - `__init__(self, FFcharge_hg, FFcharge_lg, Pedcharge_hg, Pedcharge_lg, - coefCharge_FF_Ped, SPE_resolution, *args, **kwargs)`: Constructor method to - initialize the `PhotoStatisticMaker` instance with charge data and other - parameters. - - `create_from_chargeContainer(cls, FFcharge, Pedcharge, coefCharge_FF_Ped, - SPE_result, **kwargs)`: Class method to create an instance of - `PhotoStatisticMaker` from charge containers. - - `create_from_run_numbers(cls, FFrun, Pedrun, SPE_result, **kwargs)`: Class - method to create an instance of `PhotoStatisticMaker` from run numbers. - - `__readSPE(SPEresults) -> tuple`: Static method to read SPE resolution from a - file and return the resolution and pixel IDs. - - `__get_charges_FF_Ped_reshaped(FFcharge, Pedcharge, SPE_resolution, - SPE_pixels_id) -> dict`: Static method to reshape the charge data based on - the intersection of pixel IDs and return a dictionary of reshaped data. - - `__readFF(FFRun, **kwargs) -> dict`: Static method to read FF data from a - file and return the FF charge data and coefficient. - - `__readPed(PedRun, **kwargs) -> dict`: Static method to read Ped data from - a file and return the Ped charge data. - - `__check_shape(self) -> None`: Method to check the shape of the - charge data arrays. - - `make(self, **kwargs) -> None`: Method to run the photo statistic method - and store the results. - - `plot_correlation(photoStat_gain, SPE_gain) -> fig`: Static method to plot - the correlation between photo statistic gain and SPE gain. - - Fields: - - `SPE_resolution`: Property to get the SPE resolution. - - `sigmaPedHG`: Property to get the standard deviation of Pedcharge_hg. - - `sigmaChargeHG`: Property to get the standard deviation of FFcharge_hg - - meanPedHG. - - `meanPedHG`: Property to get the mean of Pedcharge_hg. - - `meanChargeHG`: Property to get the mean of FFcharge_hg - meanPedHG. - - `BHG`: Property to calculate the BHG value. - - `gainHG`: Property to calculate the gain for high gain. - - `sigmaPedLG`: Property to get the standard deviation of Pedcharge_lg. - - `sigmaChargeLG`: Property to get the standard deviation of FFcharge_lg - - meanPedLG. - - `meanPedLG`: Property to get the mean of Pedcharge_lg. - - `meanChargeLG`: Property to get the mean of FFcharge_lg - meanPedLG. - - `BLG`: Property to calculate the BLG value. - - `gainLG`: Property to calculate the gain for low gain. - """ - - _reduced_name = "PhotoStatistic" - - # constructors - def __init__( - self, - FFcharge_hg: np.ndarray, - FFcharge_lg: np.ndarray, - Pedcharge_hg: np.ndarray, - Pedcharge_lg: np.ndarray, - coefCharge_FF_Ped: float, - SPE_resolution, - *args, - **kwargs, - ) -> None: - """ - Initializes the instance of the PhotoStatisticMaker class with charge data - and other parameters. - - Args: - FFcharge_hg (np.ndarray): Array of charge data for high gain in the FF ( - Flat Field) image. - FFcharge_lg (np.ndarray): Array of charge data for low gain in the FF image. - Pedcharge_hg (np.ndarray): Array of charge data for high gain in the Ped - (Pedestal) image. - Pedcharge_lg (np.ndarray): Array of charge data for low gain in the Ped - image. - coefCharge_FF_Ped (float): Coefficient to convert FF charge to Ped charge. - SPE_resolution: Array-like of single photoelectron (SPE) resolutions for - each pixel, or single value to use the same for each pixel. - - Raises: - TypeError: If SPE_resolution is not provided in a valid format. - - Returns: - None - """ - super().__init__(*args, **kwargs) - - self.__coefCharge_FF_Ped = coefCharge_FF_Ped - - self.__FFcharge_hg = FFcharge_hg - self.__FFcharge_lg = FFcharge_lg - - self.__Pedcharge_hg = Pedcharge_hg - self.__Pedcharge_lg = Pedcharge_lg - - if ( - isinstance(SPE_resolution, np.ndarray) - and len(SPE_resolution) == self.npixels - ): - self.__SPE_resolution = SPE_resolution - elif isinstance(SPE_resolution, list) and len(SPE_resolution) == self.npixels: - self.__SPE_resolution = np.array(SPE_resolution) - elif isinstance(SPE_resolution, float): - self.__SPE_resolution = SPE_resolution * np.ones((self.npixels)) - else: - e = TypeError( - "SPE_resolution must be a float, a numpy.ndarray or list instance" - ) - raise e - - self.__check_shape() - - @classmethod - def create_from_chargeContainer( - cls, - FFcharge: ChargesContainer, - Pedcharge: ChargesContainer, - coefCharge_FF_Ped: float, - SPE_result, - **kwargs, - ): - """ - Create an instance of the PhotoStatisticMaker class from Pedestal and - Flatfield runs stored in ChargesContainer. - - Args: - FFcharge (ChargesContainer): Array of charge data for the FF image. - Pedcharge (ChargesContainer): Array of charge data for the Ped image. - coefCharge_FF_Ped (float): Coefficient to convert FF charge to Ped charge. - SPE_result (str or Path): Path to the SPE result file (optional). - **kwargs: Additional keyword arguments for initializing the - PhotoStatisticMaker instance. - - Returns: - PhotoStatisticMaker: An instance of the PhotoStatisticMaker class created - from the ChargesContainer instances. - """ - if isinstance(SPE_result, str) or isinstance(SPE_result, Path): - SPE_resolution, SPE_pixels_id = __class__.__readSPE(SPE_result) - else: - SPE_pixels_id = None - - kwargs_init = __class__.__get_charges_FF_Ped_reshaped( - FFcharge, Pedcharge, SPE_resolution, SPE_pixels_id - ) - - kwargs.update(kwargs_init) - return cls(coefCharge_FF_Ped=coefCharge_FF_Ped, **kwargs) - - @classmethod - def create_from_run_numbers( - cls, FFrun: int, Pedrun: int, SPE_result: str, **kwargs - ): - """ - Create an instance of the PhotoStatisticMaker class by reading the FF (Flat - Field) and Ped (Pedestal) charge data from run numbers. - - Args: - FFrun (int): The run number for the FF charge data. - Pedrun (int): The run number for the Ped charge data. - SPE_result (str): The path to the SPE result file. - **kwargs: Additional keyword arguments. - - Returns: - PhotoStatisticMaker: An instance of the PhotoStatisticMaker class created - from the FF and Ped charge data and the SPE result file. - """ - FFkwargs = __class__.__readFF(FFrun, **kwargs) - Pedkwargs = __class__.__readPed(Pedrun, **kwargs) - kwargs.update(FFkwargs) - kwargs.update(Pedkwargs) - return cls.create_from_chargeContainer(SPE_result=SPE_result, **kwargs) - - # methods - @staticmethod - def __readSPE(SPEresults) -> tuple: - """ - Reads the SPE resolution from a file and returns the resolution values and - corresponding pixel IDs. - - Args: - SPEresults (str): The file path to the SPE results file. - - Returns: - tuple: A tuple containing the SPE resolution values and corresponding - pixel IDs. - """ - log.info(f"reading SPE resolution from {SPEresults}") - table = QTable.read(SPEresults) - table.sort("pixels_id") - return ( - table["resolution"][table["is_valid"]].value, - table["pixels_id"][table["is_valid"]].value, - ) - - @staticmethod - def __get_charges_FF_Ped_reshaped( - FFcharge: ChargesContainer, - Pedcharge: ChargesContainer, - SPE_resolution: np.ndarray, - SPE_pixels_id: np.ndarray, - ) -> dict: - """ - Reshapes the FF (Flat Field) and Ped (Pedestal) charges based on the - intersection of pixel IDs between the two charges. - Selects the charges for the high-gain and low-gain channels and returns them - along with the common pixel IDs. - - Args: - FFcharge (ChargesContainer): The charges container for the Flat Field data. - Pedcharge (ChargesContainer): The charges container for the Pedestal data. - SPE_resolution (np.ndarray): An array containing the SPE resolutions. - SPE_pixels_id (np.ndarray): An array containing the pixel IDs for the SPE - data. - - Returns: - dict: A dictionary containing the reshaped data, including the common - pixel IDs, SPE resolution (if provided), and selected charges for the - high-gain and low-gain channels. - """ - log.info("reshape of SPE, Ped and FF data with intersection of pixel ids") - out = {} - - FFped_intersection = np.intersect1d(Pedcharge.pixels_id, FFcharge.pixels_id) - if not (SPE_pixels_id is None): - SPEFFPed_intersection = np.intersect1d(FFped_intersection, SPE_pixels_id) - mask_SPE = np.array( - [ - SPE_pixels_id[i] in SPEFFPed_intersection - for i in range(len(SPE_pixels_id)) - ], - dtype=bool, - ) - out["SPE_resolution"] = SPE_resolution[mask_SPE] - - out["pixels_id"] = SPEFFPed_intersection - out["FFcharge_hg"] = ChargesMaker.select_charges_hg( - FFcharge, SPEFFPed_intersection - ) - out["FFcharge_lg"] = ChargesMaker.select_charges_lg( - FFcharge, SPEFFPed_intersection - ) - out["Pedcharge_hg"] = ChargesMaker.select_charges_hg( - Pedcharge, SPEFFPed_intersection - ) - out["Pedcharge_lg"] = ChargesMaker.select_charges_lg( - Pedcharge, SPEFFPed_intersection - ) - - log.info(f"data have {len(SPEFFPed_intersection)} pixels in common") - return out - - @staticmethod - def __readFF(FFRun: int, **kwargs) -> dict: - """ - Reads FF charge data from a FITS file. - Args: - - FFRun (int): The run number for the FF data. - - kwargs (optional): Additional keyword arguments. - Returns: - - dict: A dictionary containing the FF charge data (`FFcharge`) and the - coefficient for the FF charge (`coefCharge_FF_Ped`). - """ - log.info("reading FF data") - method = kwargs.get("method", "FullWaveformSum") - FFchargeExtractorWindowLength = kwargs.get( - "FFchargeExtractorWindowLength", None - ) - if method != "FullWaveformSum": - if FFchargeExtractorWindowLength is None: - e = Exception( - "we have to specify FFchargeExtractorWindowLength argument if " - "charge extractor method is not FullwaveformSum" - ) - log.error(e, exc_info=True) - raise e - else: - coefCharge_FF_Ped = FFchargeExtractorWindowLength / constants.N_SAMPLES - else: - coefCharge_FF_Ped = 1 - if isinstance(FFRun, int): - try: - FFcharge = ChargesContainerIO.load( - f"{os.environ['NECTARCAMDATA']}/charges/{method}", FFRun - ) - log.info(f"charges have ever been computed for FF run {FFRun}") - except Exception as e: - log.error("charge have not been yet computed") - raise e - else: - e = TypeError("FFRun must be int") - log.error(e, exc_info=True) - raise e - return {"FFcharge": FFcharge, "coefCharge_FF_Ped": coefCharge_FF_Ped} - - @staticmethod - def __readPed(PedRun: int, **kwargs) -> dict: - """ - Reads Ped charge data from a FITS file. - Args: - - PedRun (int): The run number for the Ped data. - - kwargs (optional): Additional keyword arguments. - Returns: - - dict: A dictionary containing the Ped charge data (`Pedcharge`). - """ - log.info("reading Ped data") - method = "FullWaveformSum" # kwargs.get('method','std') - if isinstance(PedRun, int): - try: - Pedcharge = ChargesContainerIO.load( - f"{os.environ['NECTARCAMDATA']}/charges/{method}", PedRun - ) - log.info(f"charges have ever been computed for Ped run {PedRun}") - except Exception as e: - log.error("charge have not been yet computed") - raise e - else: - e = TypeError("PedRun must be int") - log.error(e, exc_info=True) - raise e - return {"Pedcharge": Pedcharge} - - def __check_shape(self) -> None: - """ - Checks the shape of certain attributes and raises an exception if the shape - is not as expected. - """ - try: - self.__FFcharge_hg[0] * self.__FFcharge_lg[0] * self.__Pedcharge_hg[ - 0 - ] * self.__Pedcharge_lg[0] * self.__SPE_resolution * self._pixels_id - except Exception as e: - log.error(e, exc_info=True) - raise e - - def make(self, **kwargs) -> None: - """ - Runs the photo statistic method and assigns values to the high_gain and - low_gain keys in the _results dictionary. - - Args: - **kwargs: Additional keyword arguments (not used in this method). - - Returns: - None - """ - log.info("running photo statistic method") - self._results["high_gain"] = self.gainHG - self._results["low_gain"] = self.gainLG - # self._results["is_valid"] = self._SPEvalid - - def plot_correlation( - photoStat_gain: np.ndarray, SPE_gain: np.ndarray - ) -> plt.Figure: - """ - Plot the correlation between the photo statistic gain and the single - photoelectron (SPE) gain. - - Args: - photoStat_gain (np.ndarray): Array of photo statistic gain values. - SPE_gain (np.ndarray): Array of SPE gain values. - - Returns: - fig (plt.Figure): The figure object containing the scatter plot and the - linear fit line. - """ - - # Create a mask to filter the data points based on certain criteria - mask = (photoStat_gain > 20) * (SPE_gain > 0) * (photoStat_gain < 80) - - # Perform a linear regression analysis on the filtered data points - a, b, r, p_value, std_err = linregress( - photoStat_gain[mask], SPE_gain[mask], "greater" - ) - - # Generate a range of x-values for the linear fit line - x = np.linspace(photoStat_gain[mask].min(), photoStat_gain[mask].max(), 1000) - - # Define a lambda function for the linear fit line - def y(x): - return a * x + b - - with quantity_support(): - # Create a scatter plot of the filtered data points - fig, ax = plt.subplots(1, 1, figsize=(8, 6)) - ax.scatter(photoStat_gain[mask], SPE_gain[mask], marker=".") - - # Plot the linear fit line using the x-values and the lambda function - ax.plot( - x, - y(x), - color="red", - label=f"linear fit,\n a = {a:.2e},\n b = {b:.2e},\n r = {r:.2e}," - f"\n p_value = {p_value:.2e},\n std_err = {std_err:.2e}", - ) - - # Plot the line y = x - ax.plot(x, x, color="black", label="y = x") - - ax.set_xlabel("Gain Photo stat (ADC)", size=15) - ax.set_ylabel("Gain SPE fit (ADC)", size=15) - ax.legend(fontsize=15) - - return fig - - -@property -def SPE_resolution(self) -> float: - """ - Returns a deep copy of the SPE resolution. - - Returns: - float: The SPE resolution. - """ - return copy.deepcopy(self.__SPE_resolution) - - -@property -def sigmaPedHG(self) -> float: - """ - Calculates and returns the standard deviation of Pedcharge_hg multiplied by the - square root of coefCharge_FF_Ped. - - Returns: - float: The standard deviation of Pedcharge_hg. - """ - return np.std(self.__Pedcharge_hg, axis=0) * np.sqrt(self.__coefCharge_FF_Ped) - - -@property -def sigmaChargeHG(self) -> float: - """ - Calculates and returns the standard deviation of FFcharge_hg minus meanPedHG. - - Returns: - float: The standard deviation of FFcharge_hg minus meanPedHG. - """ - return np.std(self.__FFcharge_hg - self.meanPedHG, axis=0) - - -@property -def meanPedHG(self) -> float: - """ - Calculates and returns the mean of Pedcharge_hg multiplied by coefCharge_FF_Ped. - - Returns: - float: The mean of Pedcharge_hg. - """ - return np.mean(self.__Pedcharge_hg, axis=0) * self.__coefCharge_FF_Ped - - -@property -def meanChargeHG(self) -> float: - """ - Calculates and returns the mean of FFcharge_hg minus meanPedHG. - - Returns: - float: The mean of FFcharge_hg minus meanPedHG. - """ - return np.mean(self.__FFcharge_hg - self.meanPedHG, axis=0) - - -@property -def BHG(self) -> float: - """ - Calculates and returns the BHG value. - - Returns: - float: The BHG value. - """ - min_events = np.min((self.__FFcharge_hg.shape[0], self.__Pedcharge_hg.shape[0])) - upper = ( - np.power( - self.__FFcharge_hg.mean(axis=1)[:min_events] - - self.__Pedcharge_hg.mean(axis=1)[:min_events] * self.__coefCharge_FF_Ped - - self.meanChargeHG.mean(), - 2, - ) - ).mean(axis=0) - lower = np.power(self.meanChargeHG.mean(), 2) - return np.sqrt(upper / lower) - - -@property -def gainHG(self) -> float: - """ - Calculates and returns the gain for high gain charge data. - - Returns: - float: The gain for high gain charge data. - """ - return ( - np.power(self.sigmaChargeHG, 2) - - np.power(self.sigmaPedHG, 2) - - np.power(self.BHG * self.meanChargeHG, 2) - ) / (self.meanChargeHG * (1 + np.power(self.SPE_resolution, 2))) - - -@property -def sigmaPedLG(self) -> float: - """ - Calculates and returns the standard deviation of Pedcharge_lg multiplied by the - square root of coefCharge_FF_Ped. - - Returns: - float: The standard deviation of Pedcharge_lg. - """ - return np.std(self.__Pedcharge_lg, axis=0) * np.sqrt(self.__coefCharge_FF_Ped) - - -@property -def sigmaChargeLG(self) -> float: - """ - Calculates and returns the standard deviation of FFcharge_lg minus meanPedLG. - - Returns: - float: The standard deviation of FFcharge_lg minus meanPedLG. - """ - return np.std(self.__FFcharge_lg - self.meanPedLG, axis=0) - - -@property -def meanPedLG(self) -> float: - """ - Calculates and returns the mean of Pedcharge_lg multiplied by coefCharge_FF_Ped. - - Returns: - float: The mean of Pedcharge_lg. - """ - return np.mean(self.__Pedcharge_lg, axis=0) * self.__coefCharge_FF_Ped - - -@property -def meanChargeLG(self) -> float: - """ - Calculates and returns the mean of FFcharge_lg minus meanPedLG. - - Returns: - float: The mean of FFcharge_lg minus meanPedLG. - """ - return np.mean(self.__FFcharge_lg - self.meanPedLG, axis=0) - - -@property -def BLG(self) -> float: - """ - Calculates and returns the BLG value. - - Returns: - float: The BLG value. - """ - min_events = np.min((self.__FFcharge_lg.shape[0], self.__Pedcharge_lg.shape[0])) - upper = ( - np.power( - self.__FFcharge_lg.mean(axis=1)[:min_events] - - self.__Pedcharge_lg.mean(axis=1)[:min_events] * self.__coefCharge_FF_Ped - - self.meanChargeLG.mean(), - 2, - ) - ).mean(axis=0) - lower = np.power(self.meanChargeLG.mean(), 2) - return np.sqrt(upper / lower) - - -@property -def gainLG(self) -> float: - """ - Calculates and returns the gain for low gain charge data. - - Returns: - float: The gain for low gain charge data. - """ - return ( - np.power(self.sigmaChargeLG, 2) - - np.power(self.sigmaPedLG, 2) - - np.power(self.BLG * self.meanChargeLG, 2) - ) / (self.meanChargeLG * (1 + np.power(self.SPE_resolution, 2))) diff --git a/src/nectarchain/makers/calibration/gain/utils/__init__.py b/src/nectarchain/makers/calibration/gain/utils/__init__.py deleted file mode 100644 index 6f3aa3a4..00000000 --- a/src/nectarchain/makers/calibration/gain/utils/__init__.py +++ /dev/null @@ -1,2 +0,0 @@ -from .error import * -from .utils import * diff --git a/src/nectarchain/makers/calibration/pedestal_makers.py b/src/nectarchain/makers/calibration/pedestalMakers.py similarity index 51% rename from src/nectarchain/makers/calibration/pedestal_makers.py rename to src/nectarchain/makers/calibration/pedestalMakers.py index 582ccc7f..ec0be630 100644 --- a/src/nectarchain/makers/calibration/pedestal_makers.py +++ b/src/nectarchain/makers/calibration/pedestalMakers.py @@ -1,20 +1,16 @@ import logging -from .core import CalibrationMaker - logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") log = logging.getLogger(__name__) log.handlers = logging.getLogger("__main__").handlers -__all__ = ["PedestalMaker"] +from .core import NectarCAMCalibrationTool +__all__ = ["PedestalNectarCAMCalibrationTool"] -class PedestalMaker(CalibrationMaker): - def __init__(self, *args, **kwargs) -> None: - super().__init__(*args, **kwargs) - def make(self): +class PedestalNectarCAMCalibrationTool(NectarCAMCalibrationTool): + def start(self): raise NotImplementedError( - "The computation of the pedestal calibration is not yet implemented, " - "feel free to contribute !:)" + "The computation of the pedestal calibration is not yet implemented, feel free to contribute !:)" ) diff --git a/src/nectarchain/makers/chargesMakers.py b/src/nectarchain/makers/chargesMakers.py new file mode 100644 index 00000000..97ff9cb2 --- /dev/null +++ b/src/nectarchain/makers/chargesMakers.py @@ -0,0 +1,144 @@ +import logging + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + +import os +import pathlib + +import numpy as np +from ctapipe.core.traits import Bool, ComponentNameList +from ctapipe.image.extractor import ( + BaselineSubtractedNeighborPeakWindowSum, + FixedWindowSum, + FullWaveformSum, + GlobalPeakWindowSum, + LocalPeakWindowSum, + NeighborPeakWindowSum, + SlidingWindowMaxSum, + TwoPassWindowSum, +) + +from ..data.container import ChargesContainers, WaveformsContainer, WaveformsContainers +from ..data.management import DataManagement +from .component import ChargesComponent, NectarCAMComponent +from .core import EventsLoopNectarCAMCalibrationTool +from .extractor.utils import CtapipeExtractor + +__all__ = ["ChargesNectarCAMCalibrationTool"] + + +class ChargesNectarCAMCalibrationTool(EventsLoopNectarCAMCalibrationTool): + """class use to make the waveform extraction from event read from r0 data""" + + name = "ChargesNectarCAMCalibration" + + componentsList = ComponentNameList( + NectarCAMComponent, + default_value=["ChargesComponent"], + help="List of Component names to be apply, the order will be respected", + ).tag(config=True) + + from_computed_waveforms = Bool( + default_value=False, + help="a flag to compute charge from waveforms stored on disk", + ) + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + def _init_output_path(self): + str_extractor_kwargs = CtapipeExtractor.get_extractor_kwargs_str( + self.extractor_kwargs + ) + if self.max_events is None: + filename = f"{self.name}_run{self.run_number}_{self.method}_{str_extractor_kwargs}.h5" + else: + filename = f"{self.name}_run{self.run_number}_maxevents{self.max_events}_{self.method}_{str_extractor_kwargs}.h5" + + self.output_path = pathlib.Path( + f"{os.environ.get('NECTARCAMDATA','/tmp')}/runs/charges/{filename}" + ) + + def start( + self, + n_events=np.inf, + # trigger_type: list = None, + restart_from_begining: bool = False, + *args, + **kwargs, + ): + if self.from_computed_waveforms: + files = DataManagement.find_waveforms( + run_number=self.run_number, max_events=self.max_events + ) + if len(files) != 1: + self.log.info( + f"{len(files)} computed wavforms files found with max_events >= {self.max_events} for run {self.run_number}, reload waveforms from event loop" + ) + super().start( + n_events=n_events, + restart_from_begining=restart_from_begining, + *args, + **kwargs, + ) + else: + self.log.info( + f"{files[0]} is the computed wavforms files found with max_events >= {self.max_events} for run {self.run_number}" + ) + waveformsContainers = WaveformsContainer.from_hdf5(files[0]) + if not (isinstance(waveformsContainers, WaveformsContainer)): + chargesContainers = ChargesContainers() + if isinstance(waveformsContainers, WaveformsContainers): + self.log.debug( + "WaveformsContainer file container multiple trigger type" + ) + self._init_writer(sliced=False) + chargesContainers = ( + ChargesComponent._create_from_waveforms_looping_eventType( + waveformsContainers=waveformsContainers, + subarray=self.event_source.subarray, + method=self.method, + **self.extractor_kwargs, + ) + ) + self._write_container(container=chargesContainers) + else: + self.log.debug( + "WaveformsContainer file container multiple slices of the run events" + ) + for slice_index, _waveformsContainers in enumerate( + waveformsContainers + ): + self._init_writer(sliced=True, slice_index=slice_index) + chargesContainers = ChargesComponent._create_from_waveforms_looping_eventType( + waveformsContainers=_waveformsContainers, + subarray=self.event_source.subarray, + method=self.method, + **self.extractor_kwargs, + ) + self._write_container(container=chargesContainers) + else: + self.log.debug( + "WaveformsContainer file container is a simple WaveformsContainer (not mapped)" + ) + self._init_writer(sliced=False) + chargesContainers = ChargesComponent.create_from_waveforms( + waveformsContainer=waveformsContainers, + subarray=self.event_source.subarray, + method=self.method, + **self.extractor_kwargs, + ) + self._write_container(container=chargesContainers) + + self.writer.close() + # self.__still_finished = True + + else: + super().start( + n_events=n_events, + restart_from_begining=restart_from_begining, + *args, + **kwargs, + ) diff --git a/src/nectarchain/makers/charges_makers.py b/src/nectarchain/makers/charges_makers.py deleted file mode 100644 index 9575935d..00000000 --- a/src/nectarchain/makers/charges_makers.py +++ /dev/null @@ -1,22 +0,0 @@ -import logging - -from ctapipe.core.traits import ComponentNameList - -from .component.core import NectarCAMComponent -from .core import EventsLoopNectarCAMCalibrationTool - -logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") -log = logging.getLogger(__name__) -log.handlers = logging.getLogger("__main__").handlers - -__all__ = ["ChargesNectarCAMCalibrationTool"] - - -class ChargesNectarCAMCalibrationTool(EventsLoopNectarCAMCalibrationTool): - """class use to make the waveform extraction from event read from r0 data""" - - componentsList = ComponentNameList( - NectarCAMComponent, - default_value=["ChargesComponent"], - help="List of Component names to be apply, the order will be respected", - ).tag(config=True) diff --git a/src/nectarchain/makers/component/FlatFieldSPEComponent.py b/src/nectarchain/makers/component/FlatFieldSPEComponent.py new file mode 100644 index 00000000..4c72452f --- /dev/null +++ b/src/nectarchain/makers/component/FlatFieldSPEComponent.py @@ -0,0 +1,219 @@ +import logging + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + +import copy + +import numpy as np +from ctapipe.core.traits import List, Unicode +from ctapipe_io_nectarcam.containers import NectarCAMDataContainer + +from ...data.container import merge_map_ArrayDataContainer +from ...utils import ComponentUtils +from .chargesComponent import ChargesComponent +from .gainComponent import GainNectarCAMComponent +from .spe import ( + SPECombinedalgorithm, + SPEHHValgorithm, + SPEHHVStdalgorithm, + SPEnominalalgorithm, + SPEnominalStdalgorithm, +) + +__all__ = [ + "FlatFieldSingleNominalSPEStdNectarCAMComponent", + "FlatFieldSingleNominalSPENectarCAMComponent", + "FlatFieldSingleHHVSPENectarCAMComponent", + "FlatFieldSingleHHVSPEStdNectarCAMComponent", + "FlatFieldCombinedSPEStdNectarCAMComponent", +] + + +class FlatFieldSingleNominalSPENectarCAMComponent(GainNectarCAMComponent): + SPEfitalgorithm = Unicode( + "SPEnominalalgorithm", + help="The Spe fit method to be use", + read_only=True, + ).tag(config=True) + + SubComponents = copy.deepcopy(GainNectarCAMComponent.SubComponents) + SubComponents.default_value = [ + "ChargesComponent", + f"{SPEfitalgorithm.default_value}", + ] + SubComponents.read_only = True + + # Windows_lenght = Integer(40, + # read_only = True, + # help = "The windows leght used for the savgol filter algorithm", + # ).tag(config = True) + # + # Order = Integer(2, + # read_only = True, + # help = "The order of the polynome used in the savgol filter algorithm", + # ).tag(config = True) + + asked_pixels_id = List( + default_value=None, + allow_none=True, + help="The pixels id where we want to perform the SPE fit", + ).tag(config=True) + + # nproc = Integer(8, + # help = "The Number of cpu used for SPE fit", + # ).tag(config = True) + # + # chunksize = Integer(1, + # help = "The chunk size for multi-processing", + # ).tag(config = True) + # + # multiproc = Bool(True, + # help = "flag to active multi-processing", + # ).tag(config = True) + + # method = Unicode(default_value = "FullWaveformSum", + # help = "the charge extraction method", + # + # ).tag(config = True) + # + # extractor_kwargs = Dict(default_value = {}, + # help = "The kwargs to be pass to the charge extractor method", + # ).tag(config = True) + + # constructor + def __init__(self, subarray, config=None, parent=None, *args, **kwargs) -> None: + chargesComponent_kwargs = {} + self._SPEfitalgorithm_kwargs = {} + other_kwargs = {} + chargesComponent_configurable_traits = ComponentUtils.get_configurable_traits( + ChargesComponent + ) + SPEfitalgorithm_configurable_traits = ComponentUtils.get_configurable_traits( + eval(self.SPEfitalgorithm) + ) + + for key in kwargs.keys(): + if key in chargesComponent_configurable_traits.keys(): + chargesComponent_kwargs[key] = kwargs[key] + elif key in SPEfitalgorithm_configurable_traits.keys(): + self._SPEfitalgorithm_kwargs[key] = kwargs[key] + else: + other_kwargs[key] = kwargs[key] + + super().__init__( + subarray=subarray, config=config, parent=parent, *args, **other_kwargs + ) + self.chargesComponent = ChargesComponent( + subarray=subarray, + config=config, + parent=parent, + *args, + **chargesComponent_kwargs, + ) + self._chargesContainers = None + + def __call__(self, event: NectarCAMDataContainer, *args, **kwargs): + self.chargesComponent(event=event, *args, **kwargs) + + def finish(self, *args, **kwargs): + is_empty = False + if self._chargesContainers is None: + self._chargesContainers = self.chargesComponent.finish(*args, **kwargs) + is_empty = self._chargesContainers.is_empty() + if is_empty: + log.warning("empty chargesContainer in output") + else: + self._chargesContainers = merge_map_ArrayDataContainer( + self._chargesContainers + ) + + if not (is_empty): + spe_fit = eval(self.SPEfitalgorithm).create_from_chargesContainer( + self._chargesContainers, parent=self, **self._SPEfitalgorithm_kwargs + ) + fit_output = spe_fit.run(pixels_id=self.asked_pixels_id, *args, **kwargs) + n_asked_pix = ( + len(self._chargesContainers.pixels_id) + if self.asked_pixels_id is None + else len(self.asked_pixels_id) + ) + conv_rate = np.sum(spe_fit.results.is_valid) / n_asked_pix + self.log.info(f"convergence rate : {conv_rate}") + return spe_fit.results + else: + return None + + +class FlatFieldSingleNominalSPEStdNectarCAMComponent( + FlatFieldSingleNominalSPENectarCAMComponent +): + SPEfitalgorithm = Unicode( + "SPEnominalStdalgorithm", + help="The Spe fit method to be use", + read_only=True, + ).tag(config=True) + + SubComponents = copy.deepcopy(GainNectarCAMComponent.SubComponents) + SubComponents.default_value = [ + "ChargesComponent", + f"{SPEfitalgorithm.default_value}", + ] + SubComponents.read_only = True + + +class FlatFieldSingleHHVSPENectarCAMComponent( + FlatFieldSingleNominalSPENectarCAMComponent +): + SPEfitalgorithm = Unicode( + "SPEHHValgorithm", + help="The Spe fit method to be use", + read_only=True, + ).tag(config=True) + + SubComponents = copy.deepcopy(GainNectarCAMComponent.SubComponents) + SubComponents.default_value = [ + "ChargesComponent", + f"{SPEfitalgorithm.default_value}", + ] + SubComponents.read_only = True + + +class FlatFieldSingleHHVSPEStdNectarCAMComponent( + FlatFieldSingleNominalSPENectarCAMComponent +): + SPEfitalgorithm = Unicode( + "SPEHHVStdalgorithm", + help="The Spe fit method to be use", + read_only=True, + ).tag(config=True) + + SubComponents = copy.deepcopy(GainNectarCAMComponent.SubComponents) + SubComponents.default_value = [ + "ChargesComponent", + f"{SPEfitalgorithm.default_value}", + ] + SubComponents.read_only = True + + +class FlatFieldCombinedSPEStdNectarCAMComponent( + FlatFieldSingleNominalSPENectarCAMComponent +): + SPEfitalgorithm = Unicode( + "SPECombinedalgorithm", + help="The Spe fit method to be use", + read_only=True, + ).tag(config=True) + + SubComponents = copy.deepcopy(GainNectarCAMComponent.SubComponents) + SubComponents.default_value = [ + "ChargesComponent", + f"{SPEfitalgorithm.default_value}", + ] + SubComponents.read_only = True + + def __init__(self, subarray, config=None, parent=None, *args, **kwargs) -> None: + super().__init__( + subarray=subarray, config=config, parent=parent, *args, **kwargs + ) diff --git a/src/nectarchain/makers/component/__init__.py b/src/nectarchain/makers/component/__init__.py index e69de29b..8ab6a9a6 100644 --- a/src/nectarchain/makers/component/__init__.py +++ b/src/nectarchain/makers/component/__init__.py @@ -0,0 +1,25 @@ +from .chargesComponent import * +from .core import * +from .FlatFieldSPEComponent import * +from .gainComponent import * +from .photostatistic_algorithm import * +from .photostatistic_component import * +from .spe import * +from .waveformsComponent import * + +__all__ = [ + "ArrayDataComponent", + "NectarCAMComponent", + "SPEHHValgorithm", + "SPEHHVStdalgorithm", + "SPECombinedalgorithm", + "FlatFieldSingleHHVSPENectarCAMComponent", + "FlatFieldSingleHHVSPEStdNectarCAMComponent", + "FlatFieldSingleNominalSPENectarCAMComponent", + "FlatFieldSingleNominalSPEStdNectarCAMComponent", + "FlatFieldCombinedSPEStdNectarCAMComponent", + "ChargesComponent", + "WaveformsComponent", + "PhotoStatisticNectarCAMComponent", + "PhotoStatisticAlgorithm", +] diff --git a/src/nectarchain/makers/component/charges_component.py b/src/nectarchain/makers/component/chargesComponent.py similarity index 80% rename from src/nectarchain/makers/component/charges_component.py rename to src/nectarchain/makers/component/chargesComponent.py index 458e17be..9a62b11e 100644 --- a/src/nectarchain/makers/component/charges_component.py +++ b/src/nectarchain/makers/component/chargesComponent.py @@ -1,4 +1,10 @@ import logging + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + +import copy import time from argparse import ArgumentError @@ -6,19 +12,30 @@ import numpy.ma as ma from ctapipe.containers import EventType from ctapipe.core.traits import Dict, Unicode +from ctapipe.image.extractor import ( + BaselineSubtractedNeighborPeakWindowSum, + FixedWindowSum, + FullWaveformSum, + GlobalPeakWindowSum, + LocalPeakWindowSum, + NeighborPeakWindowSum, + SlidingWindowMaxSum, + TwoPassWindowSum, +) from ctapipe.instrument import SubarrayDescription from ctapipe_io_nectarcam import constants from ctapipe_io_nectarcam.containers import NectarCAMDataContainer from numba import bool_, float64, guvectorize, int64 -from ...data.container import ChargesContainer, ChargesContainers, WaveformsContainer +from ...data.container import ( + ChargesContainer, + ChargesContainers, + WaveformsContainer, + WaveformsContainers, +) from ..extractor.utils import CtapipeExtractor from .core import ArrayDataComponent -logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") -log = logging.getLogger(__name__) -log.handlers = logging.getLogger("__main__").handlers - __all__ = ["ChargesComponent"] list_ctapipe_charge_extractor = [ @@ -105,6 +122,9 @@ class ChargesComponent(ArrayDataComponent): help="The kwargs to be pass to the charge extractor method", ).tag(config=True) + SubComponents = copy.deepcopy(ArrayDataComponent.SubComponents) + SubComponents.read_only = True + def __init__(self, subarray, config=None, parent=None, *args, **kwargs): super().__init__( subarray=subarray, config=config, parent=parent, *args, **kwargs @@ -171,41 +191,43 @@ def __call__( self.__charges_lg[f"{name}"].append(__image[0].tolist()) self.__peak_lg[f"{name}"].append(__image[1].tolist()) + @staticmethod + def _get_extractor_kwargs_from_method_and_kwargs(method: str, kwargs: dict): + extractor_kwargs = {} + for key in eval(method).class_own_traits().keys(): + if key in kwargs.keys(): + extractor_kwargs[key] = kwargs[key] + if ( + "apply_integration_correction" in eval(method).class_own_traits().keys() + ): # to change the default behavior of ctapipe extractor + extractor_kwargs["apply_integration_correction"] = kwargs.get( + "apply_integration_correction", False + ) + return extractor_kwargs + @staticmethod def _get_imageExtractor(method: str, subarray: SubarrayDescription, **kwargs): """ - Create an instance of a charge extraction method based on the provided method - name and subarray description. + Create an instance of a charge extraction method based on the provided method name and subarray description. Args: method (str): The name of the charge extraction method. subarray (SubarrayDescription): The description of the subarray. - **kwargs (dict): Additional keyword arguments for the charge extraction - method. + **kwargs (dict): Additional keyword arguments for the charge extraction method. Returns: - imageExtractor: An instance of the charge extraction method specified by - `method` with the provided subarray description and keyword arguments. + imageExtractor: An instance of the charge extraction method specified by `method` with the provided subarray description and keyword arguments. """ if not ( method in list_ctapipe_charge_extractor or method in list_nectarchain_charge_extractor ): raise ArgumentError( - f"method must be in {list_ctapipe_charge_extractor} or" - f" {list_nectarchain_charge_extractor}" - ) - extractor_kwargs = {} - for key in eval(method).class_own_traits().keys(): - if key in kwargs.keys(): - extractor_kwargs[key] = kwargs[key] - if ( - "apply_integration_correction" in eval(method).class_own_traits().keys() - ): # to change the default behavior of ctapipe extractor - extractor_kwargs["apply_integration_correction"] = kwargs.get( - "apply_integration_correction", False + f"method must be in {list_ctapipe_charge_extractor} or {list_nectarchain_charge_extractor}" ) + extractor_kwargs = __class__._get_extractor_kwargs_from_method_and_kwargs( + method=method, kwargs=kwargs + ) log.debug( - f"Extracting charges with method {method} and extractor_kwargs" - f" {extractor_kwargs}" + f"Extracting charges with method {method} and extractor_kwargs {extractor_kwargs}" ) imageExtractor = eval(method)(subarray, **extractor_kwargs) return imageExtractor @@ -224,10 +246,12 @@ def finish(self, *args, **kwargs): output = ChargesContainers() for i, trigger in enumerate(self.trigger_list): chargesContainer = ChargesContainer( - run_number=self._run_number, - npixels=self._npixels, + run_number=ChargesContainer.fields["run_number"].type(self._run_number), + npixels=ChargesContainer.fields["npixels"].type(self._npixels), camera=self.CAMERA_NAME, - pixels_id=self._pixels_id, + pixels_id=ChargesContainer.fields["pixels_id"].dtype.type( + self._pixels_id + ), method=self.method, nevents=self.nevents(trigger), charges_hg=self.charges_hg(trigger), @@ -253,12 +277,10 @@ def sort(chargesContainer: ChargesContainer, method: str = "event_id"): """ Sorts the charges in a ChargesContainer object based on the specified method. Args: - chargesContainer (ChargesContainer): The ChargesContainer object to be - sorted. + chargesContainer (ChargesContainer): The ChargesContainer object to be sorted. method (str, optional): The sorting method. Defaults to 'event_id'. Returns: - ChargesContainer: A new ChargesContainer object with the charges sorted - based on the specified method. + ChargesContainer: A new ChargesContainer object with the charges sorted based on the specified method. Raises: ArgumentError: If the specified method is not valid. @@ -293,14 +315,12 @@ def sort(chargesContainer: ChargesContainer, method: str = "event_id"): @staticmethod def select_charges_hg(chargesContainer: ChargesContainer, pixel_id: np.ndarray): """ - Selects the charges from the ChargesContainer object for the given pixel_id and - returns the result transposed. + Selects the charges from the ChargesContainer object for the given pixel_id and returns the result transposed. Args: chargesContainer (ChargesContainer): The ChargesContainer object. pixel_id (np.ndarray): An array of pixel IDs. Returns: - np.ndarray: The selected charges from the ChargesContainer object for the - given pixel_id, transposed. + np.ndarray: The selected charges from the ChargesContainer object for the given pixel_id, transposed. """ res = __class__.select_container_array_field( container=chargesContainer, pixel_id=pixel_id, field="charges_hg" @@ -311,14 +331,12 @@ def select_charges_hg(chargesContainer: ChargesContainer, pixel_id: np.ndarray): @staticmethod def select_charges_lg(chargesContainer: ChargesContainer, pixel_id: np.ndarray): """ - Selects the charges from the ChargesContainer object for the given pixel_id - and returns the result transposed. + Selects the charges from the ChargesContainer object for the given pixel_id and returns the result transposed. Args: chargesContainer (ChargesContainer): The ChargesContainer object. pixel_id (np.ndarray): An array of pixel IDs. Returns: - np.ndarray: The selected charges from the ChargesContainer object for the - given pixel_id, transposed. + np.ndarray: The selected charges from the ChargesContainer object for the given pixel_id, transposed. """ res = __class__.select_container_array_field( container=chargesContainer, pixel_id=pixel_id, field="charges_lg" @@ -328,89 +346,103 @@ def select_charges_lg(chargesContainer: ChargesContainer, pixel_id: np.ndarray): def charges_hg(self, trigger: EventType): """ - Returns the charges for a specific trigger type as a NumPy array of unsigned - 16-bit integers. + Returns the charges for a specific trigger type as a NumPy array of unsigned 16-bit integers. Args: trigger (EventType): The specific trigger type. Returns: np.ndarray: The charges for the specific trigger type. """ return np.array( - self.__charges_hg[__class__._get_name_trigger(trigger)], dtype=np.uint16 + self.__charges_hg[__class__._get_name_trigger(trigger)], + dtype=ChargesContainer.fields["charges_hg"].dtype, ) def charges_lg(self, trigger: EventType): """ - Returns the charges for a specific trigger type as a NumPy array of unsigned - 16-bit integers. + Returns the charges for a specific trigger type as a NumPy array of unsigned 16-bit integers. Args: trigger (EventType): The specific trigger type. Returns: np.ndarray: The charges for the specific trigger type. """ return np.array( - self.__charges_lg[__class__._get_name_trigger(trigger)], dtype=np.uint16 + self.__charges_lg[__class__._get_name_trigger(trigger)], + dtype=ChargesContainer.fields["charges_lg"].dtype, ) def peak_hg(self, trigger: EventType): """ - Returns the peak charges for a specific trigger type as a NumPy array of - unsigned 16-bit integers. + Returns the peak charges for a specific trigger type as a NumPy array of unsigned 16-bit integers. Args: trigger (EventType): The specific trigger type. Returns: np.ndarray: The peak charges for the specific trigger type. """ return np.array( - self.__peak_hg[__class__._get_name_trigger(trigger)], dtype=np.uint16 + self.__peak_hg[__class__._get_name_trigger(trigger)], + dtype=ChargesContainer.fields["peak_hg"].dtype, ) def peak_lg(self, trigger: EventType): """ - Returns the peak charges for a specific trigger type as a NumPy array of - unsigned 16-bit integers. + Returns the peak charges for a specific trigger type as a NumPy array of unsigned 16-bit integers. Args: trigger (EventType): The specific trigger type. Returns: np.ndarray: The peak charges for the specific trigger type. """ return np.array( - self.__peak_lg[__class__._get_name_trigger(trigger)], dtype=np.uint16 + self.__peak_lg[__class__._get_name_trigger(trigger)], + dtype=ChargesContainer.fields["peak_lg"].dtype, ) + @staticmethod + def _create_from_waveforms_looping_eventType( + waveformsContainers: WaveformsContainers, **kwargs + ): + chargesContainers = ChargesContainers() + for key in waveformsContainers.containers.keys(): + chargesContainers.containers[key] = ChargesComponent.create_from_waveforms( + waveformsContainer=waveformsContainers.containers[key], **kwargs + ) + return chargesContainers + @staticmethod def create_from_waveforms( waveformsContainer: WaveformsContainer, + subarray=SubarrayDescription, method: str = "FullWaveformSum", **kwargs, ) -> ChargesContainer: """ - Create a ChargesContainer object from waveforms using the specified charge - extraction method. + Create a ChargesContainer object from waveforms using the specified charge extraction method. Args: waveformsContainer (WaveformsContainer): The waveforms container object. - method (str, optional): The charge extraction method to use (default is - "FullWaveformSum"). - **kwargs: Additional keyword arguments to pass to the charge extraction - method. + method (str, optional): The charge extraction method to use (default is "FullWaveformSum"). + **kwargs: Additional keyword arguments to pass to the charge extraction method. Returns: - ChargesContainer: The charges container object containing the computed - charges and peak times. + ChargesContainer: The charges container object containing the computed charges and peak times. """ chargesContainer = ChargesContainer() for field in waveformsContainer.keys(): if not (field in ["nsamples", "wfs_hg", "wfs_lg"]): chargesContainer[field] = waveformsContainer[field] log.info(f"computing hg charge with {method} method") - charges_hg, peak_hg = __class__.compute_charge( - waveformsContainer, constants.HIGH_GAIN, method, **kwargs + charges_hg, peak_hg = __class__.compute_charges( + waveformsContainer=waveformsContainer, + channel=constants.HIGH_GAIN, + subarray=subarray, + method=method, + **kwargs, ) - charges_hg = np.array(charges_hg, dtype=np.uint16) log.info(f"computing lg charge with {method} method") - charges_lg, peak_lg = __class__.compute_charge( - waveformsContainer, constants.LOW_GAIN, method, **kwargs + charges_lg, peak_lg = __class__.compute_charges( + waveformsContainer=waveformsContainer, + channel=constants.LOW_GAIN, + subarray=subarray, + method=method, + **kwargs, ) - charges_lg = np.array(charges_lg, dtype=np.uint16) chargesContainer.charges_hg = charges_hg chargesContainer.charges_lg = charges_lg chargesContainer.peak_hg = peak_hg @@ -419,8 +451,8 @@ def create_from_waveforms( return chargesContainer @staticmethod - def compute_charge( - waveformContainer: WaveformsContainer, + def compute_charges( + waveformsContainer: WaveformsContainer, channel: int, subarray: SubarrayDescription, method: str = "FullWaveformSum", @@ -432,18 +464,15 @@ def compute_charge( Args: waveformContainer (WaveformsContainer): The waveforms container object. channel (int): The channel to compute charges for. - method (str, optional): The charge extraction method to use (default is - "FullWaveformSum"). - **kwargs: Additional keyword arguments to pass to the charge extraction - method. + method (str, optional): The charge extraction method to use (default is "FullWaveformSum"). + **kwargs: Additional keyword arguments to pass to the charge extraction method. Raises: ArgumentError: If the extraction method is unknown. ArgumentError: If the channel is unknown. Returns: tuple: A tuple containing the computed charges and peak times. """ - # import is here for fix issue with pytest (TypeError : inference is not - # possible with python <3.9 (Numba conflict bc there is no inference...)) + # import is here for fix issue with pytest (TypeError : inference is not possible with python <3.9 (Numba conflict bc there is no inference...)) from ..extractor.utils import CtapipeExtractor if tel_id is None: @@ -457,31 +486,35 @@ def compute_charge( [ CtapipeExtractor.get_image_peak_time( imageExtractor( - waveformContainer.wfs_hg[i], + waveformsContainer.wfs_hg[i], tel_id, channel, - waveformContainer.broken_pixels_hg, + waveformsContainer.broken_pixels_hg[i], ) ) - for i in range(len(waveformContainer.wfs_hg)) + for i in range(len(waveformsContainer.wfs_hg)) ] ).transpose(1, 0, 2) - return out[0], out[1] + return ChargesContainer.fields["charges_hg"].dtype.type( + out[0] + ), ChargesContainer.fields["peak_hg"].dtype.type(out[1]) elif channel == constants.LOW_GAIN: out = np.array( [ CtapipeExtractor.get_image_peak_time( imageExtractor( - waveformContainer.wfs_lg[i], + waveformsContainer.wfs_lg[i], tel_id, channel, - waveformContainer.broken_pixels_lg, + waveformsContainer.broken_pixels_lg[i], ) ) - for i in range(len(waveformContainer.wfs_lg)) + for i in range(len(waveformsContainer.wfs_lg)) ] ).transpose(1, 0, 2) - return out[0], out[1] + return ChargesContainer.fields["charges_lg"].dtype.type( + out[0] + ), ChargesContainer.fields["peak_lg"].dtype.type(out[1]) else: raise ArgumentError( f"channel must be {constants.LOW_GAIN} or {constants.HIGH_GAIN}" @@ -495,17 +528,12 @@ def histo_hg( Computes histogram of high gain charges from a ChargesContainer object. Args: - chargesContainer (ChargesContainer): A ChargesContainer object that holds - information about charges from a specific run. - n_bins (int, optional): The number of bins in the charge histogram. Defaults - to 1000. - autoscale (bool, optional): Whether to automatically detect the number of - bins based on the pixel data. Defaults to True. + chargesContainer (ChargesContainer): A ChargesContainer object that holds information about charges from a specific run. + n_bins (int, optional): The number of bins in the charge histogram. Defaults to 1000. + autoscale (bool, optional): Whether to automatically detect the number of bins based on the pixel data. Defaults to True. Returns: - ma.masked_array: A masked array representing the charge histogram, - where each row corresponds to an event and each column corresponds to a - bin in the histogram. + ma.masked_array: A masked array representing the charge histogram, where each row corresponds to an event and each column corresponds to a bin in the histogram. """ return __class__._histo( chargesContainer=chargesContainer, @@ -522,17 +550,12 @@ def histo_lg( Computes histogram of low gain charges from a ChargesContainer object. Args: - chargesContainer (ChargesContainer): A ChargesContainer object that holds - information about charges from a specific run. - n_bins (int, optional): The number of bins in the charge histogram. Defaults - to 1000. - autoscale (bool, optional): Whether to automatically detect the number of - bins based on the pixel data. Defaults to True. + chargesContainer (ChargesContainer): A ChargesContainer object that holds information about charges from a specific run. + n_bins (int, optional): The number of bins in the charge histogram. Defaults to 1000. + autoscale (bool, optional): Whether to automatically detect the number of bins based on the pixel data. Defaults to True. Returns: - ma.masked_array: A masked array representing the charge histogram, - where each row corresponds to an event and each column corresponds to a - bin in the histogram. + ma.masked_array: A masked array representing the charge histogram, where each row corresponds to an event and each column corresponds to a bin in the histogram. """ return __class__._histo( chargesContainer=chargesContainer, @@ -553,18 +576,13 @@ def _histo( Numba is used to compute histograms in a vectorized way. Args: - chargesContainer (ChargesContainer): A ChargesContainer object that holds - information about charges from a specific run. + chargesContainer (ChargesContainer): A ChargesContainer object that holds information about charges from a specific run. field (str): The field name for which the histogram is computed. - n_bins (int, optional): The number of bins in the charge histogram. Defaults - to 1000. - autoscale (bool, optional): Whether to automatically detect the number of - bins based on the pixel data. Defaults to True. + n_bins (int, optional): The number of bins in the charge histogram. Defaults to 1000. + autoscale (bool, optional): Whether to automatically detect the number of bins based on the pixel data. Defaults to True. Returns: - ma.masked_array: A masked array representing the charge histogram, - where each row corresponds to an event and each column corresponds to a - bin in the histogram. + ma.masked_array: A masked array representing the charge histogram, where each row corresponds to an event and each column corresponds to a bin in the histogram. """ mask_broken_pix = np.array( (chargesContainer[field] == chargesContainer[field].mean(axis=0)).mean( @@ -573,8 +591,7 @@ def _histo( dtype=bool, ) log.debug( - f"there are {mask_broken_pix.sum()} broken pixels (charge stays at same " - f"level for each events)" + f"there are {mask_broken_pix.sum()} broken pixels (charge stays at same level for each events)" ) if autoscale: diff --git a/src/nectarchain/makers/component/core.py b/src/nectarchain/makers/component/core.py index 08219e46..ee5b1470 100644 --- a/src/nectarchain/makers/component/core.py +++ b/src/nectarchain/makers/component/core.py @@ -1,27 +1,37 @@ -import copy import logging + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + +import copy from abc import abstractmethod +from collections.abc import Iterable import numpy as np from ctapipe.containers import EventType from ctapipe.core import Component, TelescopeComponent -from ctapipe.core.traits import Integer, Unicode +from ctapipe.core.traits import ComponentNameList, Integer, Unicode from ctapipe.instrument import CameraGeometry +from ctapipe.io import HDF5TableReader from ctapipe_io_nectarcam import constants from ctapipe_io_nectarcam.containers import NectarCAMDataContainer +from ...data.container import ( + ArrayDataContainer, + ChargesContainer, + ChargesContainers, + NectarCAMContainer, + TriggerMapContainer, + WaveformsContainer, + WaveformsContainers, +) from ...data.container.core import ArrayDataContainer -logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") -log = logging.getLogger(__name__) -log.handlers = logging.getLogger("__main__").handlers - __all__ = [ "ArrayDataComponent", "NectarCAMComponent", "get_valid_component", - "get_specific_traits", - "get_configurable_traits", ] @@ -29,25 +39,17 @@ def get_valid_component(): return NectarCAMComponent.non_abstract_subclasses() -def get_specific_traits(component: Component): - traits_dict = component.class_traits() - traits_dict.pop("config", True) - traits_dict.pop("parent", True) - return traits_dict - - -def get_configurable_traits(component: Component): - traits_dict = get_specific_traits(component) - output_traits_dict = traits_dict.copy() - for key, item in traits_dict.items(): - if item.read_only: - output_traits_dict.pop(key) - return output_traits_dict - - class NectarCAMComponent(TelescopeComponent): """The base class for NectarCAM components""" + SubComponents = ComponentNameList( + Component, + default_value=None, + allow_none=True, + read_only=True, + help="List of Component names that are used insite current component, this is used to resolve recusively the configurable traits defined in sub-components", + ).tag(config=True) + def __init__(self, subarray, config=None, parent=None, *args, **kwargs): super().__init__( subarray=subarray, config=config, parent=parent, *args, **kwargs @@ -85,11 +87,6 @@ def npixels(self): return copy.deepcopy(self.__npixels) -# class NectarCAMTelescopeComponent(TelescopeComponent) : -# """The base class for NectarCAM telescope component""" -# pass - - class ArrayDataComponent(NectarCAMComponent): TEL_ID = Integer( default_value=0, @@ -133,8 +130,7 @@ def _init_trigger_type(self, trigger: EventType, **kwargs): Initializes empty lists for different trigger types in the ArrayDataMaker class. Args: - trigger (EventType): The trigger type for which the lists are being - initialized. + trigger (EventType): The trigger type for which the lists are being initialized. Returns: None. The method only initializes the empty lists for the trigger type. @@ -159,10 +155,9 @@ def _compute_broken_pixels(wfs_hg, wfs_lg, **kwargs): wfs_lg (ndarray): Low gain waveforms. **kwargs: Additional keyword arguments. Returns: - tuple: Two arrays of zeros with the same shape as `wfs_hg` (or `wfs_lg`) - but without the last dimension. + tuple: Two arrays of zeros with the same shape as `wfs_hg` (or `wfs_lg`) but without the last dimension. """ - log.warning("computation of broken pixels is not yet implemented") + log.debug("computation of broken pixels is not yet implemented") return np.zeros((wfs_hg.shape[:-1]), dtype=bool), np.zeros( (wfs_hg.shape[:-1]), dtype=bool ) @@ -180,7 +175,7 @@ def _compute_broken_pixels_event( Returns: tuple: Two arrays of zeros with the length of `pixels_id`. """ - log.warning("computation of broken pixels is not yet implemented") + log.debug("computation of broken pixels is not yet implemented") return np.zeros((len(pixels_id)), dtype=bool), np.zeros( (len(pixels_id)), dtype=bool ) @@ -211,8 +206,7 @@ def __call__(self, event: NectarCAMDataContainer, *args, **kwargs): **kwargs: Additional keyword arguments that can be passed to the method. Returns: - If the return_wfs keyword argument is True, the method returns the high - and low gain waveforms from the event. + If the return_wfs keyword argument is True, the method returns the high and low gain waveforms from the event. """ name = __class__._get_name_trigger(event.trigger.event_type) @@ -248,14 +242,11 @@ def select_container_array_field( container: ArrayDataContainer, pixel_id: np.ndarray, field: str ) -> np.ndarray: """ - Selects specific fields from an ArrayDataContainer object based on a given - list of pixel IDs. + Selects specific fields from an ArrayDataContainer object based on a given list of pixel IDs. Args: - container (ArrayDataContainer): An object of type ArrayDataContainer that - contains the data. - pixel_id (ndarray): An array of pixel IDs for which the data needs to be - selected. + container (ArrayDataContainer): An object of type ArrayDataContainer that contains the data. + pixel_id (ndarray): An array of pixel IDs for which the data needs to be selected. field (str): The name of the field to be selected from the container. Returns: @@ -266,8 +257,7 @@ def select_container_array_field( ) for pixel in pixel_id[~mask_contain_pixels_id]: log.warning( - f"You asked for pixel_id {pixel} but it is not present in this " - f"container, skip this one" + f"You asked for pixel_id {pixel} but it is not present in this container, skip this one" ) res = np.array( [ @@ -279,10 +269,29 @@ def select_container_array_field( for pixel in pixel_id[mask_contain_pixels_id] ] ) - # could be nice to return np.ma.masked_array(data = res, mask = - # container.broken_pixels_hg.transpose(res.shape[1],res.shape[0],res.shape[2])) + ####could be nice to return np.ma.masked_array(data = res, mask = container.broken_pixels_hg.transpose(res.shape[1],res.shape[0],res.shape[2])) return res + @staticmethod + def merge_along_slices( + containers_generator: Iterable, + ) -> ArrayDataContainer: + for i, container in enumerate(containers_generator): + if i == 0: + merged_containers = copy.deepcopy(container) + else: + for trigger in container.containers.keys(): + if trigger in merged_containers.containers.keys(): + merged_containers.containers[trigger] = __class__.merge( + merged_containers.containers[trigger], + container.containers[trigger], + ) + else: + merged_containers.containers[trigger] = copy.deepcopy( + container.containers[trigger] + ) + return merged_containers + @staticmethod def merge( container_a: ArrayDataContainer, container_b: ArrayDataContainer @@ -295,26 +304,31 @@ def merge( if type(container_a) != type(container_b): raise Exception("The containers have to be instnace of the same class") - if np.array_equal(container_a.pixels_id, container_b.pixels_id): + if not (np.array_equal(container_a.pixels_id, container_b.pixels_id)): raise Exception("The containers have not the same pixels ids") - merged_container = container_a.__class__.__new__() + merged_container = container_a.__class__() for field in container_a.keys(): - if ~isinstance(container_a[field], np.ndarray): - if container_a[field] != container_b[field]: + if not (isinstance(container_a[field], np.ndarray)): + if field != "nevents" and (container_a[field] != container_b[field]): raise Exception( - f"merge impossible because of {field} field (values are" - f" {container_a[field]} and {container_b[field]}" + f"merge impossible because of {field} filed (values are {container_a[field]} and {container_b[field]}" ) for field in container_a.keys(): if isinstance(container_a[field], np.ndarray): - merged_container[field] = np.concatenate( - container_a[field], container_a[field], axis=0 - ) + if field != "pixels_id": + merged_container[field] = np.concatenate( + (container_a[field], container_b[field]), axis=0 + ) + else: + merged_container[field] = container_a[field] else: - merged_container[field] = container_a[field] + if field == "nevents": + merged_container[field] = container_a[field] + container_b[field] + else: + merged_container[field] = container_a[field] return merged_container @@ -343,13 +357,14 @@ def nevents(self, trigger: EventType): Returns the number of events for the specified trigger type. Args: - trigger (EventType): The trigger type for which the number of events is - requested. + trigger (EventType): The trigger type for which the number of events is requested. Returns: int: The number of events for the specified trigger type. """ - return len(self.__event_id[__class__._get_name_trigger(trigger)]) + return ArrayDataContainer.fields["nevents"].type( + len(self.__event_id[__class__._get_name_trigger(trigger)]) + ) @property def _broken_pixels_hg(self): @@ -366,15 +381,14 @@ def broken_pixels_hg(self, trigger: EventType): Returns an array of broken pixels for high gain for the specified trigger type. Args: - trigger (EventType): The trigger type for which the broken pixels for - high gain are requested. + trigger (EventType): The trigger type for which the broken pixels for high gain are requested. Returns: - np.ndarray: An array of broken pixels for high gain for the specified - trigger type. + np.ndarray: An array of broken pixels for high gain for the specified trigger type. """ return np.array( - self.__broken_pixels_hg[__class__._get_name_trigger(trigger)], dtype=bool + self.__broken_pixels_hg[__class__._get_name_trigger(trigger)], + dtype=ArrayDataContainer.fields["broken_pixels_hg"].dtype, ) @property @@ -392,15 +406,14 @@ def broken_pixels_lg(self, trigger: EventType): Returns an array of broken pixels for low gain for the specified trigger type. Args: - trigger (EventType): The trigger type for which the broken pixels for low - gain are requested. + trigger (EventType): The trigger type for which the broken pixels for low gain are requested. Returns: - np.ndarray: An array of broken pixels for low gain for the specified - trigger type. + np.ndarray: An array of broken pixels for low gain for the specified trigger type. """ return np.array( - self.__broken_pixels_lg[__class__._get_name_trigger(trigger)], dtype=bool + self.__broken_pixels_lg[__class__._get_name_trigger(trigger)], + dtype=ArrayDataContainer.fields["broken_pixels_lg"].dtype, ) def ucts_timestamp(self, trigger: EventType): @@ -408,14 +421,14 @@ def ucts_timestamp(self, trigger: EventType): Returns an array of UCTS timestamps for the specified trigger type. Args: - trigger (EventType): The trigger type for which the UCTS timestamps are - requested. + trigger (EventType): The trigger type for which the UCTS timestamps are requested. Returns: np.ndarray: An array of UCTS timestamps for the specified trigger type. """ return np.array( - self.__ucts_timestamp[__class__._get_name_trigger(trigger)], dtype=np.uint64 + self.__ucts_timestamp[__class__._get_name_trigger(trigger)], + dtype=ArrayDataContainer.fields["ucts_timestamp"].dtype, ) def ucts_busy_counter(self, trigger: EventType): @@ -423,15 +436,14 @@ def ucts_busy_counter(self, trigger: EventType): Returns an array of UCTS busy counters for the specified trigger type. Args: - trigger (EventType): The trigger type for which the UCTS busy counters - are requested. + trigger (EventType): The trigger type for which the UCTS busy counters are requested. Returns: np.ndarray: An array of UCTS busy counters for the specified trigger type. """ return np.array( self.__ucts_busy_counter[__class__._get_name_trigger(trigger)], - dtype=np.uint32, + dtype=ArrayDataContainer.fields["ucts_busy_counter"].dtype, ) def ucts_event_counter(self, trigger: EventType): @@ -439,15 +451,14 @@ def ucts_event_counter(self, trigger: EventType): Returns an array of UCTS event counters for the specified trigger type. Args: - trigger (EventType): The trigger type for which the UCTS event counters - are requested. + trigger (EventType): The trigger type for which the UCTS event counters are requested. Returns: np.ndarray: An array of UCTS event counters for the specified trigger type. """ return np.array( self.__ucts_event_counter[__class__._get_name_trigger(trigger)], - dtype=np.uint32, + dtype=ArrayDataContainer.fields["ucts_event_counter"].dtype, ) def event_type(self, trigger: EventType): @@ -455,14 +466,14 @@ def event_type(self, trigger: EventType): Returns an array of event types for the specified trigger type. Args: - trigger (EventType): The trigger type for which the event types are - requested. + trigger (EventType): The trigger type for which the event types are requested. Returns: np.ndarray: An array of event types for the specified trigger type. """ return np.array( - self.__event_type[__class__._get_name_trigger(trigger)], dtype=np.uint8 + self.__event_type[__class__._get_name_trigger(trigger)], + dtype=ArrayDataContainer.fields["event_type"].dtype, ) def event_id(self, trigger: EventType): @@ -476,7 +487,8 @@ def event_id(self, trigger: EventType): np.ndarray: An array of event IDs for the specified trigger type. """ return np.array( - self.__event_id[__class__._get_name_trigger(trigger)], dtype=np.uint32 + self.__event_id[__class__._get_name_trigger(trigger)], + dtype=ArrayDataContainer.fields["event_id"].dtype, ) def multiplicity(self, trigger: EventType): @@ -484,8 +496,7 @@ def multiplicity(self, trigger: EventType): Returns an array of multiplicities for the specified trigger type. Args: - trigger (EventType): The trigger type for which the multiplicities are - requested. + trigger (EventType): The trigger type for which the multiplicities are requested. Returns: np.ndarray: An array of multiplicities for the specified trigger type. @@ -494,15 +505,16 @@ def multiplicity(self, trigger: EventType): if len(tmp) == 0: return np.array([]) else: - return np.uint16(np.count_nonzero(tmp, axis=1)) + return ArrayDataContainer.fields["multiplicity"].dtype.type( + np.count_nonzero(tmp, axis=1) + ) def trig_pattern(self, trigger: EventType): """ Returns an array of trigger patterns for the specified trigger type. Args: - trigger (EventType): The trigger type for which the trigger patterns are - requested. + trigger (EventType): The trigger type for which the trigger patterns are requested. Returns: np.ndarray: An array of trigger patterns for the specified trigger type. @@ -515,18 +527,15 @@ def trig_pattern(self, trigger: EventType): def trig_pattern_all(self, trigger: EventType): """ - Returns an array of trigger patterns for all events for the specified trigger - type. + Returns an array of trigger patterns for all events for the specified trigger type. Args: - trigger (EventType): The trigger type for which the trigger patterns for all - events are requested. + trigger (EventType): The trigger type for which the trigger patterns for all events are requested. Returns: - np.ndarray: An array of trigger patterns for all events for the specified - trigger type. + np.ndarray: An array of trigger patterns for all events for the specified trigger type. """ return np.array( self.__trig_patter_all[f"{__class__._get_name_trigger(trigger)}"], - dtype=bool, + dtype=ArrayDataContainer.fields["trig_pattern_all"].dtype, ) diff --git a/src/nectarchain/makers/component/gainComponent.py b/src/nectarchain/makers/component/gainComponent.py new file mode 100644 index 00000000..60a7c1fd --- /dev/null +++ b/src/nectarchain/makers/component/gainComponent.py @@ -0,0 +1,17 @@ +import logging + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + +from abc import abstractmethod + +from .core import NectarCAMComponent + +__all__ = ["GainNectarCAMComponent"] + + +class GainNectarCAMComponent(NectarCAMComponent): + @abstractmethod + def finish(self): + pass diff --git a/src/nectarchain/makers/component/photostatistic_algorithm.py b/src/nectarchain/makers/component/photostatistic_algorithm.py new file mode 100644 index 00000000..2b622aa1 --- /dev/null +++ b/src/nectarchain/makers/component/photostatistic_algorithm.py @@ -0,0 +1,382 @@ +import logging +import sys + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + + +import copy +import os + +import matplotlib +import numpy as np +from astropy.visualization import quantity_support +from ctapipe.core import Component +from matplotlib import pyplot as plt +from scipy.stats import linregress + +from ...data.container import ChargesContainer, GainContainer, SPEfitContainer +from ..component import ChargesComponent + +__all__ = ["PhotoStatisticAlgorithm"] + + +class PhotoStatisticAlgorithm(Component): + def __init__( + self, + pixels_id: np.ndarray, + FFcharge_hg: np.ndarray, + FFcharge_lg: np.ndarray, + Pedcharge_hg: np.ndarray, + Pedcharge_lg: np.ndarray, + coefCharge_FF_Ped: float, + SPE_resolution: np.ndarray, + SPE_high_gain: np.ndarray, + config=None, + parent=None, + **kwargs, + ) -> None: + # constructors + super().__init__(config=config, parent=parent, **kwargs) + + self._pixels_id = pixels_id + self.__coefCharge_FF_Ped = coefCharge_FF_Ped + + self.__SPE_resolution = SPE_resolution + self.__SPE_high_gain = SPE_high_gain + + self.__FFcharge_hg = FFcharge_hg + self.__FFcharge_lg = FFcharge_lg + + self.__Pedcharge_hg = Pedcharge_hg + self.__Pedcharge_lg = Pedcharge_lg + + self.__check_shape() + + self.__results = GainContainer( + is_valid=np.zeros((self.npixels), dtype=bool), + high_gain=np.zeros((self.npixels, 3)), + low_gain=np.zeros((self.npixels, 3)), + pixels_id=self._pixels_id, + ) + + @classmethod + def create_from_chargesContainer( + cls, + FFcharge: ChargesContainer, + Pedcharge: ChargesContainer, + SPE_result: SPEfitContainer, + coefCharge_FF_Ped: float, + **kwargs, + ): + kwargs_init = __class__.__get_charges_FF_Ped_reshaped( + FFcharge, Pedcharge, SPE_result + ) + + kwargs.update(kwargs_init) + return cls(coefCharge_FF_Ped=coefCharge_FF_Ped, **kwargs) + + @staticmethod + def __get_charges_FF_Ped_reshaped( + FFcharge: ChargesContainer, + Pedcharge: ChargesContainer, + SPE_result: SPEfitContainer, + ) -> dict: + log.info("reshape of SPE, Ped and FF data with intersection of pixel ids") + out = {} + + FFped_intersection = np.intersect1d(Pedcharge.pixels_id, FFcharge.pixels_id) + SPEFFPed_intersection = np.intersect1d( + FFped_intersection, SPE_result.pixels_id[SPE_result.is_valid] + ) + mask_SPE = np.array( + [ + SPE_result.pixels_id[i] in SPEFFPed_intersection + for i in range(len(SPE_result.pixels_id)) + ], + dtype=bool, + ) + out["SPE_resolution"] = SPE_result.resolution[mask_SPE].T[0] + out["SPE_high_gain"] = SPE_result.high_gain[mask_SPE].T[0] + + out["pixels_id"] = SPEFFPed_intersection + out["FFcharge_hg"] = ChargesComponent.select_charges_hg( + FFcharge, SPEFFPed_intersection + ) + out["FFcharge_lg"] = ChargesComponent.select_charges_lg( + FFcharge, SPEFFPed_intersection + ) + out["Pedcharge_hg"] = ChargesComponent.select_charges_hg( + Pedcharge, SPEFFPed_intersection + ) + out["Pedcharge_lg"] = ChargesComponent.select_charges_lg( + Pedcharge, SPEFFPed_intersection + ) + + log.info(f"data have {len(SPEFFPed_intersection)} pixels in common") + return out + + def __check_shape(self) -> None: + """ + Checks the shape of certain attributes and raises an exception if the shape is not as expected. + """ + try: + self.__FFcharge_hg[0] * self.__FFcharge_lg[0] * self.__Pedcharge_hg[ + 0 + ] * self.__Pedcharge_lg[0] * self.__SPE_resolution * self._pixels_id + except Exception as e: + log.error(e, exc_info=True) + raise e + + def run(self, pixels_id: np.ndarray = None, **kwargs) -> None: + log.info("running photo statistic method") + if pixels_id is None: + pixels_id = self._pixels_id + mask = np.array( + [pixel_id in pixels_id for pixel_id in self._pixels_id], dtype=bool + ) + self._results.high_gain = np.array( + (self.gainHG * mask, np.zeros(self.npixels), np.zeros(self.npixels)) + ).T + self._results.low_gain = np.array( + (self.gainLG * mask, np.zeros(self.npixels), np.zeros(self.npixels)) + ).T + self._results.is_valid = mask + + figpath = kwargs.get("figpath", False) + if figpath: + os.makedirs(figpath, exist_ok=True) + fig = __class__.plot_correlation( + self._results.high_gain.T[0], self.__SPE_high_gain + ) + fig.savefig(f"{figpath}/plot_correlation_Photo_Stat_SPE.pdf") + fig.clf() + plt.close(fig) + return 0 + + @staticmethod + def plot_correlation( + photoStat_gain: np.ndarray, SPE_gain: np.ndarray + ) -> plt.Figure: + """ + Plot the correlation between the photo statistic gain and the single photoelectron (SPE) gain. + + Args: + photoStat_gain (np.ndarray): Array of photo statistic gain values. + SPE_gain (np.ndarray): Array of SPE gain values. + + Returns: + fig (plt.Figure): The figure object containing the scatter plot and the linear fit line. + """ + matplotlib.use("TkAgg") + # Create a mask to filter the data points based on certain criteria + mask = (photoStat_gain > 20) * (SPE_gain > 0) * (photoStat_gain < 80) + + if not (np.max(mask)): + log.debug("mask conditions are much strict, remove the mask") + mask = np.ones(len(mask), dtype=bool) + # Perform a linear regression analysis on the filtered data points + a, b, r, p_value, std_err = linregress( + photoStat_gain[mask], SPE_gain[mask], "greater" + ) + + # Generate a range of x-values for the linear fit line + x = np.linspace(photoStat_gain[mask].min(), photoStat_gain[mask].max(), 1000) + + # Define a lambda function for the linear fit line + y = lambda x: a * x + b + + with quantity_support(): + # Create a scatter plot of the filtered data points + fig, ax = plt.subplots(1, 1, figsize=(8, 6)) + ax.scatter(photoStat_gain[mask], SPE_gain[mask], marker=".") + + # Plot the linear fit line using the x-values and the lambda function + ax.plot( + x, + y(x), + color="red", + label=f"linear fit,\n a = {a:.2e},\n b = {b:.2e},\n r = {r:.2e},\n p_value = {p_value:.2e},\n std_err = {std_err:.2e}", + ) + + # Plot the line y = x + ax.plot(x, x, color="black", label="y = x") + + ax.set_xlabel("Gain Photo stat (ADC)", size=15) + ax.set_ylabel("Gain SPE fit (ADC)", size=15) + ax.legend(fontsize=15) + + return fig + + @property + def SPE_resolution(self) -> float: + """ + Returns a deep copy of the SPE resolution. + + Returns: + float: The SPE resolution. + """ + return copy.deepcopy(self.__SPE_resolution) + + @property + def sigmaPedHG(self) -> float: + """ + Calculates and returns the standard deviation of Pedcharge_hg multiplied by the square root of coefCharge_FF_Ped. + + Returns: + float: The standard deviation of Pedcharge_hg. + """ + return np.std(self.__Pedcharge_hg, axis=0) * np.sqrt(self.__coefCharge_FF_Ped) + + @property + def sigmaChargeHG(self) -> float: + """ + Calculates and returns the standard deviation of FFcharge_hg minus meanPedHG. + + Returns: + float: The standard deviation of FFcharge_hg minus meanPedHG. + """ + return np.std(self.__FFcharge_hg - self.meanPedHG, axis=0) + + @property + def meanPedHG(self) -> float: + """ + Calculates and returns the mean of Pedcharge_hg multiplied by coefCharge_FF_Ped. + + Returns: + float: The mean of Pedcharge_hg. + """ + return np.mean(self.__Pedcharge_hg, axis=0) * self.__coefCharge_FF_Ped + + @property + def meanChargeHG(self) -> float: + """ + Calculates and returns the mean of FFcharge_hg minus meanPedHG. + + Returns: + float: The mean of FFcharge_hg minus meanPedHG. + """ + return np.mean(self.__FFcharge_hg - self.meanPedHG, axis=0) + + @property + def BHG(self) -> float: + """ + Calculates and returns the BHG value. + + Returns: + float: The BHG value. + """ + min_events = np.min((self.__FFcharge_hg.shape[0], self.__Pedcharge_hg.shape[0])) + upper = ( + np.power( + self.__FFcharge_hg.mean(axis=1)[:min_events] + - self.__Pedcharge_hg.mean(axis=1)[:min_events] + * self.__coefCharge_FF_Ped + - self.meanChargeHG.mean(), + 2, + ) + ).mean(axis=0) + lower = np.power(self.meanChargeHG.mean(), 2) + return np.sqrt(upper / lower) + + @property + def gainHG(self) -> float: + """ + Calculates and returns the gain for high gain charge data. + + Returns: + float: The gain for high gain charge data. + """ + return ( + np.power(self.sigmaChargeHG, 2) + - np.power(self.sigmaPedHG, 2) + - np.power(self.BHG * self.meanChargeHG, 2) + ) / (self.meanChargeHG * (1 + np.power(self.SPE_resolution, 2))) + + @property + def sigmaPedLG(self) -> float: + """ + Calculates and returns the standard deviation of Pedcharge_lg multiplied by the square root of coefCharge_FF_Ped. + + Returns: + float: The standard deviation of Pedcharge_lg. + """ + return np.std(self.__Pedcharge_lg, axis=0) * np.sqrt(self.__coefCharge_FF_Ped) + + @property + def sigmaChargeLG(self) -> float: + """ + Calculates and returns the standard deviation of FFcharge_lg minus meanPedLG. + + Returns: + float: The standard deviation of FFcharge_lg minus meanPedLG. + """ + return np.std(self.__FFcharge_lg - self.meanPedLG, axis=0) + + @property + def meanPedLG(self) -> float: + """ + Calculates and returns the mean of Pedcharge_lg multiplied by coefCharge_FF_Ped. + + Returns: + float: The mean of Pedcharge_lg. + """ + return np.mean(self.__Pedcharge_lg, axis=0) * self.__coefCharge_FF_Ped + + @property + def meanChargeLG(self) -> float: + """ + Calculates and returns the mean of FFcharge_lg minus meanPedLG. + + Returns: + float: The mean of FFcharge_lg minus meanPedLG. + """ + return np.mean(self.__FFcharge_lg - self.meanPedLG, axis=0) + + @property + def BLG(self) -> float: + """ + Calculates and returns the BLG value. + + Returns: + float: The BLG value. + """ + min_events = np.min((self.__FFcharge_lg.shape[0], self.__Pedcharge_lg.shape[0])) + upper = ( + np.power( + self.__FFcharge_lg.mean(axis=1)[:min_events] + - self.__Pedcharge_lg.mean(axis=1)[:min_events] + * self.__coefCharge_FF_Ped + - self.meanChargeLG.mean(), + 2, + ) + ).mean(axis=0) + lower = np.power(self.meanChargeLG.mean(), 2) + return np.sqrt(upper / lower) + + @property + def gainLG(self) -> float: + """ + Calculates and returns the gain for low gain charge data. + + Returns: + float: The gain for low gain charge data. + """ + return ( + np.power(self.sigmaChargeLG, 2) + - np.power(self.sigmaPedLG, 2) + - np.power(self.BLG * self.meanChargeLG, 2) + ) / (self.meanChargeLG * (1 + np.power(self.SPE_resolution, 2))) + + @property + def results(self): + return copy.deepcopy(self.__results) + + @property + def _results(self): + return self.__results + + @property + def npixels(self): + return len(self._pixels_id) diff --git a/src/nectarchain/makers/component/photostatistic_component.py b/src/nectarchain/makers/component/photostatistic_component.py new file mode 100644 index 00000000..2a1d6e64 --- /dev/null +++ b/src/nectarchain/makers/component/photostatistic_component.py @@ -0,0 +1,134 @@ +import logging + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + +import copy + +from ctapipe.containers import EventType +from ctapipe.core.traits import List, Path, Unicode +from ctapipe_io_nectarcam.constants import N_SAMPLES +from ctapipe_io_nectarcam.containers import NectarCAMDataContainer + +from ...data.container import SPEfitContainer, merge_map_ArrayDataContainer +from ...utils import ComponentUtils +from .chargesComponent import ChargesComponent +from .gainComponent import GainNectarCAMComponent +from .photostatistic_algorithm import PhotoStatisticAlgorithm + +__all__ = ["PhotoStatisticNectarCAMComponent"] + + +class PhotoStatisticNectarCAMComponent(GainNectarCAMComponent): + SPE_result = Path( + help="the path of the SPE result container computed with very high voltage data", + ).tag(config=True) + PhotoStatAlgorithm = Unicode( + "PhotoStatisticAlgorithm", + help="The photo-statitic algorithm to be used", + read_only=True, + ).tag(config=True) + + SubComponents = copy.deepcopy(GainNectarCAMComponent.SubComponents) + SubComponents.default_value = [ + "ChargesComponent", + f"{PhotoStatAlgorithm.default_value}", + ] + SubComponents.read_only = True + + asked_pixels_id = List( + default_value=None, + allow_none=True, + help="The pixels id where we want to perform the SPE fit", + ).tag(config=True) + + # constructor + def __init__(self, subarray, config=None, parent=None, *args, **kwargs) -> None: + chargesComponent_kwargs = {} + self._PhotoStatAlgorithm_kwargs = {} + other_kwargs = {} + chargesComponent_configurable_traits = ComponentUtils.get_configurable_traits( + ChargesComponent + ) + PhotoStatAlgorithm_configurable_traits = ComponentUtils.get_configurable_traits( + eval(self.PhotoStatAlgorithm) + ) + + for key in kwargs.keys(): + if key in chargesComponent_configurable_traits.keys(): + chargesComponent_kwargs[key] = kwargs[key] + elif key in PhotoStatAlgorithm_configurable_traits.keys(): + self._PhotoStatAlgorithm_kwargs[key] = kwargs[key] + else: + other_kwargs[key] = kwargs[key] + + super().__init__( + subarray=subarray, config=config, parent=parent, *args, **other_kwargs + ) + self.FF_chargesComponent = ChargesComponent( + subarray=subarray, + config=config, + parent=parent, + *args, + **chargesComponent_kwargs, + ) + self._FF_chargesContainers = None + self.Ped_chargesComponent = ChargesComponent( + subarray=subarray, + config=config, + parent=parent, + *args, + **chargesComponent_kwargs, + ) + self._Ped_chargesContainers = None + + self.__coefCharge_FF_Ped = ( + int( + chargesComponent_kwargs.get("extractor_kwargs").get( + "window_width", N_SAMPLES + ) + ) + / N_SAMPLES + ) + + def __call__(self, event: NectarCAMDataContainer, *args, **kwargs): + if event.trigger.event_type == EventType.FLATFIELD: + self.FF_chargesComponent(event=event, *args, **kwargs) + elif event.trigger.event_type in [ + EventType.SKY_PEDESTAL, + EventType.DARK_PEDESTAL, + EventType.ELECTRONIC_PEDESTAL, + ]: + self.Ped_chargesComponent(event=event, *args, **kwargs) + else: + self.log.warning( + f"event {event.index.event_id} is event type {event.trigger.event_type} which is not used here" + ) + + def finish(self, *args, **kwargs): + if self._FF_chargesContainers is None: + self._FF_chargesContainers = self.FF_chargesComponent.finish( + *args, **kwargs + ) + self._FF_chargesContainers = merge_map_ArrayDataContainer( + self._FF_chargesContainers + ) + if self._Ped_chargesContainers is None: + self._Ped_chargesContainers = self.Ped_chargesComponent.finish( + *args, **kwargs + ) + self._Ped_chargesContainers = merge_map_ArrayDataContainer( + self._Ped_chargesContainers + ) + nectarGainSPEresult = SPEfitContainer.from_hdf5(self.SPE_result) + photo_stat = eval(self.PhotoStatAlgorithm).create_from_chargesContainer( + FFcharge=self._FF_chargesContainers, + Pedcharge=self._Ped_chargesContainers, + SPE_result=nectarGainSPEresult, + coefCharge_FF_Ped=self.__coefCharge_FF_Ped, + parent=self, + **self._PhotoStatAlgorithm_kwargs, + ) + fit_output = photo_stat.run(pixels_id=self.asked_pixels_id, *args, **kwargs) + return photo_stat.results diff --git a/src/nectarchain/makers/component/spe/__init__.py b/src/nectarchain/makers/component/spe/__init__.py new file mode 100644 index 00000000..4933f073 --- /dev/null +++ b/src/nectarchain/makers/component/spe/__init__.py @@ -0,0 +1,2 @@ +# from .parameters import * +from .spe_algorithm import * diff --git a/src/nectarchain/makers/calibration/gain/parameters.py b/src/nectarchain/makers/component/spe/parameters.py similarity index 100% rename from src/nectarchain/makers/calibration/gain/parameters.py rename to src/nectarchain/makers/component/spe/parameters.py diff --git a/src/nectarchain/makers/calibration/gain/parameters_signal_fromHHVFit.yaml b/src/nectarchain/makers/component/spe/parameters_SPECombined_fromHHVFit.yaml similarity index 100% rename from src/nectarchain/makers/calibration/gain/parameters_signal_fromHHVFit.yaml rename to src/nectarchain/makers/component/spe/parameters_SPECombined_fromHHVFit.yaml diff --git a/src/nectarchain/makers/calibration/gain/parameters_signal.yaml b/src/nectarchain/makers/component/spe/parameters_SPEHHV.yaml similarity index 100% rename from src/nectarchain/makers/calibration/gain/parameters_signal.yaml rename to src/nectarchain/makers/component/spe/parameters_SPEHHV.yaml diff --git a/src/nectarchain/makers/calibration/gain/parameters_signalStd.yaml b/src/nectarchain/makers/component/spe/parameters_SPEHHVStd.yaml similarity index 100% rename from src/nectarchain/makers/calibration/gain/parameters_signalStd.yaml rename to src/nectarchain/makers/component/spe/parameters_SPEHHVStd.yaml diff --git a/src/nectarchain/makers/component/spe/parameters_SPEnominal.yaml b/src/nectarchain/makers/component/spe/parameters_SPEnominal.yaml new file mode 100644 index 00000000..18b6ee31 --- /dev/null +++ b/src/nectarchain/makers/component/spe/parameters_SPEnominal.yaml @@ -0,0 +1,32 @@ +{ + pedestalWidth : { + value : 50, + min : 1, + max : 150 + }, + luminosity : { + value : 1, + min : 0.01, + max : 5.0 + }, + pp : { + value : 0.45, + min : 0.2, + max : 0.8 + }, + resolution : { + value : 0.5, + min : 0.3, + max : 0.7 + }, + n : { + value : 0.697, + min : 0.5, + max : 0.9 + }, + mean : { + value : 50, + min : 10, + max : 200 + } +} \ No newline at end of file diff --git a/src/nectarchain/makers/component/spe/parameters_SPEnominalStd.yaml b/src/nectarchain/makers/component/spe/parameters_SPEnominalStd.yaml new file mode 100644 index 00000000..3be2f77f --- /dev/null +++ b/src/nectarchain/makers/component/spe/parameters_SPEnominalStd.yaml @@ -0,0 +1,32 @@ +{ + pedestalWidth : { + value : 50, + min : 1, + max : 150 + }, + luminosity : { + value : 1, + min : 0.01, + max : 5.0 + }, + pp : { + value : 0.45, + min : .NAN, + max : .NAN + }, + resolution : { + value : 0.5, + min : 0.3, + max : 0.7 + }, + n : { + value : 0.697, + min : .NAN, + max : .NAN + }, + mean : { + value : 50, + min : 10, + max : 200 + } +} \ No newline at end of file diff --git a/src/nectarchain/makers/calibration/gain/parameters_signal_combined.yaml b/src/nectarchain/makers/component/spe/parameters_signal_combined.yaml similarity index 100% rename from src/nectarchain/makers/calibration/gain/parameters_signal_combined.yaml rename to src/nectarchain/makers/component/spe/parameters_signal_combined.yaml diff --git a/src/nectarchain/makers/calibration/gain/flatfield_spe_makers.py b/src/nectarchain/makers/component/spe/spe_algorithm.py similarity index 58% rename from src/nectarchain/makers/calibration/gain/flatfield_spe_makers.py rename to src/nectarchain/makers/component/spe/spe_algorithm.py index 74abd003..6c656a96 100644 --- a/src/nectarchain/makers/calibration/gain/flatfield_spe_makers.py +++ b/src/nectarchain/makers/component/spe/spe_algorithm.py @@ -1,5 +1,10 @@ -import copy import logging + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + +import copy import os import time from inspect import signature @@ -7,10 +12,17 @@ from typing import Tuple import astropy.units as u +import matplotlib import matplotlib.pyplot as plt +import matplotlib.style as mplstyle + +mplstyle.use("fast") + import numpy as np import yaml -from astropy.table import Column, QTable +from astropy.table import QTable +from ctapipe.core.component import Component +from ctapipe.core.traits import Bool, Float, Integer, Path, Unicode from iminuit import Minuit from matplotlib.colors import to_rgba from matplotlib.patches import Rectangle @@ -18,133 +30,96 @@ from scipy.signal import find_peaks, savgol_filter from scipy.special import gammainc -from ....data.container import ChargesContainer -from ...charges_makers import ChargesMaker -from .gain_makers import GainMaker +from ....data.container import ChargesContainer, SPEfitContainer +from ....utils import MPE2, MeanValueError, Statistics, UtilsMinuit, weight_gaussian +from ..chargesComponent import ChargesComponent from .parameters import Parameter, Parameters -from .utils import MPE2, MeanValueError, Statistics, UtilsMinuit, weight_gaussian -logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") -log = logging.getLogger(__name__) -log.handlers = logging.getLogger("__main__").handlers +__all__ = [ + "SPEHHValgorithm", + "SPEHHVStdalgorithm", + "SPEnominalStdalgorithm", + "SPEnominalalgorithm", + "SPECombinedalgorithm", +] -__all__ = ["FlatFieldSingleHHVSPEMaker", "FlatFieldSingleHHVStdSPEMaker"] - - -class FlatFieldSPEMaker(GainMaker): - - """ - The `FlatFieldSPEMaker` class is used for flat field single photoelectron (SPE) - calibration calculations on data. It inherits from the `GainMaker` class and adds - functionality specific to flat field SPE calibration. - - Example Usage: - # Create an instance of the FlatFieldSPEMaker class - flat_field_maker = FlatFieldSPEMaker() - - # Read parameters from a YAML file - flat_field_maker.read_param_from_yaml("parameters.yaml") - - # Update parameters based on data - flat_field_maker._update_parameters(parameters, charge, counts) - - # Update the result table based on the parameters - flat_field_maker._update_table_from_parameters() - - Main functionalities: - - Inherits from the `GainMaker` class and adds functionality specific to flat - field SPE calibration. - - Reads parameters from a YAML file and updates the internal parameters of the - class. - - Updates the parameters based on data, such as charge and counts. - - Updates a result table based on the parameters. - - Methods: - - `read_param_from_yaml(parameters_file, only_update)`: Reads parameters from a - YAML file and updates the internal parameters of the class. If `only_update` is - True, only the parameters that exist in the YAML file will be updated. - - `_update_parameters(parameters, charge, counts, **kwargs)`: Updates the - parameters based on data, such as charge and counts. It performs a Gaussian fit - on the data to determine the pedestal and mean values, and updates the - corresponding parameters accordingly. - - `_get_mean_gaussian_fit(charge, counts, extension, **kwargs)`: Performs a - Gaussian fit on the data to determine the pedestal and mean values. It returns - the fit coefficients. - - `_update_table_from_parameters()`: Updates a result table based on the parameters. - It adds columns to the table for each - parameter and its corresponding error. - - Attributes: - - `_Windows_length`: A class attribute that represents the length of the windows - used for smoothing the data. - - `_Order`: A class attribute that represents the order of the polynomial used - for smoothing the data. - - Members: - - `npixels`: A property that returns the number of pixels. - - `parameters`: A property that returns a deep copy of the internal parameters - of the class. - - `_parameters`: A property that returns the internal parameters of the class. - """ - - _Windows_length = 40 - _Order = 2 - - # constructors - def __init__(self, *args, **kwargs) -> None: - """ - Initializes the FlatFieldSPEMaker class. - Parameters: - *args: Variable length argument list. - **kwargs: Arbitrary keyword arguments. +class SPEalgorithm(Component): + Windows_lenght = Integer( + 40, + read_only=True, + help="The windows leght used for the savgol filter algorithm", + ).tag(config=True) - Returns: - None - """ - super().__init__(*args, **kwargs) + Order = Integer( + 2, + read_only=True, + help="The order of the polynome used in the savgol filter algorithm", + ).tag(config=True) + + def __init__( + self, pixels_id: np.ndarray, config=None, parent=None, **kwargs + ) -> None: + super().__init__(config=config, parent=parent, **kwargs) + self.__pixels_id = pixels_id + self.__pedestal = Parameter( + name="pedestal", + value=0, + ) self.__parameters = Parameters() + self.__parameters.append(self.__pedestal) + self.__results = SPEfitContainer( + is_valid=np.zeros((self.npixels), dtype=bool), + high_gain=np.empty((self.npixels, 3)), + low_gain=np.empty((self.npixels, 3)), + pixels_id=self.__pixels_id, + likelihood=np.empty((self.npixels)), + p_value=np.empty((self.npixels)), + pedestal=np.empty((self.npixels, 3)), + pedestalWidth=np.empty((self.npixels, 3)), + resolution=np.empty((self.npixels, 3)), + luminosity=np.empty((self.npixels, 3)), + mean=np.empty((self.npixels, 3)), + n=np.empty((self.npixels, 3)), + pp=np.empty((self.npixels, 3)), + ) @property - def npixels(self): - """ - Returns the number of pixels. + def parameters(self): + return copy.deepcopy(self.__parameters) - Returns: - int: The number of pixels. - """ - return len(self._pixels_id) + @property + def _parameters(self): + return self.__parameters @property - def parameters(self): - """ - Returns a deep copy of the internal parameters. + def results(self): + return copy.deepcopy(self.__results) - Returns: - dict: A deep copy of the internal parameters. - """ - return copy.deepcopy(self.__parameters) + @property + def _results(self): + return self.__results @property - def _parameters(self): - """ - Returns the internal parameters. + def pixels_id(self): + return copy.deepcopy(self.__pixels_id) - Returns: - dict: The internal parameters. - """ - return self.__parameters + @property + def _pixels_id(self): + return self.__pixels_id + + @property + def npixels(self): + return len(self.__pixels_id) # methods def read_param_from_yaml(self, parameters_file, only_update=False) -> None: """ - Reads parameters from a YAML file and updates the internal parameters of the - FlatFieldSPEMaker class. + Reads parameters from a YAML file and updates the internal parameters of the FlatFieldSPEMaker class. Args: parameters_file (str): The name of the YAML file containing the parameters. - only_update (bool, optional): If True, only the parameters that exist in - the YAML file will be updated. Default is False. + only_update (bool, optional): If True, only the parameters that exist in the YAML file will be updated. Default is False. Returns: None @@ -180,19 +155,16 @@ def _update_parameters( parameters: Parameters, charge: np.ndarray, counts: np.ndarray, **kwargs ) -> Parameters: """ - Update the parameters of the FlatFieldSPEMaker class based on the input - charge and counts data. + Update the parameters of the FlatFieldSPEMaker class based on the input charge and counts data. Args: - parameters (Parameters): An instance of the Parameters class that holds - the internal parameters of the FlatFieldSPEMaker class. + parameters (Parameters): An instance of the Parameters class that holds the internal parameters of the FlatFieldSPEMaker class. charge (np.ndarray): An array of charge values. counts (np.ndarray): An array of corresponding counts values. **kwargs: Additional keyword arguments. Returns: - Parameters: The updated parameters object with the pedestal and mean - values and their corresponding limits. + Parameters: The updated parameters object with the pedestal and mean values and their corresponding limits. """ try: coeff_ped, coeff_mean = __class__._get_mean_gaussian_fit( @@ -229,7 +201,7 @@ def _update_parameters( @staticmethod def _get_mean_gaussian_fit( - charge: np.ndarray, counts: np.ndarray, extension: str = "", **kwargs + charge: np.ndarray, counts: np.ndarray, pixel_id=None, **kwargs ) -> Tuple[np.ndarray, np.ndarray]: """ Perform a Gaussian fit on the data to determine the pedestal and mean values. @@ -237,12 +209,11 @@ def _get_mean_gaussian_fit( Args: charge (np.ndarray): An array of charge values. counts (np.ndarray): An array of corresponding counts. - extension (str, optional): An extension string. Defaults to "". + pixel_id (int): The id of the current pixel. Default to None **kwargs: Additional keyword arguments. Returns: - Tuple[np.ndarray, np.ndarray]: A tuple of fit coefficients for the - pedestal and mean. + Tuple[np.ndarray, np.ndarray]: A tuple of fit coefficients for the pedestal and mean. Example Usage: flat_field_maker = FlatFieldSPEMaker() @@ -252,9 +223,9 @@ def _get_mean_gaussian_fit( print(coeff) # Output: [norm,peak_value, peak_width] print(coeff_mean) # Output: [norm,peak_value_mean, peak_width_mean] """ - windows_length = __class__._Windows_length - order = __class__._Order - histo_smoothed = savgol_filter(counts, windows_length, order) + windows_lenght = __class__.Windows_lenght.default_value + order = __class__.Order.default_value + histo_smoothed = savgol_filter(counts, windows_lenght, order) peaks = find_peaks(histo_smoothed, 10) peak_max = np.argmax(histo_smoothed[peaks[0]]) peak_pos, peak_value = charge[peaks[0][peak_max]], counts[peaks[0][peak_max]] @@ -273,8 +244,7 @@ def _get_mean_gaussian_fit( ax.plot( charge, histo_smoothed, - label=f"smoothed data with savgol filter (windows length :" - f" {windows_length}, order : {order})", + label=f"smoothed data with savgol filter (windows lenght : {windows_lenght}, order : {order})", ) ax.plot( charge, @@ -305,9 +275,7 @@ def _get_mean_gaussian_fit( exist_ok=True, ) fig.savefig( - f"{os.environ.get('NECTARCHAIN_LOG')}/" - f"{os.getpid()}/figures/initialization_pedestal_pixel{extension}_" - f"{os.getpid()}.pdf" + f"{os.environ.get('NECTARCHAIN_LOG')}/{os.getpid()}/figures/initialization_pedestal_pixel{pixel_id}_{os.getpid()}.pdf" ) fig.clf() plt.close(fig) @@ -337,8 +305,7 @@ def _get_mean_gaussian_fit( ax.plot( charge, histo_smoothed, - label=f"smoothed data with savgol filter (windows length :" - f" {windows_length}, order : {order})", + label=f"smoothed data with savgol filter (windows lenght : {windows_lenght}, order : {order})", ) ax.plot( charge, @@ -365,28 +332,27 @@ def _get_mean_gaussian_fit( ax.set_ylabel("Events", size=15) ax.legend(fontsize=7) os.makedirs( - f"{os.environ.get('NECTARCHAIN_LOG')}/{os.getpid()}/figures/", + f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{os.getpid()}/figures/", exist_ok=True, ) fig.savefig( - f"{os.environ.get('NECTARCHAIN_LOG')}/" - f"{os.getpid()}/figures/initialization_mean_pixel{extension}_" - f"{os.getpid()}.pdf" + f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{os.getpid()}/figures/initialization_mean_pixel{pixel_id}_{os.getpid()}.pdf" ) fig.clf() plt.close(fig) del fig, ax return coeff, coeff_mean + +''' def _update_table_from_parameters(self) -> None: """ Update the result table based on the parameters of the FlatFieldSPEMaker class. - This method adds columns to the table for each parameter and its - corresponding error. + This method adds columns to the table for each parameter and its corresponding error. """ for param in self._parameters.parameters: - if not (param.name in self._results.colnames): + if not (param.name in self._results.keys()): self._results.add_column( Column( data=np.empty((self.npixels), dtype=np.float64), @@ -401,44 +367,48 @@ def _update_table_from_parameters(self) -> None: unit=param.unit, ) ) +''' -class FlatFieldSingleHHVSPEMaker(FlatFieldSPEMaker): - """ - This class represents a FlatFieldSingleHHVSPEMaker object. +class SPEnominalalgorithm(SPEalgorithm): + parameters_file = Unicode( + "parameters_SPEnominal.yaml", + read_only=True, + help="The name of the SPE fit parameters file", + ).tag(config=True) - Args: - charge (np.ma.masked_array or array-like): The charge data. - counts (np.ma.masked_array or array-like): The counts data. - *args: Additional positional arguments. - **kwargs: Additional keyword arguments. - Attributes: - __charge (np.ma.masked_array): The charge data as a masked array. - __counts (np.ma.masked_array): The counts data as a masked array. - __pedestal (Parameter): The pedestal value. - _parameters (list): List of parameters. - __parameters_file (str): The path to the parameters file. - _results (Table): Table of results. - Methods: - __init__: Initializes the FlatFieldSingleHHVSPEMaker object. - create_from_chargesContainer: Creates an instance of FlatFieldSingleHHVSPEMaker - using charge and counts data from a ChargesContainer object. - create_from_run_number(cls, run_number, **kwargs): Class method that creates an - instance from a run number. - make(self, pixels_id=None, multiproc=True, display=True, **kwargs): Method that - performs the fit on the specified pixels and returns the fit results. - display(self, pixels_id, **kwargs): Method that plots the fit for the specified - pixels. - """ - - __parameters_file = "parameters_signal.yaml" __fit_array = None - _reduced_name = "FlatFieldSingleSPE" - __nproc_default = 8 - __chunksize_default = 1 - # constructors - def __init__(self, charge, counts, *args, **kwargs) -> None: + tol = Float( + 1e-1, + help="The tolerance used for minuit", + read_only=True, + ).tag(config=True) + + nproc = Integer( + 8, + help="The Number of cpu used for SPE fit", + ).tag(config=True) + + chunksize = Integer( + 1, + help="The chunk size for multi-processing", + ).tag(config=True) + + multiproc = Bool( + True, + help="flag to active multi-processing", + ).tag(config=True) + + def __init__( + self, + pixels_id: np.ndarray, + charge: np.ndarray, + counts: np.ndarray, + config=None, + parent=None, + **kwargs, + ) -> None: """ Initializes the FlatFieldSingleHHVSPEMaker object. Args: @@ -447,7 +417,7 @@ def __init__(self, charge, counts, *args, **kwargs) -> None: *args: Additional positional arguments. **kwargs: Additional keyword arguments. """ - super().__init__(*args, **kwargs) + super().__init__(pixels_id=pixels_id, config=config, parent=parent, **kwargs) if isinstance(charge, np.ma.masked_array): self.__charge = charge else: @@ -456,56 +426,29 @@ def __init__(self, charge, counts, *args, **kwargs) -> None: self.__counts = counts else: self.__counts = np.ma.asarray(counts) - self.__pedestal = Parameter( - name="pedestal", - value=( - np.min(self.__charge) - + np.sum(self.__charge * self.__counts) / np.sum(self.__counts) - ) - / 2, - min=np.min(self.__charge), - max=np.sum(self.__charge * self.__counts) / np.sum(self.__counts), - unit=u.dimensionless_unscaled, - ) - self._parameters.append(self.__pedestal) - self.read_param_from_yaml(kwargs.get("parameters_file", self.__parameters_file)) - self._update_table_from_parameters() - self._results.add_column( - Column( - np.zeros((self.npixels), dtype=np.float64), - "likelihood", - unit=u.dimensionless_unscaled, - ) - ) - self._results.add_column( - Column( - np.zeros((self.npixels), dtype=np.float64), - "pvalue", - unit=u.dimensionless_unscaled, - ) - ) + + self.read_param_from_yaml(kwargs.get("parameters_file", self.parameters_file)) @classmethod - def create_from_chargesContainer(cls, signal: ChargesContainer, **kwargs): + def create_from_chargesContainer( + cls, signal: ChargesContainer, config=None, parent=None, **kwargs + ): """ - Creates an instance of FlatFieldSingleHHVSPEMaker using charge and counts - data from a ChargesContainer object. + Creates an instance of FlatFieldSingleHHVSPEMaker using charge and counts data from a ChargesContainer object. Args: signal (ChargesContainer): The ChargesContainer object. **kwargs: Additional keyword arguments. Returns: FlatFieldSingleHHVSPEMaker: An instance of FlatFieldSingleHHVSPEMaker. """ - histo = ChargesMaker.histo_hg(signal, autoscale=True) + histo = ChargesComponent.histo_hg(signal, autoscale=True) return cls( - charge=histo[1], counts=histo[0], pixels_id=signal.pixels_id, **kwargs - ) - - @classmethod - def create_from_run_number(cls, run_number: int, **kwargs): - raise NotImplementedError( - "Need to implement here the use of the WaveformsMaker and ChargesMaker to " - "produce the chargesContainer to be pass into the __ini__" + pixels_id=signal.pixels_id, + charge=histo[1], + counts=histo[0], + config=config, + parent=parent, + **kwargs, ) # getters and setters @@ -540,8 +483,7 @@ def _counts(self): # methods def _fill_results_table_from_dict(self, dico: dict, pixels_id: np.ndarray) -> None: """ - Populates the results table with fit values and errors for each pixel based - on the dictionary provided as input. + Populates the results table with fit values and errors for each pixel based on the dictionary provided as input. Args: dico (dict): A dictionary containing fit values and errors for each pixel. @@ -550,38 +492,36 @@ def _fill_results_table_from_dict(self, dico: dict, pixels_id: np.ndarray) -> No Returns: None """ + ########NEED TO BE OPTIMIZED!!!########### chi2_sig = signature(__class__.cost(self._charge, self._counts)) for i in range(len(pixels_id)): values = dico[i].get(f"values_{i}", None) errors = dico[i].get(f"errors_{i}", None) if not ((values is None) or (errors is None)): - index = np.argmax(self._results["pixels_id"] == pixels_id[i]) + index = np.argmax(self._results.pixels_id == pixels_id[i]) if len(values) != len(chi2_sig.parameters): e = Exception( - "the size out the minuit output parameters values array does " - "not fit the signature of the minimized cost function" + "the size out the minuit output parameters values array does not fit the signature of the minimized cost function" ) - log.error(e, exc_info=True) + self.log.error(e, exc_info=True) raise e for j, key in enumerate(chi2_sig.parameters): - self._results[key][index] = values[j] - self._results[f"{key}_error"][index] = errors[j] + self._results[key][index][0] = values[j] + self._results[key][index][1] = errors[j] + self._results[key][index][2] = errors[j] if key == "mean": - self._high_gain[index] = values[j] - self._results["high_gain_error"][index] = [ - errors[j], - errors[j], - ] - self._results["high_gain"][index] = values[j] - self._results["is_valid"][index] = True - self._results["likelihood"][index] = __class__.__fit_array[i].fcn( + self._results.high_gain[index][0] = values[j] + self._results.high_gain[index][1] = errors[j] + self._results.high_gain[index][2] = errors[j] + self._results.is_valid[index] = True + self._results.likelihood[index] = __class__.__fit_array[i].fcn( __class__.__fit_array[i].values ) ndof = ( self._counts.data[index][~self._counts.mask[index]].shape[0] - __class__.__fit_array[i].nfit ) - self._results["pvalue"][index] = Statistics.chi2_pvalue( + self._results.p_value[index] = Statistics.chi2_pvalue( ndof, __class__.__fit_array[i].fcn(__class__.__fit_array[i].values) ) @@ -627,8 +567,7 @@ def _NG_Likelihood_Chi2( @staticmethod def cost(charge: np.ndarray, counts: np.ndarray): """ - Defines a function called Chi2 that calculates the chi-square value using - the _NG_Likelihood_Chi2 method. + Defines a function called Chi2 that calculates the chi-square value using the _NG_Likelihood_Chi2 method. Parameters: charge (np.ndarray): An array of charge values. counts (np.ndarray): An array of count values. @@ -645,9 +584,7 @@ def Chi2( n: float, pedestalWidth: float, ): - # assert not(np.isnan(pp) or np.isnan(resolution) or np.isnan(mean) or - # np.isnan(n) or np.isnan(pedestal) or np.isnan(pedestalWidth) or - # np.isnan(luminosity)) + # assert not(np.isnan(pp) or np.isnan(resolution) or np.isnan(mean) or np.isnan(n) or np.isnan(pedestal) or np.isnan(pedestalWidth) or np.isnan(luminosity)) for i in range(1000): if gammainc(i + 1, luminosity) < 1e-5: ntotalPE = i @@ -668,17 +605,15 @@ def Chi2( return Chi2 - # @njit(parallel=True, nopython=True) + # @njit(parallel=True,nopython = True) def _make_fit_array_from_parameters( self, pixels_id: np.ndarray = None, **kwargs ) -> np.ndarray: """ - Create an array of Minuit fit instances based on the parameters and data for - each pixel. + Create an array of Minuit fit instances based on the parameters and data for each pixel. Args: - pixels_id (optional): An array of pixel IDs. If not provided, all pixels - will be used. + pixels_id (optional): An array of pixel IDs. If not provided, all pixels will be used. Returns: np.ndarray: An array of Minuit fit instances, one for each pixel. @@ -705,7 +640,7 @@ def _make_fit_array_from_parameters( parname: minuitParameters["values"][parname] for parname in minuitParameters["values"] } - log.info(f"creation of fit instance for pixel: {_id}") + self.log.info(f"creation of fit instance for pixel: {_id}") fit_array[i] = Minuit( __class__.cost( self._charge[index].data[~self._charge[index].mask], @@ -713,18 +648,18 @@ def _make_fit_array_from_parameters( ), **minuit_kwargs, ) - log.debug("fit created") + self.log.debug("fit created") fit_array[i].errordef = Minuit.LIKELIHOOD fit_array[i].strategy = 0 - fit_array[i].tol = 1e40 + fit_array[i].tol = self.tol fit_array[i].print_level = 1 fit_array[i].throw_nan = True UtilsMinuit.set_minuit_parameters_limits_and_errors( fit_array[i], minuitParameters ) - log.debug(fit_array[i].values) - log.debug(fit_array[i].limits) - log.debug(fit_array[i].fixed) + # self.log.debug(fit_array[i].values) + # self.log.debug(fit_array[i].limits) + # self.log.debug(fit_array[i].fixed) return fit_array @@ -737,10 +672,8 @@ def run_fit(i: int) -> dict: i (int): The index of the pixel to perform the fit on. Returns: - dict: A dictionary containing the fit values and errors for - the specified pixel. - The keys are "values_i" and "errors_i", where "i" is the index of - the pixel. + dict: A dictionary containing the fit values and errors for the specified pixel. + The keys are "values_i" and "errors_i", where "i" is the index of the pixel. """ log.info("Starting") __class__.__fit_array[i].migrad() @@ -750,70 +683,45 @@ def run_fit(i: int) -> dict: log.info("Finished") return {f"values_{i}": _values, f"errors_{i}": _errors} - def make( + def run( self, pixels_id: np.ndarray = None, - multiproc: bool = True, display: bool = True, **kwargs, ) -> np.ndarray: - """ - Perform a fit on specified pixels and return the fit results. - - Args: - pixels_id (np.ndarray, optional): An array of pixel IDs to perform the - fit on. If not provided, the fit will be performed on all pixels. - Default is None. - multiproc (bool, optional): A boolean indicating whether to use - multiprocessing for the fit. Default is True. - display (bool, optional): A boolean indicating whether to display - the fit results. Default is True. - **kwargs (optional): Additional keyword arguments. - - Returns: - np.ndarray: An array of fit instances. - - Example Usage: - # Initialize the FlatFieldSingleHHVSPEMaker object - maker = FlatFieldSingleHHVSPEMaker(charge, counts) - - # Perform the fit on all pixels and display the fit results - results = maker.make() - - # Perform the fit on specific pixels and display the fit results - results = maker.make(pixels_id=[1, 2, 3]) - """ - log.info("running maker") - log.info("checking asked pixels id") + self.log.info("running maker") + self.log.info("checking asked pixels id") if pixels_id is None: pixels_id = self.pixels_id npix = self.npixels else: - log.debug("checking that asked pixels id are in data") + self.log.debug("checking that asked pixels id are in data") pixels_id = np.asarray(pixels_id) mask = np.array([_id in self.pixels_id for _id in pixels_id], dtype=bool) if False in mask: - log.debug(f"The following pixels are not in data : {pixels_id[~mask]}") + self.log.debug( + f"The following pixels are not in data : {pixels_id[~mask]}" + ) pixels_id = pixels_id[mask] npix = len(pixels_id) if npix == 0: - log.warning("The asked pixels id are all out of the data") + self.log.warning("The asked pixels id are all out of the data") return None else: - log.info("creation of the fits instance array") + self.log.info("creation of the fits instance array") __class__.__fit_array = self._make_fit_array_from_parameters( pixels_id=pixels_id, display=display, **kwargs ) - log.info("running fits") - if multiproc: - nproc = kwargs.get("nproc", __class__.__nproc_default) + self.log.info("running fits") + if self.multiproc: + nproc = kwargs.get("nproc", self.nproc) chunksize = kwargs.get( "chunksize", - max(__class__.__chunksize_default, npix // (nproc * 10)), + max(self.chunksize, npix // (nproc * 10)), ) - log.info(f"pooling with nproc {nproc}, chunksize {chunksize}") + self.log.info(f"pooling with nproc {nproc}, chunksize {chunksize}") t = time.time() with Pool(nproc) as pool: @@ -826,33 +734,109 @@ def make( try: res = result.get() except Exception as e: - log.error(e, exc_info=True) + self.log.error(e, exc_info=True) raise e - log.debug(res) - log.info( - f"time for multiproc with starmap_async execution is" - f" {time.time() - t:.2e} sec" + self.log.debug(str(res)) + self.log.info( + f"time for multiproc with starmap_async execution is {time.time() - t:.2e} sec" ) else: - log.info("running in mono-cpu") + self.log.info("running in mono-cpu") t = time.time() res = [__class__.run_fit(i) for i in range(npix)] - log.debug(res) - log.info(f"time for singleproc execution is {time.time() - t:.2e} sec") + self.log.debug(res) + self.log.info( + f"time for singleproc execution is {time.time() - t:.2e} sec" + ) - log.info("filling result table from fits results") + self.log.info("filling result table from fits results") self._fill_results_table_from_dict(res, pixels_id) output = copy.copy(__class__.__fit_array) __class__.__fit_array = None if display: - log.info("plotting") + self.log.info("plotting") + t = time.time() self.display(pixels_id, **kwargs) + log.info( + f"time for plotting {len(pixels_id)} pixels : {time.time() - t:.2e} sec" + ) return output - def plot_single( + def plot_single_pyqtgraph( + pixel_id: int, + charge: np.ndarray, + counts: np.ndarray, + pp: float, + resolution: float, + gain: float, + gain_error: float, + n: float, + pedestal: float, + pedestalWidth: float, + luminosity: float, + likelihood: float, + ) -> tuple: + import pyqtgraph as pg + from pyqtgraph.Qt import QtCore, QtGui + + # from pyqtgraph.Qt import QtGui + + app = pg.mkQApp(name="minimal") + # + ## Create a window + win = pg.GraphicsLayoutWidget(show=False) + win.setWindowTitle(f"SPE fit pixel id : {pixel_id}") + + # Add a plot to the window + plot = win.addPlot(title=f"SPE fit pixel id : {pixel_id}") + + plot.addLegend() + error = pg.ErrorBarItem( + x=charge, y=counts, top=np.sqrt(counts), bottom=np.sqrt(counts), beam=0.5 + ) + plot.addItem(error) + plot.plot( + x=charge, + y=np.trapz(counts, charge) + * MPE2( + charge, + pp, + resolution, + gain, + n, + pedestal, + pedestalWidth, + luminosity, + ), + name=f"SPE model fit", + ) + legend = pg.TextItem( + f"SPE model fit gain : {gain - gain_error:.2f} < {gain:.2f} < {gain + gain_error:.2f} ADC/pe,\n likelihood : {likelihood:.2f}", + color=(200, 200, 200), + ) + legend.setPos(pedestal, np.max(counts) / 2) + font = QtGui.QFont() + font.setPointSize(12) + legend.setFont(font) + legend.setTextWidth(500) + plot.addItem(legend) + + label_style = {"color": "#EEE", "font-size": "12pt"} + plot.setLabel("bottom", "Charge (ADC)", **label_style) + plot.setLabel("left", "Events", **label_style) + plot.setRange( + xRange=[ + pedestal - 6 * pedestalWidth, + np.quantile(charge.data[~charge.mask], 0.84), + ] + ) + # ax.legend(fontsize=18) + return win + + def plot_single_matplotlib( pixel_id: int, charge: np.ndarray, counts: np.ndarray, @@ -872,8 +856,7 @@ def plot_single( Args: pixel_id (int): The ID of the pixel for which the plot is generated. charge (np.ndarray): An array of charge values. - counts (np.ndarray): An array of event counts corresponding to the - charge values. + counts (np.ndarray): An array of event counts corresponding to the charge values. pp (float): The value of the `pp` parameter. resolution (float): The value of the `resolution` parameter. gain (float): The value of the `gain` parameter. @@ -885,8 +868,7 @@ def plot_single( likelihood (float): The value of the `likelihood` parameter. Returns: - tuple: A tuple containing the generated plot figure and the axes of the - plot. + tuple: A tuple containing the generated plot figure and the axes of the plot. """ fig, ax = plt.subplots(1, 1, figsize=(8, 8)) ax.errorbar(charge, counts, np.sqrt(counts), zorder=0, fmt=".", label="data") @@ -905,57 +887,117 @@ def plot_single( ), zorder=1, linewidth=2, - label=f"SPE model fit \n gain : {gain - gain_error:.2f} < {gain:.2f} <" - f" {gain + gain_error:.2f} ADC/pe,\n likelihood : {likelihood:.2f}", + label=f"SPE model fit \n gain : {gain - gain_error:.2f} < {gain:.2f} < {gain + gain_error:.2f} ADC/pe,\n likelihood : {likelihood:.2f}", ) ax.set_xlabel("Charge (ADC)", size=15) ax.set_ylabel("Events", size=15) ax.set_title(f"SPE fit pixel id : {pixel_id}") - ax.set_xlim([pedestal - 6 * pedestalWidth, None]) + ax.set_xlim([pedestal - 6 * pedestalWidth, np.max(charge)]) ax.legend(fontsize=18) return fig, ax - def display(self, pixels_id: np.ndarray, **kwargs) -> None: + def display(self, pixels_id: np.ndarray, package="pyqtgraph", **kwargs) -> None: """ Display and save the plot for each specified pixel ID. Args: pixels_id (np.ndarray): An array of pixel IDs. + package (str): the package use to plot, can be matplotlib or pyqtgraph. Default to pyqtgraph **kwargs: Additional keyword arguments. - figpath (str): The path to save the generated plot figures. Defaults - to "/tmp/NectarGain_pid{os.getpid()}". + figpath (str): The path to save the generated plot figures. Defaults to "/tmp/NectarGain_pid{os.getpid()}". """ - figpath = kwargs.get("figpath", f"/tmp/NectarGain_pid{os.getpid()}") + figpath = kwargs.get( + "figpath", + f"{os.environ.get('NECTARCHAIN_FIGURES','/tmp')}/NectarGain_pid{os.getpid()}", + ) + self.log.debug(f"saving figures in {figpath}") os.makedirs(figpath, exist_ok=True) - for _id in pixels_id: - index = np.argmax(self._results["pixels_id"] == _id) - fig, ax = __class__.plot_single( - _id, - self._charge[index], - self._counts[index], - self._results["pp"][index].value, - self._results["resolution"][index].value, - self._results["high_gain"][index].value, - self._results["high_gain_error"][index].value.mean(), - self._results["n"][index].value, - self._results["pedestal"][index].value, - self._results["pedestalWidth"][index].value, - self._results["luminosity"][index].value, - self._results["likelihood"][index], - ) - fig.savefig(f"{figpath}/fit_SPE_pixel{_id}.pdf") - fig.clf() - plt.close(fig) - del fig, ax + if package == "matplotlib": + matplotlib.use("TkAgg") + for _id in pixels_id: + index = np.argmax(self._results.pixels_id == _id) + fig, ax = __class__.plot_single_matplotlib( + _id, + self._charge[index], + self._counts[index], + self._results.pp[index][0], + self._results.resolution[index][0], + self._results.high_gain[index][0], + self._results.high_gain[index][1:].mean(), + self._results.n[index][0], + self._results.pedestal[index][0], + self._results.pedestalWidth[index][0], + self._results.luminosity[index][0], + self._results.likelihood[index], + ) + fig.savefig(f"{figpath}/fit_SPE_pixel{_id}.pdf") + fig.clf() + plt.close(fig) + del fig, ax + elif package == "pyqtgraph": + import pyqtgraph as pg + import pyqtgraph.exporters + + for _id in pixels_id: + index = np.argmax(self._results.pixels_id == _id) + try: + widget = None + widget = __class__.plot_single_pyqtgraph( + _id, + self._charge[index], + self._counts[index], + self._results.pp[index][0], + self._results.resolution[index][0], + self._results.high_gain[index][0], + self._results.high_gain[index][1:].mean(), + self._results.n[index][0], + self._results.pedestal[index][0], + self._results.pedestalWidth[index][0], + self._results.luminosity[index][0], + self._results.likelihood[index], + ) + exporter = pg.exporters.ImageExporter(widget.getItem(0, 0)) + exporter.parameters()["width"] = 1000 + exporter.export(f"{figpath}/fit_SPE_pixel{_id}.png") + except Exception as e: + log.warning(e, exc_info=True) + finally: + del widget + + +class SPEHHValgorithm(SPEnominalalgorithm): + """class to perform fit of the SPE HHV signal with n and pp free""" + parameters_file = Unicode( + "parameters_SPEHHV.yaml", + read_only=True, + help="The name of the SPE fit parameters file", + ).tag(config=True) + tol = Float( + 1e40, + help="The tolerance used for minuit", + read_only=True, + ).tag(config=True) -class FlatFieldSingleHHVStdSPEMaker(FlatFieldSingleHHVSPEMaker): + +class SPEnominalStdalgorithm(SPEnominalalgorithm): """class to perform fit of the SPE signal with n and pp fixed""" - __parameters_file = "parameters_signalStd.yaml" - _reduced_name = "FlatFieldSingleStdSPE" + parameters_file = Unicode( + "parameters_SPEnominalStd.yaml", + read_only=True, + help="The name of the SPE fit parameters file", + ).tag(config=True) - def __init__(self, charge: np.ndarray, counts: np.ndarray, *args, **kwargs) -> None: + def __init__( + self, + pixels_id: np.ndarray, + charge: np.ndarray, + counts: np.ndarray, + config=None, + parent=None, + **kwargs, + ) -> None: """ Initializes a new instance of the FlatFieldSingleHHVStdSPEMaker class. @@ -965,161 +1007,133 @@ def __init__(self, charge: np.ndarray, counts: np.ndarray, *args, **kwargs) -> N *args: Additional positional arguments. **kwargs: Additional keyword arguments. """ - super().__init__(charge, counts, *args, **kwargs) + super().__init__( + pixels_id=pixels_id, + charge=charge, + counts=counts, + config=config, + parent=parent, + **kwargs, + ) self.__fix_parameters() def __fix_parameters(self) -> None: """ - Fixes the values of the n and pp parameters by setting their frozen attribute - to True. + Fixes the values of the n and pp parameters by setting their frozen attribute to True. """ - log.info("updating parameters by fixing pp and n") + self.log.info("updating parameters by fixing pp and n") pp = self._parameters["pp"] pp.frozen = True n = self._parameters["n"] n.frozen = True -class FlatFieldSingleNominalSPEMaker(FlatFieldSingleHHVSPEMaker): - """ - A class to perform a fit of the single photoelectron (SPE) signal at nominal - voltage using fitted data obtained from a 1400V run. - Inherits from FlatFieldSingleHHVSPEMaker. - Fixes the parameters n, pp, and res. - Optionally fixes the luminosity parameter. - - Args: - charge (np.ndarray): The charge values. - counts (np.ndarray): The counts values. - nectarGainSPEresult (str): The path to the fitted data obtained from a 1400V - run. - same_luminosity (bool, optional): Whether to fix the luminosity parameter. - Defaults to False. - *args: Variable length argument list. - **kwargs: Arbitrary keyword arguments. - - Attributes: - __parameters_file (str): The path to the parameters file for the fit at - nominal voltage. - _reduced_name (str): The name of the reduced data for the fit at nominal - voltage. - __same_luminosity (bool): Whether the luminosity parameter should be fixed. - __nectarGainSPEresult (QTable): The fitted data obtained from a 1400V run, - filtered for valid pixels. - - Example Usage: - # Create an instance of FlatFieldSingleNominalSPEMaker - maker = FlatFieldSingleNominalSPEMaker(charge, counts, - nectarGainSPEresult='fit_result.txt', same_luminosity=True) - - # Perform the fit on the specified pixels and return the fit results - results = maker.make(pixels_id=[1, 2, 3]) - - # Plot the fit for the specified pixels - maker.display(pixels_id=[1, 2, 3]) - """ - - __parameters_file = "parameters_signal_fromHHVFit.yaml" - _reduced_name = "FlatFieldSingleNominalSPE" +class SPEHHVStdalgorithm(SPEnominalStdalgorithm): + parameters_file = Unicode( + "parameters_SPEHHVStd.yaml", + read_only=True, + help="The name of the SPE fit parameters file", + ).tag(config=True) + tol = Float( + 1e40, + help="The tolerance used for minuit", + read_only=True, + ).tag(config=True) + + +class SPECombinedalgorithm(SPEnominalalgorithm): + parameters_file = Unicode( + "parameters_SPECombined_fromHHVFit.yaml", + read_only=True, + help="The name of the SPE fit parameters file", + ).tag(config=True) + + tol = Float( + 1e5, + help="The tolerance used for minuit", + read_only=True, + ).tag(config=True) + + SPE_result = Path( + help="the path of the SPE result container computed with very high voltage data", + ).tag(config=True) + + same_luminosity = Bool( + help="if the luminosity is the same between high voltage and low voltage runs", + default_value=False, + ).tag(config=True) def __init__( self, + pixels_id: np.ndarray, charge: np.ndarray, counts: np.ndarray, - nectarGainSPEresult: str, - same_luminosity: bool = False, - *args, + config=None, + parent=None, **kwargs, ) -> None: """ - Initializes an instance of FlatFieldSingleNominalSPEMaker. + Initializes a new instance of the FlatFieldSingleHHVStdSPEMaker class. Args: - charge (np.ndarray): The charge values. - counts (np.ndarray): The counts values. - nectarGainSPEresult (str): The path to the fitted data obtained from a - 1400V run. - same_luminosity (bool, optional): Whether to fix the luminosity - parameter. Defaults to False. - *args: Variable length argument list. - **kwargs: Arbitrary keyword arguments. - """ - super().__init__(charge, counts, *args, **kwargs) - self.__fix_parameters(same_luminosity) - self.__same_luminosity = same_luminosity - self.__nectarGainSPEresult = self._read_SPEresult(nectarGainSPEresult) - if len(self.__nectarGainSPEresult) == 0: - log.warning( - "The intersection between pixels id from the data and those valid " - "from the SPE fit result is empty" - ) - - @property - def nectarGainSPEresult(self): - """ - QTable: The fitted data obtained from a 1400V run, filtered for valid pixels. - """ - return copy.deepcopy(self.__nectarGainSPEresult) - - @property - def same_luminosity(self): - """ - bool: Whether the luminosity parameter should be fixed. - """ - return copy.deepcopy(self.__same_luminosity) - - def _read_SPEresult(self, nectarGainSPEresult: str): + charge (np.ndarray): The charge data. + counts (np.ndarray): The counts data. + *args: Additional positional arguments. + **kwargs: Additional keyword arguments. """ - Reads the fitted data obtained from a 1400V run and returns a filtered table - of valid pixels. - - Args: - nectarGainSPEresult (str): The path to the fitted data obtained from a - 1400V run. + super().__init__( + pixels_id=pixels_id, + charge=charge, + counts=counts, + config=config, + parent=parent, + **kwargs, + ) - Returns: - QTable: The filtered table of valid pixels. - """ - table = QTable.read(nectarGainSPEresult, format="ascii.ecsv") - table = table[table["is_valid"]] - argsort = [] - mask = [] - for _id in self._pixels_id: - if _id in table["pixels_id"]: - argsort.append(np.where(_id == table["pixels_id"])[0][0]) - mask.append(True) - else: - mask.append(False) - self._pixels_id = self._pixels_id[np.array(mask)] - return table[np.array(argsort)] + self.__fix_parameters() + self._nectarGainSPEresult = SPEfitContainer.from_hdf5(self.SPE_result) + if ( + len( + pixels_id[ + np.in1d( + pixels_id, + self._nectarGainSPEresult.pixels_id[ + self._nectarGainSPEresult.is_valid + ], + ) + ] + ) + == 0 + ): + self.log.warning( + "The intersection between pixels id from the data and those valid from the SPE fit result is empty" + ) - def __fix_parameters(self, same_luminosity: bool) -> None: + def __fix_parameters(self) -> None: """ Fixes the parameters n, pp, res, and possibly luminosity. Args: same_luminosity (bool): Whether to fix the luminosity parameter. """ - log.info("updating parameters by fixing pp, n and res") + self.log.info("updating parameters by fixing pp, n and res") pp = self._parameters["pp"] pp.frozen = True n = self._parameters["n"] n.frozen = True resolution = self._parameters["resolution"] resolution.frozen = True - if same_luminosity: - log.info("fixing luminosity") + if self.same_luminosity: + self.log.info("fixing luminosity") luminosity = self._parameters["luminosity"] luminosity.frozen = True def _make_fit_array_from_parameters(self, pixels_id=None, **kwargs): """ - Generates the fit array from the fixed parameters and the fitted data - obtained from a 1400V run. + Generates the fit array from the fixed parameters and the fitted data obtained from a 1400V run. Args: - pixels_id (array-like, optional): The pixels to generate the fit array - for. Defaults to None. + pixels_id (array-like, optional): The pixels to generate the fit array for. Defaults to None. **kwargs: Arbitrary keyword arguments. Returns: @@ -1127,7 +1141,7 @@ def _make_fit_array_from_parameters(self, pixels_id=None, **kwargs): """ return super()._make_fit_array_from_parameters( pixels_id=pixels_id, - nectarGainSPEresult=self.__nectarGainSPEresult, + nectarGainSPEresult=self._nectarGainSPEresult, **kwargs, ) @@ -1141,8 +1155,7 @@ def _update_parameters( **kwargs, ): """ - Updates the parameters with the fixed values from the fitted data obtained - from a 1400V run. + Updates the parameters with the fixed values from the fitted data obtained from a 1400V run. Args: parameters (Parameters): The parameters to update. @@ -1155,18 +1168,41 @@ def _update_parameters( Returns: dict: The updated parameters. """ - param = super()._update_parameters(parameters, charge, counts, **kwargs) + param = super(__class__, __class__)._update_parameters( + parameters, charge, counts, **kwargs + ) luminosity = param["luminosity"] resolution = param["resolution"] pp = param["pp"] n = param["n"] - index = np.where(pixel_id == nectarGainSPEresult["pixels_id"])[0][0] + index = np.where(pixel_id == nectarGainSPEresult.pixels_id)[0][0] - resolution.value = nectarGainSPEresult[index]["resolution"].value - pp.value = nectarGainSPEresult[index]["pp"].value - n.value = nectarGainSPEresult[index]["n"].value + resolution.value = nectarGainSPEresult.resolution[index] + pp.value = nectarGainSPEresult.pp[index] + n.value = nectarGainSPEresult.n[index]["n"] if luminosity.frozen: - luminosity.value = nectarGainSPEresult[index]["luminosity"].value + luminosity.value = nectarGainSPEresult.luminosity[index].value return param + + def run( + self, + pixels_id: np.ndarray = None, + display: bool = True, + **kwargs, + ) -> np.ndarray: + if pixels_id is None: + pixels_id = self._nectarGainSPEresult.pixels_id[ + self._nectarGainSPEresult.is_valid + ] + else: + pixels_id = np.asarray(pixels_id)[ + np.in1d( + pixels_id, + self._nectarGainSPEresult.pixels_id[ + self._nectarGainSPEresult.is_valid + ], + ) + ] + return super().run(pixels_id=pixels_id, display=display, **kwargs) diff --git a/src/nectarchain/makers/component/waveforms_component.py b/src/nectarchain/makers/component/waveformsComponent.py similarity index 89% rename from src/nectarchain/makers/component/waveforms_component.py rename to src/nectarchain/makers/component/waveformsComponent.py index f6bf50b2..c8dd4a5f 100644 --- a/src/nectarchain/makers/component/waveforms_component.py +++ b/src/nectarchain/makers/component/waveformsComponent.py @@ -1,5 +1,11 @@ -import copy import logging + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + + +import copy from argparse import ArgumentError import numpy as np @@ -12,14 +18,13 @@ from ...data.container import WaveformsContainer, WaveformsContainers from .core import ArrayDataComponent -logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") -log = logging.getLogger(__name__) -log.handlers = logging.getLogger("__main__").handlers - __all__ = ["WaveformsComponent"] class WaveformsComponent(ArrayDataComponent): + SubComponents = copy.deepcopy(ArrayDataComponent.SubComponents) + SubComponents.read_only = True + def __init__(self, subarray, config=None, parent=None, *args, **kwargs): super().__init__( subarray=subarray, config=config, parent=parent, *args, **kwargs @@ -32,9 +37,9 @@ def __init__(self, subarray, config=None, parent=None, *args, **kwargs): @staticmethod def create_from_events_list( events_list: list, - run_number: int, - npixels: int, - nsamples: int, + run_number: np.uint16, + npixels: np.uint16, + nsamples: np.uint8, subarray: SubarrayDescription, pixels_id: int, tel_id: int = None, @@ -42,8 +47,7 @@ def create_from_events_list( """Create a container for the extracted waveforms from a list of events. Args: - events_list (list[NectarCAMDataContainer]): A list of events to extract - waveforms from. + events_list (list[NectarCAMDataContainer]): A list of events to extract waveforms from. run_number (int): The ID of the run to be loaded. npixels (int): The number of pixels in the waveforms. nsamples (int): The number of samples in the waveforms. @@ -51,8 +55,7 @@ def create_from_events_list( pixels_id (int): The ID of the pixels to extract waveforms from. Returns: - WaveformsContainer: A container object that contains the extracted - waveforms and other relevant information. + WaveformsContainer: A container object that contains the extracted waveforms and other relevant information. """ if tel_id is None: tel_id = __class__.TEL_ID.default_value @@ -127,8 +130,7 @@ def __call__(self, event: NectarCAMDataContainer, *args, **kwargs): """Process an event and extract waveforms. Args: - event (NectarCAMDataContainer): The event from which to process and extract - waveforms. + event (NectarCAMDataContainer): The event to process and extract waveforms from. trigger (EventType): The type of trigger for the event. """ @@ -156,18 +158,21 @@ def finish(self, *args, **kwargs): trigger_type (EventType): The selected trigger types. Returns: - list[WaveformsContainer]: A list of output containers for the selected - trigger types. + list[WaveformsContainer]: A list of output containers for the selected trigger types. """ output = WaveformsContainers() for i, trigger in enumerate(self.trigger_list): waveformsContainer = WaveformsContainer( - run_number=self._run_number, - npixels=self._npixels, - nsamples=self._nsamples, + run_number=WaveformsContainer.fields["run_number"].type( + self._run_number + ), + npixels=WaveformsContainer.fields["npixels"].type(self._npixels), + nsamples=WaveformsContainer.fields["nsamples"].type(self._nsamples), # subarray=self.subarray, camera=self.CAMERA_NAME, - pixels_id=self._pixels_id, + pixels_id=WaveformsContainer.fields["pixels_id"].dtype.type( + self._pixels_id + ), nevents=self.nevents(trigger), wfs_hg=self.wfs_hg(trigger), wfs_lg=self.wfs_lg(trigger), @@ -190,8 +195,7 @@ def sort(waveformsContainer: WaveformsContainer, method: str = "event_id"): """Sort the waveformsContainer based on a specified method. Args: - waveformsContainer (WaveformsContainer): The waveformsContainer to be - sorted. + waveformsContainer (WaveformsContainer): The waveformsContainer to be sorted. method (str, optional): The sorting method. Defaults to 'event_id'. Returns: @@ -201,7 +205,6 @@ def sort(waveformsContainer: WaveformsContainer, method: str = "event_id"): run_number=waveformsContainer.run_number, npixels=waveformsContainer.npixels, nsamples=waveformsContainer.nsamples, - subarray=waveformsContainer.subarray, camera=waveformsContainer.camera, pixels_id=waveformsContainer.pixels_id, nevents=waveformsContainer.nevents, @@ -214,7 +217,6 @@ def sort(waveformsContainer: WaveformsContainer, method: str = "event_id"): in [ "run_number", "npixels", - "nsamples", "subarray", "camera", "pixels_id", @@ -233,10 +235,8 @@ def select_waveforms_hg( """Select HIGH GAIN waveforms from the container. Args: - waveformsContainer (WaveformsContainer): The container object that - contains the waveforms. - pixel_id (np.ndarray): An array of pixel IDs to select specific waveforms - from the container. + waveformsContainer (WaveformsContainer): The container object that contains the waveforms. + pixel_id (np.ndarray): An array of pixel IDs to select specific waveforms from the container. Returns: np.ndarray: An array of selected waveforms from the container. @@ -254,10 +254,8 @@ def select_waveforms_lg( """Select LOW GAIN waveforms from the container. Args: - waveformsContainer (WaveformsContainer): The container object that - contains the waveforms. - pixel_id (np.ndarray): An array of pixel IDs to select specific waveforms - from the container. + waveformsContainer (WaveformsContainer): The container object that contains the waveforms. + pixel_id (np.ndarray): An array of pixel IDs to select specific waveforms from the container. Returns: np.ndarray: An array of selected waveforms from the container. @@ -292,29 +290,27 @@ def wfs_hg(self, trigger: EventType): Returns the waveform data for the specified trigger type. Args: - trigger (EventType): The type of trigger for which the waveform data is - requested. + trigger (EventType): The type of trigger for which the waveform data is requested. Returns: An array of waveform data for the specified trigger type. """ return np.array( - self.__wfs_hg[__class__._get_name_trigger(trigger)], dtype=np.uint16 + self.__wfs_hg[__class__._get_name_trigger(trigger)], + dtype=WaveformsContainer.fields["wfs_hg"].dtype, ) def wfs_lg(self, trigger: EventType): """ - Returns the waveform data for the specified trigger type in the low gain - channel. + Returns the waveform data for the specified trigger type in the low gain channel. Args: - trigger (EventType): The type of trigger for which the waveform data is - requested. + trigger (EventType): The type of trigger for which the waveform data is requested. Returns: - An array of waveform data for the specified trigger type in the low gain - channel. + An array of waveform data for the specified trigger type in the low gain channel. """ return np.array( - self.__wfs_lg[__class__._get_name_trigger(trigger)], dtype=np.uint16 + self.__wfs_lg[__class__._get_name_trigger(trigger)], + dtype=WaveformsContainer.fields["wfs_lg"].dtype, ) diff --git a/src/nectarchain/makers/core.py b/src/nectarchain/makers/core.py index 7ef1cfc4..dec9a9f3 100644 --- a/src/nectarchain/makers/core.py +++ b/src/nectarchain/makers/core.py @@ -1,9 +1,17 @@ -import copy import logging + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + +import copy +import os import pathlib import numpy as np +from ctapipe.containers import Container from ctapipe.core import Component, Tool +from ctapipe.core.container import FieldValidationError from ctapipe.core.traits import ( Bool, ComponentNameList, @@ -14,32 +22,25 @@ ) from ctapipe.io import HDF5TableWriter from ctapipe.io.datawriter import DATA_MODEL_VERSION -from ctapipe_io_nectarcam import NectarCAMEventSource +from tables.exceptions import HDF5ExtError from tqdm.auto import tqdm from ..data import DataManagement +from ..data.container import LightNectarCAMEventSource from ..data.container.core import NectarCAMContainer, TriggerMapContainer -from .component.core import ( - NectarCAMComponent, - get_configurable_traits, - get_valid_component, -) - -logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") -log = logging.getLogger(__name__) -log.handlers = logging.getLogger("__main__").handlers +from ..utils import ComponentUtils +from .component import * +from .component import NectarCAMComponent, get_valid_component __all__ = ["EventsLoopNectarCAMCalibrationTool"] -"""The code snippet is a part of a class hierarchy for data processing. -It includes the `BaseMaker` abstract class, the `EventsLoopMaker` and -`ArrayDataMaker` subclasses. +"""The code snippet is a part of a class hierarchy for data processing. +It includes the `BaseMaker` abstract class, the `EventsLoopMaker` and `ArrayDataMaker` subclasses. These classes are used to perform computations on data from a specific run.""" class BaseNectarCAMCalibrationTool(Tool): - """Mother class for all the makers, the role of makers is to do computation on - the data.""" + """Mother class for all the makers, the role of makers is to do computation on the data.""" name = "BaseNectarCAMCalibration" @@ -50,29 +51,26 @@ class BaseNectarCAMCalibrationTool(Tool): @staticmethod def load_run( run_number: int, max_events: int = None, run_file: str = None - ) -> NectarCAMEventSource: - """Static method to load from $NECTARCAMDATA directory data - for specified run with max_events. + ) -> LightNectarCAMEventSource: + """Static method to load from $NECTARCAMDATA directory data for specified run with max_events Args:self.__run_number = run_number run_number (int): run_id - maxevents (int, optional): max of events to be loaded. Defaults to -1, - to load everything. + maxevents (int, optional): max of events to be loaded. Defaults to -1, to load everythings. run_file (optional) : if provided, will load this run file Returns: - List[ctapipe_io_nectarcam.NectarCAMEventSource]: List of EventSource for - each run files + List[ctapipe_io_nectarcam.LightNectarCAMEventSource]: List of EventSource for each run files """ # Load the data from the run file. if run_file is None: generic_filename, _ = DataManagement.findrun(run_number) log.info(f"{str(generic_filename)} will be loaded") - eventsource = NectarCAMEventSource( + eventsource = LightNectarCAMEventSource( input_url=generic_filename, max_events=max_events ) else: log.info(f"{run_file} will be loaded") - eventsource = NectarCAMEventSource( + eventsource = LightNectarCAMEventSource( input_url=run_file, max_events=max_events ) return eventsource @@ -84,8 +82,7 @@ class EventsLoopNectarCAMCalibrationTool(BaseNectarCAMCalibrationTool): Args: run_number (int): The ID of the run to be processed. - max_events (int, optional): The maximum number of events to be loaded. - Defaults to None. + max_events (int, optional): The maximum number of events to be loaded. Defaults to None. run_file (optional): The specific run file to be loaded. Example Usage: @@ -107,16 +104,17 @@ class EventsLoopNectarCAMCalibrationTool(BaseNectarCAMCalibrationTool): ("r", "run-number"): "EventsLoopNectarCAMCalibrationTool.run_number", ("m", "max-events"): "EventsLoopNectarCAMCalibrationTool.max_events", ("o", "output"): "EventsLoopNectarCAMCalibrationTool.output_path", + "events-per-slice": "EventsLoopNectarCAMCalibrationTool.events_per_slice", } flags = { "overwrite": ( - {"HDF5TableWriter": {"overwrite": True}}, + {"HDF5TableWriter": {"overwrite": False}}, "Overwrite output file if it exists", ), **flag( "progress", - "ProcessorTool.progress_bar", + "EventsLoopNectarCAMCalibrationTool.progress_bar", "show a progress bar during event processing", "don't show a progress bar during event processing", ), @@ -126,19 +124,21 @@ class EventsLoopNectarCAMCalibrationTool(BaseNectarCAMCalibrationTool): [ HDF5TableWriter, ] - + classes_with_traits(NectarCAMEventSource) + + classes_with_traits(LightNectarCAMEventSource) + classes_with_traits(NectarCAMComponent) ) - output_path = Path( - help="output filename", - default_value=pathlib.Path("/tmp/EventsLoopNectarCAMCalibrationTool.h5"), - ).tag(config=True) - run_number = Integer(help="run number to be treated", default_value=-1).tag( config=True ) + output_path = Path( + help="output filename", + default_value=pathlib.Path( + f"{os.environ.get('NECTARCAMDATA','/tmp')}/runs/{name}_run{run_number.default_value}.h5" + ), + ).tag(config=True) + max_events = Integer( help="maximum number of events to be loaded", default_value=None, @@ -156,48 +156,153 @@ class EventsLoopNectarCAMCalibrationTool(BaseNectarCAMCalibrationTool): help="List of Component names to be apply, the order will be respected", ).tag(config=True) + events_per_slice = Integer( + help="number of events that will be treat before to pull the buffer and write to disk, if None, all the events will be loaded", + default_value=None, + allow_none=True, + ).tag(config=True) + def __new__(cls, *args, **kwargs): - """This method is used to pass to the current instance of Tool the traits - defined in the components provided in the componentsList trait. - WARNING : This method is maybe not the best way to do it, need to discuss - with ctapipe developers. + """This method is used to pass to the current instance of Tool the traits defined + in the components provided in the componentsList trait. + WARNING : This method is maybe not the best way to do it, need to discuss with ctapipe developpers. """ _cls = super(EventsLoopNectarCAMCalibrationTool, cls).__new__( cls, *args, **kwargs ) log.warning( - "the componentName in componentsList must be defined in the " - "nectarchain.makers.component module, otherwise the import of the " - "componentName will raise an error" + "the componentName in componentsList must be defined in the nectarchain.makers.component module, otherwise the import of the componentName will raise an error" ) for componentName in _cls.componentsList: - configurable_traits = get_configurable_traits(eval(componentName)) + class_name = ComponentUtils.get_class_name_from_ComponentName(componentName) + configurable_traits = ComponentUtils.get_configurable_traits(class_name) _cls.add_traits(**configurable_traits) _cls.aliases.update( {key: f"{componentName}.{key}" for key in configurable_traits.keys()} ) return _cls - def _load_eventsource(self): + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + if not ("output_path" in kwargs.keys()): + self._init_output_path() + + def _init_output_path(self): + self.output_path = pathlib.Path( + f"{os.environ.get('NECTARCAMDATA','/tmp')}/runs/{self.name}_run{self.run_number}.h5" + ) + + def _load_eventsource(self, *args, **kwargs): + self.log.debug("loading event source") self.event_source = self.enter_context( self.load_run(self.run_number, self.max_events, run_file=self.run_file) ) def _get_provided_component_kwargs(self, componentName: str): - component_kwargs = get_configurable_traits(eval(componentName)) + class_name = ComponentUtils.get_class_name_from_ComponentName(componentName) + component_kwargs = ComponentUtils.get_configurable_traits(class_name) output_component_kwargs = {} for key in component_kwargs.keys(): if hasattr(self, key): output_component_kwargs[key] = getattr(self, key) return output_component_kwargs + def _init_writer(self, sliced: bool = False, slice_index: int = 0): + if hasattr(self, "writer"): + self.writer.close() + + if sliced: + if slice_index == 0: + if self.overwrite: + try: + log.info( + f"overwrite set to true, trying to remove file {self.output_path}" + ) + os.remove(self.output_path) + log.info(f"{self.output_path} removed") + except OSError: + pass + else: + if os.path.exists(self.output_path): + raise Exception( + f"file {self.output_path} does exist,\n set overwrite to True if you want to overwrite" + ) + self.log.info( + f"initilization of writter in sliced mode (slice index = {slice_index})" + ) + # self.output_path = pathlib.Path(f"{self.output_path.parent)}/{self.output_path.stem}_1{self.output_path.suffix}" + group_name = f"data_{slice_index}" + mode = "a" + else: + self.log.info("initilization of writter in full mode") + if self.overwrite: + try: + log.info( + f"overwrite set to true, trying to remove file {self.output_path}" + ) + os.remove(self.output_path) + log.info(f"{self.output_path} removed") + except OSError: + pass + else: + if os.path.exists(self.output_path): + raise Exception( + f"file {self.output_path} does exist,\n set overwrite to True if you want to overwrite" + ) + group_name = "data" + mode = "w" + try: + os.makedirs(self.output_path.parent, exist_ok=True) + self.writer = self.enter_context( + HDF5TableWriter( + filename=self.output_path, # pathlib.Path(f"{os.environ.get('NECTARCAMDATA',self.output_path.parent)}/{self.output_path.stem}_run{self.run_number}{self.output_path.suffix}"), + parent=self, + mode=mode, + group_name=group_name, + ) + ) + except HDF5ExtError as err: + self.log.warning(err.args[0], exc_info=True) + self.log.warning("retry with w mode instead of a") + self.writer = self.enter_context( + HDF5TableWriter( + filename=self.output_path, # pathlib.Path(f"{os.environ.get('NECTARCAMDATA',self.output_path.parent)}/{self.output_path.stem}_run{self.run_number}{self.output_path.suffix}"), + parent=self, + mode="w", + group_name=group_name, + ) + ) + except Exception as err: + self.log.error(err, exc_info=True) + raise err + def setup(self, *args, **kwargs): + self.log.info("setup of the Tool") if self.run_number == -1: raise Exception("run_number need to be set up") - self._load_eventsource() + self._setup_eventsource(*args, **kwargs) + + self._setup_components(*args, **kwargs) + + if self.output_path.exists() and self.overwrite: + os.remove(self.output_path) + + self._init_writer(sliced=not (self.events_per_slice is None), slice_index=1) + + self._n_traited_events = 0 + + # self.comp = MyComponent(parent=self) + # self.comp2 = SecondaryMyComponent(parent=self) + # self.comp3 = TelescopeWiseComponent(parent=self, subarray=subarray) + # self.advanced = AdvancedComponent(parent=self) + + def _setup_eventsource(self, *args, **kwargs): + self._load_eventsource(*args, **kwargs) self.__npixels = self._event_source.camera_config.num_pixels self.__pixels_id = self._event_source.camera_config.expected_pixels_id + def _setup_components(self, *args, **kwargs): + self.log.info("setup of components") self.components = [] for componentName in self.componentsList: if componentName in get_valid_component(): @@ -213,22 +318,6 @@ def setup(self, *args, **kwargs): # ) ) - self.writer = self.enter_context( - HDF5TableWriter( - filename=pathlib.Path( - f"{self.output_path.parent}/{self.output_path.stem}" - f"_{self.run_number}{self.output_path.suffix}" - ), - parent=self, - group_name="data", - ) - ) - - # self.comp = MyComponent(parent=self) - # self.comp2 = SecondaryMyComponent(parent=self) - # self.comp3 = TelescopeWiseComponent(parent=self, subarray=subarray) - # self.advanced = AdvancedComponent(parent=self) - def start( self, n_events=np.inf, @@ -241,19 +330,17 @@ def start( Method to extract data from the EventSource. Args: - n_events (int, optional): The maximum number of events to process. - Default is np.inf. - restart_from_begining (bool, optional): Whether to restart the event - source reader. Default is False. + n_events (int, optional): The maximum number of events to process. Default is np.inf. + restart_from_begining (bool, optional): Whether to restart the event source reader. Default is False. *args: Additional arguments that can be passed to the method. **kwargs: Additional keyword arguments that can be passed to the method. Returns: The output container created by the _make_output_container method. """ - if ~np.isfinite(n_events): + if ~np.isfinite(n_events) and (self.events_per_slice is None): self.log.warning( - "no needed events number specified, it may cause a memory error" + "neither needed events number specified or events per slice, it may cause a memory error" ) # if isinstance(trigger_type, EventType) or trigger_type is None: # trigger_type = [trigger_type] @@ -264,54 +351,85 @@ def start( self.log.debug("restart from begining : creation of the EventSource reader") self._load_eventsource() - n_traited_events = 0 + n_events_in_slice = 0 + slice_index = 1 for i, event in enumerate( tqdm( self._event_source, desc=self._event_source.__class__.__name__, - total=min(self._event_source.max_events, n_events), + total=len(self._event_source) + if self._event_source.max_events is None + else int(np.min((self._event_source.max_events, n_events))), unit="ev", disable=not self.progress_bar, ) ): - if i % 100 == 0: - self.log.info(f"reading event number {i}") + # if i % 100 == 0: + # self.log.info(f"reading event number {i}") for component in self.components: component(event, *args, **kwargs) - n_traited_events += 1 - if n_traited_events >= n_events: + self._n_traited_events += 1 + n_events_in_slice += 1 + if self._n_traited_events >= n_events: break + if ( + not (self.events_per_slice is None) + and n_events_in_slice >= self.events_per_slice + ): + self.log.info(f"slice number {slice_index} is full, pulling buffer") + self._finish_components(*args, **kwargs) + self.writer.close() + slice_index += 1 + self._init_writer(sliced=True, slice_index=slice_index) + self._setup_components() + n_events_in_slice = 0 + + def finish(self, return_output_component=False, *args, **kwargs): + self.log.info("finishing Tool") + + output = self._finish_components(*args, **kwargs) - def finish(self, *args, **kwargs): - # self.write = self.enter_context( - # HDF5TableWriter(filename=filename, parent=self) - # ) + self.writer.close() + super().finish() + self.log.warning("Shutting down.") + if return_output_component: + return output + + def _finish_components(self, *args, **kwargs): + self.log.info("finishing components and writting to output file") output = [] for component in self.components: output.append(component.finish(*args, **kwargs)) log.info(output) - for _output in output: - if isinstance(_output, NectarCAMContainer): + for i, _output in enumerate(output): + if not (_output is None): + self._write_container(_output, i) + return output + + def _write_container(self, container: Container, index_component: int = 0) -> None: + try: + container.validate() + if isinstance(container, NectarCAMContainer): self.writer.write( - table_name=str(_output.__class__.__name__), - containers=_output, + table_name=str(container.__class__.__name__), + containers=container, ) - elif isinstance(_output, TriggerMapContainer): - for i, key in enumerate(_output.containers.keys()): + elif isinstance(container, TriggerMapContainer): + for key in container.containers.keys(): self.writer.write( - table_name=f"{_output.containers[key].__class__.__name__}_" - f"{i}/ {key.name}", - containers=_output.containers[key], + table_name=f"{container.containers[key].__class__.__name__}_{index_component}/{key.name}", + containers=container.containers[key], ) else: raise TypeError( - "component output must be an instance of TriggerMapContainer or " - "NectarCAMContainer" + "component output must be an instance of TriggerMapContainer or NectarCAMContainer" ) - - self.writer.close() - super().finish() - self.log.warning("Shutting down.") + except FieldValidationError as e: + log.warning(e, exc_info=True) + self.log.warning("the container has not been written") + except Exception as e: + self.log.error(e, exc_info=True) + raise e @property def event_source(self): @@ -323,15 +441,15 @@ def event_source(self): @event_source.setter def event_source(self, value): """ - Setter method to set a new NectarCAMEventSource to the _reader attribute. + Setter method to set a new LightNectarCAMEventSource to the _reader attribute. Args: - value: a NectarCAMEventSource instance. + value: a LightNectarCAMEventSource instance. """ - if isinstance(value, NectarCAMEventSource): + if isinstance(value, LightNectarCAMEventSource): self._event_source = value else: - raise TypeError("The reader must be a NectarCAMEventSource") + raise TypeError("The reader must be a LightNectarCAMEventSource") @property def _npixels(self): diff --git a/src/nectarchain/makers/extractor/utils.py b/src/nectarchain/makers/extractor/utils.py index 0fab1212..0fc1804a 100644 --- a/src/nectarchain/makers/extractor/utils.py +++ b/src/nectarchain/makers/extractor/utils.py @@ -12,7 +12,7 @@ class CtapipeExtractor: A class to extract the image and peak time from a DL1CameraContainer object. """ - def get_image_peak_time(cameraContainer): + def get_image_peak_time(cameraContainer: DL1CameraContainer): """ Extracts the image and peak time from a DL1CameraContainer object. @@ -23,3 +23,15 @@ def get_image_peak_time(cameraContainer): tuple: A tuple containing the image and peak time values from the container. """ return cameraContainer.image, cameraContainer.peak_time + + def get_extractor_kwargs_str(extractor_kwargs): + if len(extractor_kwargs) == 0: + str_extractor_kwargs = "" + else: + extractor_kwargs_list = [ + f"{key}_{value}" for key, value in extractor_kwargs.items() + ] + str_extractor_kwargs = extractor_kwargs_list[0] + for item in extractor_kwargs_list[1:]: + str_extractor_kwargs += f"_{item}" + return str_extractor_kwargs diff --git a/src/nectarchain/makers/waveformsMakers.py b/src/nectarchain/makers/waveformsMakers.py new file mode 100644 index 00000000..fb0f94cb --- /dev/null +++ b/src/nectarchain/makers/waveformsMakers.py @@ -0,0 +1,39 @@ +import logging + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + +import os +import pathlib + +from ctapipe.core.traits import ComponentNameList + +from .component import NectarCAMComponent +from .core import EventsLoopNectarCAMCalibrationTool + +__all__ = ["WaveformsNectarCAMCalibrationTool"] + + +class WaveformsNectarCAMCalibrationTool(EventsLoopNectarCAMCalibrationTool): + """class use to make the waveform extraction from event read from r0 data""" + + name = "WaveformsNectarCAMCalibration" + + componentsList = ComponentNameList( + NectarCAMComponent, + default_value=["WaveformsComponent"], + help="List of Component names to be apply, the order will be respected", + ).tag(config=True) + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + def _init_output_path(self): + if self.max_events is None: + filename = f"{self.name}_run{self.run_number}.h5" + else: + filename = f"{self.name}_run{self.run_number}_maxevents{self.max_events}.h5" + self.output_path = pathlib.Path( + f"{os.environ.get('NECTARCAMDATA','/tmp')}/runs/waveforms/{filename}" + ) diff --git a/src/nectarchain/makers/waveforms_makers.py b/src/nectarchain/makers/waveforms_makers.py deleted file mode 100644 index d82c87c9..00000000 --- a/src/nectarchain/makers/waveforms_makers.py +++ /dev/null @@ -1,22 +0,0 @@ -import logging - -from ctapipe.core.traits import ComponentNameList - -from .component import NectarCAMComponent -from .core import EventsLoopNectarCAMCalibrationTool - -logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") -log = logging.getLogger(__name__) -log.handlers = logging.getLogger("__main__").handlers - -__all__ = ["WaveformsNectarCAMCalibrationTool"] - - -class WaveformsNectarCAMCalibrationTool(EventsLoopNectarCAMCalibrationTool): - """class used to make the waveform extraction from event read from r0 data""" - - componentsList = ComponentNameList( - NectarCAMComponent, - default_value=["WaveformsComponent"], - help="List of Component names to be applied, the order will be respected", - ).tag(config=True) diff --git a/src/nectarchain/user_scripts/ggrolleron/gain_PhotoStat_computation.py b/src/nectarchain/user_scripts/ggrolleron/gain_PhotoStat_computation.py index f5a9bc1d..582fb06f 100644 --- a/src/nectarchain/user_scripts/ggrolleron/gain_PhotoStat_computation.py +++ b/src/nectarchain/user_scripts/ggrolleron/gain_PhotoStat_computation.py @@ -1,73 +1,85 @@ +import json import logging import os import sys from pathlib import Path -import numpy as np - -os.makedirs(os.environ.get("NECTARCHAIN_LOG"), exist_ok=True) - # to quiet numba -logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +os.makedirs(os.environ.get("NECTARCHAIN_LOG"), exist_ok=True) logging.getLogger("numba").setLevel(logging.WARNING) +logging.basicConfig( + format="%(asctime)s %(name)s %(levelname)s %(message)s", + level=logging.DEBUG, + filename=f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{Path(__file__).stem}_{os.getpid()}.log", +) +log = logging.getLogger(__name__) import argparse +import copy +import glob -from astropy.table import QTable - -from nectarchain.makers.calibration.gain.photostatistic_makers import ( - PhotoStatisticMaker, -) +from nectarchain.data.management import DataManagement +from nectarchain.makers.calibration import PhotoStatisticNectarCAMCalibrationTool +from nectarchain.makers.extractor.utils import CtapipeExtractor parser = argparse.ArgumentParser( - prog="gain_PhotoStat_computation.py", - description="compute high gain and low gain with Photo-statistic method, it need a pedestal run and a FF run with a SPE fit results (for resolution value needed in this method). Output data will be saved in $NECTARCAMDATA/../PhotoStat/data/PhotoStat-FF{FF_run_number}-ped{ped_run_number}-SPEres{SPE_fit_results_tag}-{chargeExtractorPath}/", + prog="gain_SPEfit_computation.py", + description=f"compute high and low gain with the Photo-statistic method, output data are saved in $NECTARCAMDATA/../PhotoStat/", ) - # run numbers -parser.add_argument("-p", "--ped_run_number", help="ped run", required=True, type=int) -parser.add_argument("-f", "--FF_run_number", help="FF run", required=True, type=int) parser.add_argument( - "--SPE_fit_results", - help="SPE fit results path for accessing SPE resolution", - type=str, - required=True, + "--FF_run_number", nargs="+", default=[], help="run(s) list", type=int +) +parser.add_argument( + "--Ped_run_number", nargs="+", default=[], help="run(s) list", type=int ) -# tag for SPE fit results propagation +# max events to be loaded parser.add_argument( - "--SPE_fit_results_tag", - help="SPE fit results tag for propagate the SPE result to output, this tag will be used to setup the path where output data will be saved, see help for description", - type=str, - default="", + "-m", + "--max_events", + nargs="+", + # default=[], + help="max events to be load", + type=int, ) +# boolean arguments parser.add_argument( - "--overwrite", + "--reload_events", action="store_true", default=False, - help="to force overwrite files on disk", + help="to force re-computation of waveforms from fits.fz files", ) -# for plotting correlation parser.add_argument( - "--correlation", + "--overwrite", action="store_true", - default=True, - help="to plot correlation between SPE gain computation and Photo-statistic gain resluts", + default=False, + help="to force overwrite files on disk", ) + # extractor arguments parser.add_argument( - "--chargeExtractorPath", - help="charge extractor path where charges are saved", + "--method", + choices=[ + "FullWaveformSum", + "FixedWindowSum", + "GlobalPeakWindowSum", + "LocalPeakWindowSum", + "SlidingWindowMaxSum", + "TwoPassWindowSum", + ], + default="LocalPeakWindowSum", + help="charge extractor method", type=str, ) - parser.add_argument( - "--FFchargeExtractorWindowLength", - help="charge extractor window length in ns", - type=int, + "--extractor_kwargs", + default={"window_width": 16, "window_shift": 4}, + help="charge extractor kwargs", + type=json.loads, ) # verbosity argument @@ -75,76 +87,118 @@ "-v", "--verbosity", help="set the verbosity level of logger", - default="info", - choices=["fatal", "debug", "info", "warning"], + default="INFO", + choices=["DEBUG", "INFO", "WARN", "ERROR", "CRITICAL"], type=str, ) +# output figures path +parser.add_argument( + "--figpath", + type=str, + default=f"{os.environ.get('NECTARCHAIN_FIGURES','/tmp')}/", + help="output figure path", +) + +# pixels selected +parser.add_argument( + "-p", "--asked_pixels_id", nargs="+", default=None, help="pixels selected", type=int +) + +parser.add_argument( + "--HHV_run_number", + help="HHV run number of which the SPE fit has ever been performed", + type=int, +) + +args = parser.parse_args() + -def main(args): - figpath = os.environ.get("NECTARCHAIN_FIGURES") +def main( + log, + **kwargs, +): + FF_run_number = kwargs.pop("FF_run_number") + Ped_run_number = kwargs.pop("Ped_run_number") - photoStat_FFandPed = PhotoStatisticMaker.create_from_run_numbers( - FFrun=args.FF_run_number, - Pedrun=args.ped_run_number, - SPE_resolution=args.SPE_fit_results, - method=args.chargeExtractorPath, - FFchargeExtractorWindowLength=args.FFchargeExtractorWindowLength, + if len(FF_run_number) != len(Ped_run_number): + raise Exception("The number of FF and Ped runs must be the same") + + max_events = kwargs.pop("max_events", [None for i in range(len(FF_run_number))]) + if max_events is None: + max_events = [None for i in range(len(FF_run_number))] + + log.info(f"max_events : {max_events}") + + figpath = args.figpath + + str_extractor_kwargs = CtapipeExtractor.get_extractor_kwargs_str( + args.extractor_kwargs ) - photoStat_FFandPed.make() - photoStat_FFandPed.save( - f"{os.environ.get('NECTARCAMDATA')}/../PhotoStat/data/PhotoStat-FF{args.FF_run_number}-ped{args.ped_run_number}-SPEres{args.SPE_fit_results_tag}-{args.chargeExtractorPath}/", - overwrite=args.overwrite, + path = DataManagement.find_SPE_HHV( + run_number=args.HHV_run_number, + method=args.method, + str_extractor_kwargs=str_extractor_kwargs, ) - log.info(f"BF^2 HG : {np.power(np.mean(photoStat_FFandPed.BHG),2)}") - log.info(f"BF^2 LG : {np.power(np.mean(photoStat_FFandPed.BLG),2)}") - - if args.correlation: - table = QTable.read(args.SPE_fit_results, format="ascii.ecsv") - table.sort("pixels_id") - mask = np.array( - [pix in photoStat_FFandPed.pixels_id for pix in table["pixels_id"].value], - dtype=bool, - ) - fig = PhotoStatisticMaker.plot_correlation( - photoStat_FFandPed.results["high_gain"], table["high_gain"][mask] - ) - os.makedirs( - f"{figpath}/PhotoStat-FF{args.FF_run_number}-ped{args.ped_run_number}-{args.chargeExtractorPath}/", - exist_ok=True, - ) - fig.savefig( - f"{figpath}/PhotoStat-FF{args.FF_run_number}-ped{args.ped_run_number}-{args.chargeExtractorPath}/correlation_PhotoStat_SPE{args.SPE_fit_results_tag}.pdf" + if len(path) == 1: + log.info( + f"{path[0]} found associated to HHV run {args.HHV_run_number}, method {args.method} and extractor kwargs {str_extractor_kwargs}" ) + else: + _text = f"no file found in $NECTARCAM_DATA/../SPEfit associated to HHV run {args.HHV_run_number}, method {args.method} and extractor kwargs {str_extractor_kwargs}" + log.error(_text) + raise FileNotFoundError(_text) + for _FF_run_number, _Ped_run_number, _max_events in zip( + FF_run_number, Ped_run_number, max_events + ): + try: + tool = PhotoStatisticNectarCAMCalibrationTool( + progress_bar=True, + run_number=_FF_run_number, + max_events=_max_events, + Ped_run_number=_Ped_run_number, + SPE_result=path[0], + **kwargs, + ) + tool.setup() + if args.reload_events and not (_max_events is None): + _figpath = f"{figpath}/{tool.name}_run{tool.run_number}_maxevents{_max_events}_{tool.method}_{str_extractor_kwargs}" + else: + _figpath = f"{figpath}/{tool.name}_run{tool.run_number}_{tool.method}_{str_extractor_kwargs}" + tool.start() + tool.finish(figpath=_figpath) + except Exception as e: + log.warning(e, exc_info=True) if __name__ == "__main__": args = parser.parse_args() - logginglevel = logging.FATAL - if args.verbosity == "warning": - logginglevel = logging.WARNING - elif args.verbosity == "info": - logginglevel = logging.INFO - elif args.verbosity == "debug": - logginglevel = logging.DEBUG - - os.makedirs(f"{os.environ.get('NECTARCHAIN_LOG')}/{os.getpid()}/figures") + kwargs = copy.deepcopy(vars(args)) + + kwargs["log_level"] = args.verbosity + + os.makedirs(f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{os.getpid()}/figures") logging.basicConfig( format="%(asctime)s %(name)s %(levelname)s %(message)s", force=True, - level=logginglevel, - filename=f"{os.environ.get('NECTARCHAIN_LOG')}/{os.getpid()}/{Path(__file__).stem}_{os.getpid()}.log", + level=args.verbosity, + filename=f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{os.getpid()}/{Path(__file__).stem}_{os.getpid()}.log", ) log = logging.getLogger(__name__) - log.setLevel(logginglevel) + log.setLevel(args.verbosity) ##tips to add message to stdout handler = logging.StreamHandler(sys.stdout) - handler.setLevel(logginglevel) + handler.setLevel(args.verbosity) formatter = logging.Formatter( "%(asctime)s - %(name)s - %(levelname)s - %(message)s" ) handler.setFormatter(formatter) log.addHandler(handler) - main(args) + kwargs.pop("verbosity") + kwargs.pop("figpath") + kwargs.pop("HHV_run_number") + + log.info(f"arguments passed to main are : {kwargs}") + main(log=log, **kwargs) diff --git a/src/nectarchain/user_scripts/ggrolleron/gain_PhotoStat_computation.sh b/src/nectarchain/user_scripts/ggrolleron/gain_PhotoStat_computation.sh index e6809ed3..5ed2fe8e 100644 --- a/src/nectarchain/user_scripts/ggrolleron/gain_PhotoStat_computation.sh +++ b/src/nectarchain/user_scripts/ggrolleron/gain_PhotoStat_computation.sh @@ -1,5 +1,4 @@ #HOW TO : #to perform photo-statistic high and low gain computation with pedestal run 2609, flat field run 2609 and SPE fit result from run 2634 (1400V run) -python gain_PhotoStat_computation.py -p 2609 -f 2608 --chargeExtractorPath LocalPeakWindowSum_4-12 --FFchargeExtractorWindowLength 16 --correlation --overwrite --SPE_fit_results "/data/users/ggroller/NECTARCAM/SPEfit/data/MULTI-1400V-SPEStd-2634-LocalPeakWindowSum_4-12/output_table.ecsv" --SPE_fit_results_tag VVH2634 +python gain_PhotoStat_computation.py --FF_run_number 3937 --Ped_run_number 3938 --HHV_run_number 3942 --max_events 100 --method LocalPeakWindowSum --extractor_kwargs '{"window_width":16,"window_shift":4}' --overwrite -v INFO --reload_events -python gain_PhotoStat_computation.py -p 3938 -f 3937 --chargeExtractorPath LocalPeakWindowSum_4-12 --FFchargeExtractorWindowLength 16 --correlation --overwrite --SPE_fit_results "/data/users/ggroller/NECTARCAM/SPEfit/data/MULTI-1400V-SPEStd-3942-LocalPeakWindowSum_4-12/results_FlatFieldSingleStdSPE.ecsv" --SPE_fit_results_tag VVH3942 diff --git a/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_combined_computation.py b/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_combined_computation.py index 0e1333a0..a43706bd 100644 --- a/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_combined_computation.py +++ b/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_combined_computation.py @@ -1,176 +1,228 @@ +import json import logging import os import sys import time from pathlib import Path -os.makedirs(os.environ.get("NECTARCHAIN_LOG"), exist_ok=True) - # to quiet numba +os.makedirs(os.environ.get("NECTARCHAIN_LOG"), exist_ok=True) logging.getLogger("numba").setLevel(logging.WARNING) +logging.basicConfig( + format="%(asctime)s %(name)s %(levelname)s %(message)s", + level=logging.DEBUG, + filename=f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{Path(__file__).stem}_{os.getpid()}.log", +) +log = logging.getLogger(__name__) import argparse +import copy +import glob -# import seaborn as sns -from nectarchain.data.container import ChargeContainer -from nectarchain.makers.calibration.gain.flatfield_spe_makers import ( - FlatFieldSingleNominalSPEMaker, +from nectarchain.data.management import DataManagement +from nectarchain.makers.calibration import ( + FlatFieldSPECombinedStdNectarCAMCalibrationTool, ) +from nectarchain.makers.extractor.utils import CtapipeExtractor parser = argparse.ArgumentParser( prog="gain_SPEfit_combined_computation.py", - description="compute high gain with SPE combined fit for one run at nominal voltage. Output data will be saved in $NECTARCAMDATA/../SPEfit/data/{multipath}nominal-prefitCombinedSPE{args.SPE_fit_results_tag}-SPEStd-{args.run_number}-{args.chargeExtractorPath}/", + description=f"compute high gain with SPE fit for one run at nominal voltage from a SPE result from a run at 1400V. By default, output data are saved in $NECTARCAMDATA/../SPEfit/data/", ) - # run numbers -parser.add_argument("-r", "--run_number", help="spe run", type=int) +parser.add_argument( + "-r", "--run_number", nargs="+", default=[], help="run(s) list", type=int +) +# max events to be loaded parser.add_argument( - "--overwrite", - action="store_true", - default=False, - help="to force overwrite files on disk", + "-m", + "--max_events", + nargs="+", + # default=[], + help="max events to be load", + type=int, ) -# output figures and path extension +# boolean arguments parser.add_argument( - "--display", action="store_true", default=False, help="whether to save plot or not" + "--reload_events", + action="store_true", + default=False, + help="to force re-computation of waveforms from fits.fz files", ) parser.add_argument( - "--output_fig_tag", type=str, default="", help="tag to set output figure path" + "--overwrite", + action="store_true", + default=False, + help="to force overwrite files on disk", ) -# pixels selected parser.add_argument( - "-p", "--pixels", nargs="+", default=None, help="pixels selected", type=int + "--events_per_slice", + type=int, + default=None, + help="will slplit the raw data in fits.fz file with events_per_slice events per slices", ) - -# multiprocessing args +# extractor arguments parser.add_argument( - "--multiproc", action="store_true", default=False, help="to use multiprocessing" + "--method", + choices=[ + "FullWaveformSum", + "FixedWindowSum", + "GlobalPeakWindowSum", + "LocalPeakWindowSum", + "SlidingWindowMaxSum", + "TwoPassWindowSum", + ], + default="LocalPeakWindowSum", + help="charge extractor method", + type=str, ) -parser.add_argument("--nproc", help="nproc used for multiprocessing", type=int) parser.add_argument( - "--chunksize", - help="chunksize used for multiprocessing, with multiproccesing, create one process per pixels is not optimal, we rather prefer to group quite a few pixels in same process, chunksize is used to set the number of pixels we use for one process, for example if you want to perform the gain computation of the whole camera with 1855 pixels on 6 CPU, a chunksize of 20 seems to be quite optimal ", - type=int, + "--extractor_kwargs", + default={"window_width": 16, "window_shift": 4}, + help="charge extractor kwargs", + type=json.loads, ) - -# extractor arguments +# verbosity argument parser.add_argument( - "--chargeExtractorPath", - help="charge extractor path where charges are saved", + "-v", + "--verbosity", + help="set the verbosity level of logger", + default="INFO", + choices=["DEBUG", "INFO", "WARN", "ERROR", "CRITICAL"], type=str, ) -# for VVH combined fit +# output figures path parser.add_argument( - "--combined", + "--display", action="store_true", default=False, - help="if True : perform a combined fit of VVH and nominal data, if False : perform a nominal fit with SPE resoltion fixed from VVH fitted data", + help="to plot the SPE histogram for each pixel", ) parser.add_argument( - "--VVH_fitted_results", - help="previoulsy fitted VVH data path for nominal SPE fit by fixing some shared parameters", + "--figpath", type=str, + default=f"{os.environ.get('NECTARCHAIN_FIGURES','/tmp')}/", + help="output figure path", ) -# tag for SPE fit results propagation + +# pixels selected parser.add_argument( - "--SPE_fit_results_tag", - help="SPE fit results tag for propagate the SPE result to output, this tag will be used to setup the path where output data will be saved, see help for description", - type=str, - default="", + "-p", "--asked_pixels_id", nargs="+", default=None, help="pixels selected", type=int +) + +# multiprocessing args +parser.add_argument( + "--multiproc", action="store_true", default=False, help="to use multiprocessing" ) +parser.add_argument( + "--nproc", help="nproc used for multiprocessing", default=8, type=int +) +parser.add_argument( + "--chunksize", help="chunksize used for multiprocessing", default=1, type=int +) + +# combined fit parser.add_argument( "--same_luminosity", action="store_true", default=False, help="if luminosity for VVH and nominal data is the same", ) - -# verbosity argument parser.add_argument( - "-v", - "--verbosity", - help="set the verbosity level of logger", - default="info", - choices=["fatal", "debug", "info", "warning"], - type=str, + "--HHV_run_number", + help="HHV run number of which the SPE fit has ever been performed", + type=int, ) +args = parser.parse_args() -def main(args): - figpath = f"{os.environ.get('NECTARCHAIN_FIGURES',f'/tmp/nectarchain_log/{os.getpid()}/figure')}/" - figpath_ext = "" if args.output_fig_tag == "" else f"-{args.output_fig_tag}" - multipath = "MULTI-" if args.multiproc else "" +def main( + log, + **kwargs, +): + run_number = kwargs.pop("run_number") + max_events = kwargs.pop("max_events", [None for i in range(len(run_number))]) + if max_events is None: + max_events = [None for i in range(len(run_number))] - charge_run = ChargeContainer.from_file( - f"{os.environ.get('NECTARCAMDATA')}/charges/{args.chargeExtractorPath}/", - args.run_number, - ) + log.info(f"max_events : {max_events}") - if args.combined: - raise NotImplementedError("combined fit not implemented yet") - else: - gain_Std = FlatFieldSingleNominalSPEMaker.create_from_chargeContainer( - signal=charge_run, - nectarGainSPEresult=args.VVH_fitted_results, - same_luminosity=args.same_luminosity, - ) - t = time.time() - gain_Std.make( - pixels_id=args.pixels, - multiproc=args.multiproc, - display=args.display, - nproc=args.nproc, - chunksize=args.chunksize, - figpath=figpath - + f"/{multipath}nominal-prefitCombinedSPE{args.SPE_fit_results_tag}-SPEStd-{args.run_number}-{args.chargeExtractorPath}{figpath_ext}", - ) + figpath = args.figpath - log.info(f"fit time = {time.time() - t } sec") - gain_Std.save( - f"{os.environ.get('NECTARCAMDATA')}/../SPEfit/data/{multipath}nominal-prefitCombinedSPE{args.SPE_fit_results_tag}-SPEStd-{args.run_number}-{args.chargeExtractorPath}/", - overwrite=args.overwrite, - ) + str_extractor_kwargs = CtapipeExtractor.get_extractor_kwargs_str( + args.extractor_kwargs + ) + path = DataManagement.find_SPE_HHV( + run_number=args.HHV_run_number, + method=args.method, + str_extractor_kwargs=str_extractor_kwargs, + ) + if len(path) == 1: log.info( - f"convergence rate : {len(gain_Std._results[gain_Std._results['is_valid']])/gain_Std.npixels}" + f"{path[0]} found associated to HHV run {args.HHV_run_number}, method {args.method} and extractor kwargs {str_extractor_kwargs}" ) + else: + _text = f"no file found in $NECTARCAM_DATA/../SPEfit associated to HHV run {args.HHV_run_number}, method {args.method} and extractor kwargs {str_extractor_kwargs}" + log.error(_text) + raise FileNotFoundError(_text) + for _run_number, _max_events in zip(run_number, max_events): + try: + tool = FlatFieldSPECombinedStdNectarCAMCalibrationTool( + progress_bar=True, + run_number=_run_number, + max_events=_max_events, + SPE_result=path[0], + **kwargs, + ) + tool.setup() + tool.start() + if args.reload_events and not (_max_events is None): + _figpath = f"{figpath}/{tool.name}_run{tool.run_number}_maxevents{_max_events}_{tool.method}_{str_extractor_kwargs}" + else: + _figpath = f"{figpath}/{tool.name}_run{tool.run_number}_{tool.method}_{str_extractor_kwargs}" + tool.finish(figpath=_figpath, display=args.display) + except Exception as e: + log.warning(e, exc_info=True) if __name__ == "__main__": + t = time.time() args = parser.parse_args() - logginglevel = logging.FATAL - if args.verbosity == "warning": - logginglevel = logging.WARNING - elif args.verbosity == "info": - logginglevel = logging.INFO - elif args.verbosity == "debug": - logginglevel = logging.DEBUG - - os.makedirs( - f"{os.environ.get('NECTARCHAIN_LOG','/tmp/nectarchain_log')}/{os.getpid()}/figures" - ) + kwargs = copy.deepcopy(vars(args)) + + kwargs["log_level"] = args.verbosity + + os.makedirs(f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{os.getpid()}/figures") logging.basicConfig( format="%(asctime)s %(name)s %(levelname)s %(message)s", force=True, - level=logginglevel, - filename=f"{os.environ.get('NECTARCHAIN_LOG')}/{os.getpid()}/{Path(__file__).stem}_{os.getpid()}.log", + level=args.verbosity, + filename=f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{os.getpid()}/{Path(__file__).stem}_{os.getpid()}.log", ) log = logging.getLogger(__name__) - log.setLevel(logginglevel) + log.setLevel(args.verbosity) ##tips to add message to stdout handler = logging.StreamHandler(sys.stdout) - handler.setLevel(logginglevel) + handler.setLevel(args.verbosity) formatter = logging.Formatter( "%(asctime)s - %(name)s - %(levelname)s - %(message)s" ) handler.setFormatter(formatter) log.addHandler(handler) - main(args) + kwargs.pop("verbosity") + kwargs.pop("figpath") + kwargs.pop("display") + kwargs.pop("HHV_run_number") + + log.info(f"arguments passed to main are : {kwargs}") + main(log=log, **kwargs) + log.info(f"time for execution is {time.time() - t:.2e} sec") diff --git a/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_combined_computation.sh b/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_combined_computation.sh index 609a4138..b1ae5be8 100644 --- a/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_combined_computation.sh +++ b/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_combined_computation.sh @@ -1,9 +1,4 @@ #HOW TO : -#to perform SPE fit of run 2633 (supposed to be taken at nominal voltage) from SPE fit performed on 2634 (taken at 1400V) -python gain_SPEfit_combined_computation.py -r 2633 --display --SPE_fit_results_tag 2634 --chargeExtractorPath LocalPeakWindowSum_4-12 --overwrite --multiproc --nproc 6 --chunksize 50 --same_luminosity --VVH_fitted_results "/data/users/ggroller/NECTARCAM/SPEfit/data/MULTI-1400V-SPEStd-2634-LocalPeakWindowSum_4-12/output_table.ecsv" - -#to perform a joint SPE fit of run 2633 (supposed to be taken at nominal voltage) and run 2634 (1400V). -python gain_SPEfit_combined_computation.py -r 2633 --display --SPE_fit_results_tag 2634 --chargeExtractorPath LocalPeakWindowSum_4-12 --overwrite --multiproc --nproc 6 --chunksize 50 --same_luminosity --combined #to perform SPE fit of run 3936 (supposed to be taken at nominal voltage) from SPE fit performed on 3942 (taken at 1400V) (luminosity can not be supposed the same) -python gain_SPEfit_combined_computation.py -r 3936 --display --SPE_fit_results_tag 3942 --chargeExtractorPath LocalPeakWindowSum_4-12 --overwrite --multiproc --nproc 6 --chunksize 50 --VVH_fitted_results "/data/users/ggroller/NECTARCAM/SPEfit/data/MULTI-1400V-SPEStd-3942-LocalPeakWindowSum_4-12/output_table.ecsv" +python gain_SPEfit_combined_computation.py -r 3936 --HHV_run_number 3942 --multiproc --nproc 8 --chunksize 2 -p 45 56 --max_events 100 --method LocalPeakWindowSum --extractor_kwargs '{"window_width":16,"window_shift":4}' --overwrite -v DEBUG --reload_events diff --git a/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_computation.py b/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_computation.py index c58e0cc2..0a1cc23f 100644 --- a/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_computation.py +++ b/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_computation.py @@ -1,32 +1,57 @@ +import json import logging import os import sys import time from pathlib import Path -os.makedirs(os.environ.get("NECTARCHAIN_LOG"), exist_ok=True) - # to quiet numba +os.makedirs(os.environ.get("NECTARCHAIN_LOG"), exist_ok=True) logging.getLogger("numba").setLevel(logging.WARNING) - +logging.basicConfig( + format="%(asctime)s %(name)s %(levelname)s %(message)s", + level=logging.DEBUG, + filename=f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{Path(__file__).stem}_{os.getpid()}.log", +) +log = logging.getLogger(__name__) import argparse +import copy -# import seaborn as sns -from nectarchain.data.container import ChargeContainer -from nectarchain.makers.calibration.gain.flatfield_spe_makers import ( - FlatFieldSingleHHVSPEMaker, - FlatFieldSingleHHVStdSPEMaker, +from nectarchain.makers.calibration import ( + FlatFieldSPEHHVNectarCAMCalibrationTool, + FlatFieldSPEHHVStdNectarCAMCalibrationTool, + FlatFieldSPENominalNectarCAMCalibrationTool, + FlatFieldSPENominalStdNectarCAMCalibrationTool, ) +from nectarchain.makers.extractor.utils import CtapipeExtractor parser = argparse.ArgumentParser( prog="gain_SPEfit_computation.py", - description="compute high gain with SPE fit for one run at very very high voltage (~1400V) or at nominal voltage (it can often fail). Output data are saved in $NECTARCAMDATA/../SPEfit/data/{multipath}{args.voltage_tag}-{SPEpath}-{args.run_number}-{args.chargeExtractorPath}/", + description=f"compute high gain with SPE fit for one run at very very high voltage (~1400V) or at nominal voltage (it can often fail). By default, output data are saved in $NECTARCAMDATA/../SPEfit/data/", ) - # run numbers -parser.add_argument("-r", "--run_number", help="spe run", type=int) +parser.add_argument( + "-r", "--run_number", nargs="+", default=[], help="run(s) list", type=int +) +# max events to be loaded +parser.add_argument( + "-m", + "--max_events", + nargs="+", + # default=[], + help="max events to be load", + type=int, +) + +# boolean arguments +parser.add_argument( + "--reload_events", + action="store_true", + default=False, + help="to force re-computation of waveforms from fits.fz files", +) parser.add_argument( "--overwrite", action="store_true", @@ -35,23 +60,61 @@ ) parser.add_argument( - "--voltage_tag", + "--events_per_slice", + type=int, + default=None, + help="will slplit the raw data in fits.fz file with events_per_slice events per slices", +) + +# extractor arguments +parser.add_argument( + "--method", + choices=[ + "FullWaveformSum", + "FixedWindowSum", + "GlobalPeakWindowSum", + "LocalPeakWindowSum", + "SlidingWindowMaxSum", + "TwoPassWindowSum", + ], + default="LocalPeakWindowSum", + help="charge extractor method", type=str, - default="", - help="tag for voltage specifcication (1400V or nominal), used to setup the output path. See help for more details", +) +parser.add_argument( + "--extractor_kwargs", + default={"window_width": 16, "window_shift": 4}, + help="charge extractor kwargs", + type=json.loads, ) -# output figures and path extension +# verbosity argument parser.add_argument( - "--display", action="store_true", default=False, help="whether to save plot or not" + "-v", + "--verbosity", + help="set the verbosity level of logger", + default="INFO", + choices=["DEBUG", "INFO", "WARN", "ERROR", "CRITICAL"], + type=str, ) + +# output figures path parser.add_argument( - "--output_fig_tag", type=str, default="", help="tag to set output figure path" + "--display", + action="store_true", + default=False, + help="to plot the SPE histogram for each pixel", +) +parser.add_argument( + "--figpath", + type=str, + default=f"{os.environ.get('NECTARCHAIN_FIGURES','/tmp')}/", + help="output figure path", ) # pixels selected parser.add_argument( - "-p", "--pixels", nargs="+", default=None, help="pixels selected", type=int + "-p", "--asked_pixels_id", nargs="+", default=None, help="pixels selected", type=int ) # for let free pp and n : @@ -59,104 +122,102 @@ "--free_pp_n", action="store_true", default=False, help="to let free pp and n" ) +# to say if it is a run taken at HHV +parser.add_argument( + "--HHV", + action="store_true", + default=False, + help="to say that these runs are taken at HHV, it will change the configuration file used", +) + # multiprocessing args parser.add_argument( "--multiproc", action="store_true", default=False, help="to use multiprocessing" ) -parser.add_argument("--nproc", help="nproc used for multiprocessing", type=int) -parser.add_argument("--chunksize", help="chunksize used for multiprocessing", type=int) - -# extractor arguments parser.add_argument( - "--chargeExtractorPath", - help="charge extractor path where charges are saved", - type=str, + "--nproc", help="nproc used for multiprocessing", default=8, type=int ) - -# verbosity argument parser.add_argument( - "-v", - "--verbosity", - help="set the verbosity level of logger", - default="info", - choices=["fatal", "debug", "info", "warning"], - type=str, + "--chunksize", help="chunksize used for multiprocessing", default=1, type=int ) +args = parser.parse_args() -def main(args): - figpath = f"{os.environ.get('NECTARCHAIN_FIGURES')}/" - figpath_ext = "" if args.output_fig_tag == "" else f"-{args.output_fig_tag}" +def main( + log, + **kwargs, +): + run_number = kwargs.pop("run_number") + max_events = kwargs.pop("max_events", [None for i in range(len(run_number))]) + if max_events is None: + max_events = [None for i in range(len(run_number))] - multipath = "MULTI-" if args.multiproc else "" - SPEpath = "SPE" if args.free_pp_n else "SPEStd" + log.info(f"max_events : {max_events}") - charge_run_1400V = ChargeContainer.from_file( - f"{os.environ.get('NECTARCAMDATA')}/charges/{args.chargeExtractorPath}/", - args.run_number, - ) - - if args.free_pp_n: - gain_Std = FlatFieldSingleHHVSPEMaker.create_from_chargeContainer( - signal=charge_run_1400V - ) + figpath = args.figpath + if args.HHV: + if args.free_pp_n: + _class = FlatFieldSPEHHVNectarCAMCalibrationTool + else: + _class = FlatFieldSPEHHVStdNectarCAMCalibrationTool else: - gain_Std = FlatFieldSingleHHVStdSPEMaker.create_from_chargeContainer( - signal=charge_run_1400V - ) + if args.free_pp_n: + _class = FlatFieldSPENominalNectarCAMCalibrationTool + else: + _class = FlatFieldSPENominalStdNectarCAMCalibrationTool + + for _run_number, _max_events in zip(run_number, max_events): + try: + tool = _class( + progress_bar=True, + run_number=_run_number, + max_events=_max_events, + **kwargs, + ) + tool.setup() + if args.reload_events and not (_max_events is None): + _figpath = f"{figpath}/{tool.name}_run{tool.run_number}_maxevents{_max_events}_{tool.method}_{CtapipeExtractor.get_extractor_kwargs_str(tool.extractor_kwargs)}" + else: + _figpath = f"{figpath}/{tool.name}_run{tool.run_number}_{tool.method}_{CtapipeExtractor.get_extractor_kwargs_str(tool.extractor_kwargs)}" + tool.start(figpath=_figpath, display=args.display) + tool.finish(figpath=_figpath, display=args.display) + except Exception as e: + log.warning(e, exc_info=True) + +if __name__ == "__main__": t = time.time() - gain_Std.make( - pixels_id=args.pixels, - multiproc=args.multiproc, - display=args.display, - nproc=args.nproc, - chunksize=args.chunksize, - figpath=figpath - + f"/{multipath}{args.voltage_tag}-{SPEpath}-{args.run_number}-{args.chargeExtractorPath}{figpath_ext}", - ) - log.info(f"fit time = {time.time() - t } sec") - gain_Std.save( - f"{os.environ.get('NECTARCAMDATA')}/../SPEfit/data/{multipath}{args.voltage_tag}-{SPEpath}-{args.run_number}-{args.chargeExtractorPath}/", - overwrite=args.overwrite, - ) - conv_rate = ( - len(gain_Std._results[gain_Std._results["is_valid"]]) / gain_Std.npixels - if args.pixels is None - else len(gain_Std._results[gain_Std._results["is_valid"]]) / len(args.pixels) - ) - log.info(f"convergence rate : {conv_rate}") + args = parser.parse_args() + kwargs = copy.deepcopy(vars(args)) + kwargs["log_level"] = args.verbosity -if __name__ == "__main__": - args = parser.parse_args() - logginglevel = logging.FATAL - if args.verbosity == "warning": - logginglevel = logging.WARNING - elif args.verbosity == "info": - logginglevel = logging.INFO - elif args.verbosity == "debug": - logginglevel = logging.DEBUG - - os.makedirs(f"{os.environ.get('NECTARCHAIN_LOG')}/{os.getpid()}/figures") + os.makedirs(f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{os.getpid()}/figures") logging.basicConfig( format="%(asctime)s %(name)s %(levelname)s %(message)s", force=True, - level=logginglevel, - filename=f"{os.environ.get('NECTARCHAIN_LOG')}/{os.getpid()}/{Path(__file__).stem}_{os.getpid()}.log", + level=args.verbosity, + filename=f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{os.getpid()}/{Path(__file__).stem}_{os.getpid()}.log", ) log = logging.getLogger(__name__) - log.setLevel(logginglevel) + log.setLevel(args.verbosity) ##tips to add message to stdout handler = logging.StreamHandler(sys.stdout) - handler.setLevel(logginglevel) + handler.setLevel(args.verbosity) formatter = logging.Formatter( "%(asctime)s - %(name)s - %(levelname)s - %(message)s" ) handler.setFormatter(formatter) log.addHandler(handler) - main(args) + kwargs.pop("verbosity") + kwargs.pop("figpath") + kwargs.pop("display") + kwargs.pop("HHV") + + log.info(f"arguments passed to main are : {kwargs}") + main(log=log, **kwargs) + log.info(f"time for execution is {time.time() - t:.2e} sec") diff --git a/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_computation.sh b/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_computation.sh index 932d6721..8ca3fbea 100644 --- a/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_computation.sh +++ b/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_computation.sh @@ -1,13 +1,11 @@ #HOW TO : -#to perform SPE fit of a HHV run -python gain_SPEfit_computation.py -r 2634 --display --voltage_tag "1400V" --chargeExtractorPath LocalPeakWindowSum_4-12 --overwrite --multiproc --nproc 6 --chunksize 2 -p 0 1 2 3 4 5 6 7 8 9 10 11 12 - -#to perform SPE fit of a HHV run letting n and pp parameters free -python gain_SPEfit_computation.py -r 2634 --display --voltage_tag "1400V" --free_pp_n --chargeExtractorPath LocalPeakWindowSum_4-12 --overwrite --multiproc --nproc 6 --chunksize 2 -p 0 1 2 3 4 5 6 7 8 9 10 11 12 +#HVV runs : 2634,3942 +#nominal runs : 2633,3936 #to perform SPE fit of a run at nominal voltage -python gain_SPEfit_computation.py -r 3936 --display --voltage_tag "nominal" --chargeExtractorPath LocalPeakWindowSum_4-12 --overwrite --multiproc --nproc 6 --chunksize 50 - +python gain_SPEfit_computation.py -r 3942 --HHV --reload_events --multiproc --nproc 8 --chunksize 2 -p 45 56 --max_events 100 --method LocalPeakWindowSum --extractor_kwargs '{"window_width":16,"window_shift":4}' --overwrite -v DEBUG #to perform SPE fit of a HHV run -python gain_SPEfit_computation.py -r 3942 --display --output_fig_tag "testMyRes" --voltage_tag "1400V" --chargeExtractorPath LocalPeakWindowSum_4-12 --multiproc --nproc 6 --chunksize 2 -p 0 1 2 3 4 5 6 7 8 9 10 11 12 --verbosity info --overwrite \ No newline at end of file +python gain_SPEfit_computation.py -r 3936 --reload_events --multiproc --nproc 8 --chunksize 2 -p 45 56 --max_events 100 --method LocalPeakWindowSum --extractor_kwargs '{"window_width":16,"window_shift":4}' --overwrite -v DEBUG +#to perform SPE fit of a HHV run letting n and pp parameters free +python gain_SPEfit_computation.py -r 3942 --HHV --free_pp_n --reload_events --multiproc --nproc 8 --chunksize 2 -p 45 56 --max_events 100 --method LocalPeakWindowSum --extractor_kwargs '{"window_width":16,"window_shift":4}' --overwrite -v DEBUG diff --git a/src/nectarchain/user_scripts/ggrolleron/load_wfs_compute_charge.py b/src/nectarchain/user_scripts/ggrolleron/load_wfs_compute_charge.py index 24249618..d6e9939d 100644 --- a/src/nectarchain/user_scripts/ggrolleron/load_wfs_compute_charge.py +++ b/src/nectarchain/user_scripts/ggrolleron/load_wfs_compute_charge.py @@ -1,25 +1,30 @@ import argparse -import glob import json import logging import os import sys from pathlib import Path +os.makedirs(os.environ.get("NECTARCHAIN_LOG"), exist_ok=True) logging.getLogger("numba").setLevel(logging.WARNING) logging.basicConfig( format="%(asctime)s %(name)s %(levelname)s %(message)s", level=logging.DEBUG, - filename=f"{os.environ.get('NECTARCHAIN_LOG')}/{Path(__file__).stem}_{os.getpid()}.log", + filename=f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{Path(__file__).stem}_{os.getpid()}.log", ) log = logging.getLogger(__name__) +import copy from nectarchain.data.container import ( - ChargeContainer, - ChargeContainers, + ChargesContainer, + ChargesContainers, WaveformsContainer, WaveformsContainers, ) +from nectarchain.makers import ( + ChargesNectarCAMCalibrationTool, + WaveformsNectarCAMCalibrationTool, +) parser = argparse.ArgumentParser( prog="load_wfs_compute_charge", @@ -28,58 +33,16 @@ # run numbers parser.add_argument( - "-s", "--spe_run_number", nargs="+", default=[], help="spe run list", type=int -) -parser.add_argument( - "-p", "--ped_run_number", nargs="+", default=[], help="ped run list", type=int -) -parser.add_argument( - "-f", "--ff_run_number", nargs="+", default=[], help="FF run list", type=int + "-r", "--run_number", nargs="+", default=[], help="run(s) list", type=int ) # max events to be loaded parser.add_argument( - "--spe_max_events", - nargs="+", - # default=[], - help="spe max events to be load", - type=int, -) -parser.add_argument( - "--ped_max_events", - nargs="+", - # default=[], - help="ped max events to be load", - type=int, -) -parser.add_argument( - "--ff_max_events", - nargs="+", - # default=[], - help="FF max events to be load", - type=int, -) - -# n_events in runs -parser.add_argument( - "--spe_nevents", + "-m", + "--max_events", nargs="+", # default=[], - help="spe n events to be load", - type=int, -) -parser.add_argument( - "--ped_nevents", - nargs="+", - # default=[], - help="ped n events to be load", - type=int, -) -parser.add_argument( - "--ff_nevents", - nargs="+", - # default=[], - help="FF n events to be load", + help="max events to be load", type=int, ) @@ -98,15 +61,15 @@ ) parser.add_argument( - "--split", - action="store_true", - default=False, - help="split waveforms extraction with 1 file per fits.fz raw data file", + "--events_per_slice", + type=int, + default=None, + help="will split the raw data in fits.fz file with events_per_slice events per slices", ) # extractor arguments parser.add_argument( - "--extractorMethod", + "--method", choices=[ "FullWaveformSum", "FixedWindowSum", @@ -131,269 +94,80 @@ "-v", "--verbosity", help="set the verbosity level of logger", - default="info", - choices=["fatal", "debug", "info", "warning"], + default="INFO", + choices=["DEBUG", "INFO", "WARN", "ERROR", "CRITICAL"], type=str, ) args = parser.parse_args() -# control shape of arguments lists -for arg in ["spe", "ff", "ped"]: - run_number = eval(f"args.{arg}_run_number") - max_events = eval(f"args.{arg}_max_events") - nevents = eval(f"args.{arg}_nevents") - - if not (max_events is None) and len(max_events) != len(run_number): - e = Exception(f"{arg}_run_number and {arg}_max_events must have same length") - log.error(e, exc_info=True) - raise e - if not (nevents is None) and len(nevents) != len(run_number): - e = Exception(f"{arg}_run_number and {arg}_nevents must have same length") - log.error(e, exc_info=True) - raise e - - -def load_wfs_no_split(i, runs_list, max_events, nevents, overwrite): - """method to load waveforms without splitting - - Args: - i (int): index in the run list - runs_list (list): the run number list - max_events (list): max_events list - nevents (list): nevents list - overwrite (bool): to overwrite - - Returns: - WaveformsContainer: the output waveformsContainer - """ - log.info("loading wfs not splitted") - wfs = WaveformsContainer(runs_list[i], max_events=max_events[i], nevents=nevents[i]) - wfs.load_wfs() - wfs.write(f"{os.environ['NECTARCAMDATA']}/waveforms/", overwrite=overwrite) - return wfs - - -def load_wfs_charge_split( - i, runs_list, max_events, overwrite, charge_childpath, extractor_kwargs -): - """_summary_ - - Args: - i (int): index in the run list - runs_list (list): the run number list - max_events (list): max_events list - nevents (list): nevents list - overwrite (bool): to overwrite - charge_childpath (str): the extraction method - extractor_kwargs (dict): the charge extractor kwargs - - Returns: - WaveformsContainers,ChargeContainers: the output WaveformsContainers and ChargeContainers - """ - - log.info("splitting wafevorms extraction with raw data list files") - log.debug(f"creation of the WaveformsContainers") - wfs = WaveformsContainers(runs_list[i], max_events=max_events[i], init_arrays=False) - log.info(f"computation of charge with {charge_childpath}") - log.info("splitting charge computation with raw data list files") - charge = ChargeContainers() - for j in range(wfs.nWaveformsContainer): - log.debug(f"reader events for file {j}") - wfs.load_wfs(index=j) - wfs.write( - f"{os.environ['NECTARCAMDATA']}/waveforms/", index=j, overwrite=overwrite - ) - log.debug(f"computation of charge for file {j}") - charge.append( - ChargeContainer.from_waveforms( - wfs.waveformsContainer[j], method=charge_childpath, **extractor_kwargs - ) - ) - log.debug(f"deleting waveformsContainer at index {j} to free RAM") - wfs.waveformsContainer[j] = WaveformsContainer.__new__(WaveformsContainer) - - log.info("merging charge") - charge = charge.merge() - return wfs, charge - - -def load_wfs_charge_split_from_wfsFiles(wfsFiles, charge_childpath, extractor_kwargs): - """_summary_ - Args: - wfsFiles (list): list of the waveformsContainer FITS files - charge_childpath (str): the extraction method - extractor_kwargs (dict): the charge extractor kwargs - - Returns: - None,ChargeContainers: the output ChargeContainers (return tuple with None to keep same structure as load_wfs_charge_split) - """ - charge = ChargeContainers() - for j, file in enumerate(wfsFiles): - log.debug(f"loading wfs from file {file}") - wfs = WaveformsContainer.load(file) - log.debug(f"computation of charge for file {file}") - charge.append( - ChargeContainer.from_waveforms( - wfs, method=charge_childpath, **extractor_kwargs - ) - ) - log.debug(f"deleting waveformsContainer from {file} to free RAM") - del wfs.wfs_hg - del wfs.wfs_lg - del wfs.ucts_timestamp - del wfs.ucts_busy_counter - del wfs.ucts_event_counter - del wfs.event_type - del wfs.event_id - del wfs.trig_pattern_all - del wfs - # gc.collect() - - log.info("merging charge") - charge = charge.merge() - return None, charge - - -def load_wfs_compute_charge( - runs_list: list, - reload_wfs: bool = False, - overwrite: bool = False, - charge_extraction_method: str = "FullWaveformSum", +def main( + log, **kwargs, -) -> None: - """this method is used to load waveforms from zfits files and compute charge with an user specified method - - Args: - runs_list (list): list of runs for which you want to perfrom waveforms and charge extraction - reload_wfs (bool, optional): argument used to reload waveforms from pre-loaded waveforms (in fits format) or from zfits file. Defaults to False. - overwrite (bool, optional): to overwrite file on disk. Defaults to False. - charge_extraction_method (str, optional): ctapipe charge extractor. Defaults to "FullWaveformSum". - - Raises: - e : an error occurred during zfits loading from ctapipe EventSource - """ - - # print(runs_list) - # print(charge_extraction_method) - # print(overwrite) - # print(reload_wfs) +): # print(kwargs) - - max_events = kwargs.get("max_events", [None for i in range(len(runs_list))]) - nevents = kwargs.get("nevents", [-1 for i in range(len(runs_list))]) - - charge_childpath = kwargs.get("charge_childpath", charge_extraction_method) - extractor_kwargs = kwargs.get("extractor_kwargs", {}) - - split = kwargs.get("split", False) - - for i in range(len(runs_list)): - log.info(f"treating run {runs_list[i]}") - log.info("waveform computation") - if not (reload_wfs): - log.info( - f"trying to load waveforms from {os.environ['NECTARCAMDATA']}/waveforms/" - ) - try: - if split: - files = glob.glob( - f"{os.environ['NECTARCAMDATA']}/waveforms/waveforms_run{runs_list[i]}_*.fits" - ) - if len(files) == 0: - raise FileNotFoundError(f"no splitted waveforms found") - else: - wfs, charge = load_wfs_charge_split_from_wfsFiles( - files, charge_childpath, extractor_kwargs - ) - - else: - wfs = WaveformsContainer.load( - f"{os.environ['NECTARCAMDATA']}/waveforms/waveforms_run{runs_list[i]}.fits" - ) - except FileNotFoundError as e: - log.warning( - f"argument said to not reload waveforms from zfits files but computed waveforms not found at {os.environ['NECTARCAMDATA']}/waveforms/waveforms_run{runs_list[i]}.fits" + run_number = kwargs.pop("run_number") + max_events = kwargs.pop("max_events", [None for i in range(len(run_number))]) + if max_events is None: + max_events = [None for i in range(len(run_number))] + + log.info(max_events) + + tool = WaveformsNectarCAMCalibrationTool() + waveforms_kwargs = {} + for key in tool.traits().keys(): + if key in kwargs.keys(): + waveforms_kwargs[key] = kwargs[key] + + tool = ChargesNectarCAMCalibrationTool() + charges_kwargs = {} + for key in tool.traits().keys(): + if key in kwargs.keys(): + charges_kwargs[key] = kwargs[key] + + log.info(f"WaveformsNectarCAMCalibrationTool kwargs are {waveforms_kwargs}") + log.info(f"ChargesNectarCAMCalibrationTool kwargs are {charges_kwargs}") + + for _run_number, _max_events in zip(run_number, max_events): + try: + if kwargs.get("reload_wfs", False): + log.info("reloading waveforms") + tool = WaveformsNectarCAMCalibrationTool( + progress_bar=True, + run_number=_run_number, + max_events=_max_events, + **waveforms_kwargs, ) - log.warning(f"reloading from zfits files") - if split: - wfs, charge = load_wfs_charge_split( - i, - runs_list, - max_events, - overwrite, - charge_childpath, - extractor_kwargs, - ) - else: - wfs = load_wfs_no_split( - i, runs_list, max_events, nevents, overwrite - ) - except Exception as e: - log.error(e, exc_info=True) - raise e - else: - if split: - wfs, charge = load_wfs_charge_split( - i, - runs_list, - max_events, - overwrite, - charge_childpath, - extractor_kwargs, + tool.setup() + tool.start() + tool.finish() + + tool = ChargesNectarCAMCalibrationTool( + progress_bar=True, + run_number=_run_number, + max_events=_max_events, + from_computed_waveforms=True, + **charges_kwargs, ) + tool.setup() + tool.start() + tool.finish() else: - wfs = load_wfs_no_split(i, runs_list, max_events, nevents, overwrite) - - if not (split): - log.info(f"computation of charge with {charge_childpath}") - charge = ChargeContainer.from_waveforms( - wfs, method=charge_childpath, **extractor_kwargs - ) - del wfs - - charge.write( - f"{os.environ['NECTARCAMDATA']}/charges/{path}/", overwrite=overwrite - ) - del charge - - -def main( - spe_run_number: list = [], - ff_run_number: list = [], - ped_run_number: list = [], - **kwargs, -): - # print(kwargs) - - spe_nevents = kwargs.pop("spe_nevents", [-1 for i in range(len(spe_run_number))]) - ff_nevents = kwargs.pop("ff_nevents", [-1 for i in range(len(ff_run_number))]) - ped_nevents = kwargs.pop("ped_nevents", [-1 for i in range(len(ped_run_number))]) - - spe_max_events = kwargs.pop( - "spe_max_events", [None for i in range(len(spe_run_number))] - ) - ff_max_events = kwargs.pop( - "ff_max_events", [None for i in range(len(ff_run_number))] - ) - ped_max_events = kwargs.pop( - "ped_max_events", [None for i in range(len(ped_run_number))] - ) - - runs_list = spe_run_number + ff_run_number + ped_run_number - nevents = spe_nevents + ff_nevents + ped_nevents - max_events = spe_max_events + ff_max_events + ped_max_events - - charge_extraction_method = kwargs.get("extractorMethod", "FullWaveformSum") + log.info("trying to compute charges from waveforms yet extracted") + tool = ChargesNectarCAMCalibrationTool( + progress_bar=True, + run_number=_run_number, + max_events=_max_events, + from_computed_waveforms=True, + **charges_kwargs, + ) - load_wfs_compute_charge( - runs_list=runs_list, - charge_extraction_method=charge_extraction_method, - nevents=nevents, - max_events=max_events, - **kwargs, - ) + tool.setup() + tool.start() + tool.finish() + except Exception as e: + log.error(e, exc_info=True) if __name__ == "__main__": @@ -404,54 +178,31 @@ def main( # spe_nevents = [49227,49148,-1] args = parser.parse_args() - logginglevel = logging.FATAL - if args.verbosity == "warning": - logginglevel = logging.WARNING - elif args.verbosity == "info": - logginglevel = logging.INFO - elif args.verbosity == "debug": - logginglevel = logging.DEBUG + kwargs = copy.deepcopy(vars(args)) - os.makedirs(f"{os.environ.get('NECTARCHAIN_LOG')}/{os.getpid()}/figures") + kwargs["log_level"] = args.verbosity + + os.makedirs(f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{os.getpid()}/figures") logging.basicConfig( format="%(asctime)s %(name)s %(levelname)s %(message)s", force=True, - level=logginglevel, - filename=f"{os.environ.get('NECTARCHAIN_LOG')}/{os.getpid()}/{Path(__file__).stem}_{os.getpid()}.log", + level=args.verbosity, + filename=f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{os.getpid()}/{Path(__file__).stem}_{os.getpid()}.log", ) log = logging.getLogger(__name__) - log.setLevel(logginglevel) + log.setLevel(args.verbosity) ##tips to add message to stdout handler = logging.StreamHandler(sys.stdout) - handler.setLevel(logginglevel) + handler.setLevel(args.verbosity) formatter = logging.Formatter( "%(asctime)s - %(name)s - %(levelname)s - %(message)s" ) handler.setFormatter(formatter) log.addHandler(handler) - arg = vars(args) - log.info(f"arguments are : {arg}") - - key_to_pop = [] - for key in arg.keys(): - if arg[key] is None: - key_to_pop.append(key) - - for key in key_to_pop: - arg.pop(key) - - log.info(f"arguments passed to main are : {arg}") - - path = args.extractorMethod - if args.extractorMethod in ["GlobalPeakWindowSum", "LocalPeakWindowSum"]: - path += f"_{args.extractor_kwargs['window_shift']}-{args.extractor_kwargs['window_width']-args.extractor_kwargs['window_shift']}" - elif args.extractorMethod in ["SlidingWindowMaxSum"]: - path += f"_{args.extractor_kwargs['window_width']}" - elif args.extractorMethod in ["FixedWindowSum"]: - path += f"_{args.extractor_kwargs['peak_index']}_{args.extractor_kwargs['window_shift']}-{args.extractor_kwargs['window_width']-args.extractor_kwargs['window_shift']}" + kwargs.pop("verbosity") - arg["path"] = path + log.info(f"arguments passed to main are : {kwargs}") - main(**arg) + main(log=log, **kwargs) diff --git a/src/nectarchain/user_scripts/ggrolleron/load_wfs_compute_charge.sh b/src/nectarchain/user_scripts/ggrolleron/load_wfs_compute_charge.sh index 89165512..0f050238 100644 --- a/src/nectarchain/user_scripts/ggrolleron/load_wfs_compute_charge.sh +++ b/src/nectarchain/user_scripts/ggrolleron/load_wfs_compute_charge.sh @@ -1,7 +1,19 @@ #HOW TO : -#spe_run_number = [2633,2634,3784] -#ff_run_number = [2608] -#ped_run_number = [2609] -#spe_nevents = [49227,49148,-1] +#FF_spe_nominal : 2633 ,3936 -python load_wfs_compute_charge.py -s 2633 2634 3784 -f 2608 -p 2609 --spe_nevents 49227 49148 -1 --extractorMethod LocalPeakWindowSum --extractor_kwargs '{"window_width":16,"window_shift":4}' \ No newline at end of file +#FF_spe_5deg : 3731 +#FF_spe_-5deg : 3784 +#FF_spe_0deg : 3750 + +#FF_spe_HHV : 2634, 3942 + +#WT_spe_nominal : 4129 + +#ff_run_number : 2608, 3937 +#ped_run_number : 2609, 3938 + +#FF and Ped : 3943 + +python load_wfs_compute_charge.py -s 2633 2634 3784 -f 2608 -p 2609 --spe_nevents 49227 49148 -1 --method LocalPeakWindowSum --extractor_kwargs '{"window_width":16,"window_shift":4}' + +python load_wfs_compute_charge.py -r 3942 --max_events 100 --method GlobalPeakWindowSum --extractor_kwargs '{"window_width":16,"window_shift":4}' --overwrite -v debug \ No newline at end of file diff --git a/src/nectarchain/user_scripts/ggrolleron/notebooks/waveforms_charges_extraction_new_framework.py b/src/nectarchain/user_scripts/ggrolleron/notebooks/waveforms_charges_extraction_new_framework.py index 94b7fa53..5d9f8c11 100644 --- a/src/nectarchain/user_scripts/ggrolleron/notebooks/waveforms_charges_extraction_new_framework.py +++ b/src/nectarchain/user_scripts/ggrolleron/notebooks/waveforms_charges_extraction_new_framework.py @@ -1,10 +1,18 @@ # %% -from nectarchain.makers import ChargesNectarCAMCalibrationTool,WaveformsNectarCAMCalibrationTool -from nectarchain.makers.component import get_valid_component -from nectarchain.data.container import ChargesContainers,ChargesContainer,WaveformsContainer,WaveformsContainers -from ctapipe.io import HDF5TableReader from ctapipe.containers import EventType +from ctapipe.io import HDF5TableReader +from nectarchain.data.container import ( + ChargesContainer, + ChargesContainers, + WaveformsContainer, + WaveformsContainers, +) +from nectarchain.makers import ( + ChargesNectarCAMCalibrationTool, + WaveformsNectarCAMCalibrationTool, +) +from nectarchain.makers.component import get_valid_component # %% get_valid_component() @@ -13,17 +21,20 @@ run_number = 3942 # %% -tool = WaveformsNectarCAMCalibrationTool(progress_bar = True,run_number = run_number,max_events = 500,log_level = 10) +tool = WaveformsNectarCAMCalibrationTool( + progress_bar=True, run_number=run_number, max_events=500, log_level=10 +) tool # %% -tool = ChargesNectarCAMCalibrationTool(progress_bar = True, - method = 'LocalPeakWindowSum', - extractor_kwargs = {"window_width": 16, "window_shift": 4}, - run_number = run_number, - max_events = 5000, - log_level = 10 - ) +tool = ChargesNectarCAMCalibrationTool( + progress_bar=True, + method="LocalPeakWindowSum", + extractor_kwargs={"window_width": 16, "window_shift": 4}, + run_number=run_number, + max_events=5000, + log_level=10, +) tool # %% @@ -43,10 +54,15 @@ trigger_type = EventType.__members__ -with HDF5TableReader(f"/tmp/EventsLoopNectarCAMCalibrationTool_{run_number}.h5") as reader : - for key,trigger in trigger_type.items() : - try : - tableReader = reader.read(table_name = f"/data/WaveformsContainer_0/{trigger.name}", containers = WaveformsContainer) +with HDF5TableReader( + f"/tmp/EventsLoopNectarCAMCalibrationTool_{run_number}.h5" +) as reader: + for key, trigger in trigger_type.items(): + try: + tableReader = reader.read( + table_name=f"/data/WaveformsContainer_0/{trigger.name}", + containers=WaveformsContainer, + ) container.containers[trigger] = next(tableReader) except Exception as err: print(err) @@ -59,10 +75,15 @@ trigger_type = EventType.__members__ -with HDF5TableReader(f"/tmp/EventsLoopNectarCAMCalibrationTool_{run_number}.h5") as reader : - for key,trigger in trigger_type.items() : - try : - tableReader = reader.read(table_name = f"/data/ChargesContainer_0/{trigger.name}", containers = ChargesContainer) +with HDF5TableReader( + f"/tmp/EventsLoopNectarCAMCalibrationTool_{run_number}.h5" +) as reader: + for key, trigger in trigger_type.items(): + try: + tableReader = reader.read( + table_name=f"/data/ChargesContainer_0/{trigger.name}", + containers=ChargesContainer, + ) container.containers[trigger] = next(tableReader) except Exception as err: print(err) @@ -70,27 +91,20 @@ # %% container.containers -# %% -from nectarchain.makers.component import ChargesComponent import matplotlib.pyplot as plt import numpy as np # %% -counts,charge = ChargesComponent.histo_hg(container.containers[EventType.FLATFIELD]) -charge.shape,counts.shape +from nectarchain.makers.component import ChargesComponent # %% -fig,ax = plt.subplots(1,1,figsize = (7,7)) -ax.errorbar( - charge[30], - counts[30], - np.sqrt(counts[30]), - zorder=0, - fmt=".", - label="data" - ) +counts, charge = ChargesComponent.histo_hg(container.containers[EventType.FLATFIELD]) +charge.shape, counts.shape # %% +fig, ax = plt.subplots(1, 1, figsize=(7, 7)) +ax.errorbar( + charge[30], counts[30], np.sqrt(counts[30]), zorder=0, fmt=".", label="data" +) - - +# %% diff --git a/src/nectarchain/user_scripts/ggrolleron/test.py b/src/nectarchain/user_scripts/ggrolleron/test.py index 031c2841..b36e784c 100644 --- a/src/nectarchain/user_scripts/ggrolleron/test.py +++ b/src/nectarchain/user_scripts/ggrolleron/test.py @@ -68,7 +68,7 @@ def test_extractor(): # wfs = WaveformsContainer.load(f"{os.environ['NECTARCAMDATA']}/waveforms/waveforms_run{run_number[0]}.fits") - # charge = ChargeContainer.compute_charge(spe_run_1000V,1,method = "gradient_extractor",) + # charge = ChargeContainer.reload_events(spe_run_1000V,1,method = "gradient_extractor",) t = time.time() charge = ChargeContainer.from_waveforms( @@ -121,7 +121,7 @@ def test_simplecontainer(): # spe_run_1000V.load(f"{os.environ['NECTARCAMDATA']}/waveforms/waveforms_run{run_number[0]}.fits") - charge = ChargeContainer.compute_charge( + charge = ChargeContainer.reload_events( spe_run_1000V, 1, method="LocalPeakWindowSum", diff --git a/src/nectarchain/utils/__init__.py b/src/nectarchain/utils/__init__.py index e69de29b..6f3aa3a4 100644 --- a/src/nectarchain/utils/__init__.py +++ b/src/nectarchain/utils/__init__.py @@ -0,0 +1,2 @@ +from .error import * +from .utils import * diff --git a/src/nectarchain/makers/calibration/gain/utils/error.py b/src/nectarchain/utils/error.py similarity index 100% rename from src/nectarchain/makers/calibration/gain/utils/error.py rename to src/nectarchain/utils/error.py diff --git a/src/nectarchain/makers/calibration/gain/utils/utils.py b/src/nectarchain/utils/utils.py similarity index 80% rename from src/nectarchain/makers/calibration/gain/utils/utils.py rename to src/nectarchain/utils/utils.py index 60a93dae..dd53b25b 100644 --- a/src/nectarchain/makers/calibration/gain/utils/utils.py +++ b/src/nectarchain/utils/utils.py @@ -1,3 +1,4 @@ +import importlib import logging import math @@ -11,12 +12,65 @@ log = logging.getLogger(__name__) log.handlers = logging.getLogger("__main__").handlers -# __all__ = ['UtilsMinuit','Multiprocessing'] +from ctapipe.core.component import Component -from ..parameters import Parameters + +class ComponentUtils: + @staticmethod + def is_in_non_abstract_subclasses( + component: Component, motherClass="NectarCAMComponent" + ): + from nectarchain.makers.component.core import NectarCAMComponent + + # module = importlib.import_module(f'nectarchain.makers.component.core') + is_in = False + if isinstance(component, eval(f"{motherClass}")): + is_in = True + else: + for _, value in eval(motherClass).non_abstract_subclasses().items(): + is_in = np.logical_or(is_in, component == value) + return is_in + + @staticmethod + def get_specific_traits(component: Component): + importlib.import_module(f"{component.__module__}") + traits_dict = component.class_traits() + if ComponentUtils.is_in_non_abstract_subclasses( + component, "NectarCAMComponent" + ) and not (component.SubComponents.default_value is None): + for component_name in component.SubComponents.default_value: #####CPT + _class = getattr( + importlib.import_module("nectarchain.makers.component"), + component_name, + ) + traits_dict.update(ComponentUtils.get_specific_traits(_class)) + traits_dict.pop("config", True) + traits_dict.pop("parent", True) + return traits_dict + + @staticmethod + def get_configurable_traits(component: Component): + traits_dict = ComponentUtils.get_specific_traits(component) + output_traits_dict = traits_dict.copy() + for key, item in traits_dict.items(): + if item.read_only: + output_traits_dict.pop(key) + return output_traits_dict + + @staticmethod + def get_class_name_from_ComponentName(componentName: str): + from nectarchain.makers.component.core import NectarCAMComponent + + for class_name, _class in NectarCAMComponent.non_abstract_subclasses().items(): + if componentName in class_name: + return _class + + raise ValueError( + "componentName is not a valid component, this component is not known as a child of NectarCAMComponent" + ) -class Multiprocessing: +class multiprocessing: @staticmethod def custom_error_callback(error): log.error(f"Got an error: {error}") @@ -31,7 +85,7 @@ def chi2_pvalue(ndof: int, likelihood: float): class UtilsMinuit: @staticmethod - def make_minuit_par_kwargs(parameters: Parameters): + def make_minuit_par_kwargs(parameters): """Create *Parameter Keyword Arguments* for the `Minuit` constructor. updated for Minuit >2.0