' in line : - url_data = line.split("
")[0].split('FC:')[1] - log.debug(f"url_data found {url_data}") - break - - if i == len(lines)-1 : - e=Exception("lfns not found on GRID") - log.error(e,exc_info=True) - log.debug(lines) - raise e - - if output_lfns : - lfns = [] - try : - #Dirac - dirac = Dirac() - loc = f"/vo.cta.in2p3.fr/nectarcam/{url_data.split('/')[-2]}/{url_data.split('/')[-1]}" - log.debug(f"searching in Dirac filecatalog at {loc}") - res = dirac.listCatalogDirectory(loc, printOutput=True) - - for key in res['Value']['Successful'][loc]['Files'].keys(): - if str(run_number) in key and "fits.fz" in key : - lfns.append(key) - except Exception as e : - log.error(e,exc_info = True) - return lfns - else : - return url_data - - -class ChainGenerator(): - @staticmethod - def chain(a : Generator ,b : Generator) : - """generic metghod to chain 2 generators - - Args: - a (Generator): generator to chain - b (Generator): generator to chain - - Yields: - Generator: a chain of a and b - """ - yield from a - yield from b - - - @staticmethod - def chainEventSource(list : list,max_events : int = None) : #useless with ctapipe_io_nectarcam.NectarCAMEventSource - """recursive method to chain a list of ctapipe.io.EventSource (which may be associated to the list of *.fits.fz file of one run) - - Args: - list (EventSource): a list of EventSource - - Returns: - Generator: a generator which chains EventSource - """ - if len(list) == 2 : - return ChainGenerator.chain(list[0],list[1]) - else : - return ChainGenerator.chain(list[0],ChainGenerator.chainEventSource(list[1:])) - -class CtaPipeExtractor(): - def get_image_peak_time(cameraContainer : DL1CameraContainer) : - return cameraContainer.image, cameraContainer.peak_time diff --git a/src/nectarchain/data/__init__.py b/src/nectarchain/data/__init__.py new file mode 100644 index 00000000..1506a100 --- /dev/null +++ b/src/nectarchain/data/__init__.py @@ -0,0 +1,2 @@ +from .container import * +from .management import * diff --git a/src/nectarchain/data/container/__init__.py b/src/nectarchain/data/container/__init__.py new file mode 100644 index 00000000..87fa166c --- /dev/null +++ b/src/nectarchain/data/container/__init__.py @@ -0,0 +1,3 @@ +from .chargesContainer import * +from .core import * +from .waveformsContainer import * diff --git a/src/nectarchain/data/container/chargesContainer.py b/src/nectarchain/data/container/chargesContainer.py new file mode 100644 index 00000000..42fdad88 --- /dev/null +++ b/src/nectarchain/data/container/chargesContainer.py @@ -0,0 +1,38 @@ +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 +from abc import ABC +from pathlib import Path + +import numpy as np +from astropy.io import fits +from ctapipe.containers import Field,partial,Map + +from .core import ArrayDataContainer,TriggerMapContainer + +class ChargesContainer(ArrayDataContainer): + """ + A container that holds information about charges from a specific run. + + Fields: + charges_hg (np.ndarray): An array of high gain charges. + charges_lg (np.ndarray): An array of low gain charges. + peak_hg (np.ndarray): An array of high gain peak time. + peak_lg (np.ndarray): An array of low gain peak time. + method (str): The charge extraction method used. + """ + + charges_hg = Field(type=np.ndarray, dtype = np.uint16, ndim = 2, description="The high gain charges") + charges_lg = Field(type=np.ndarray, dtype = np.uint16, ndim = 2, description="The low gain charges") + peak_hg = Field(type=np.ndarray, dtype = np.uint16, ndim = 2, description="The high gain peak time") + peak_lg = Field(type=np.ndarray, dtype = np.uint16, ndim = 2, description="The low gain peak time") + method = Field(type=str, description="The charge extraction method used") + +class ChargesContainers(TriggerMapContainer): + containers = Field(default_factory=partial(Map, ChargesContainer), + description = "trigger mapping of ChargesContainer" + ) \ No newline at end of file diff --git a/src/nectarchain/data/container/core.py b/src/nectarchain/data/container/core.py new file mode 100644 index 00000000..ca6dffd4 --- /dev/null +++ b/src/nectarchain/data/container/core.py @@ -0,0 +1,70 @@ +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 numpy as np +from ctapipe.containers import Field,Container,partial,Map + + +__all__ = ["ArrayDataContainer","TriggerMapContainer"] + +class NectarCAMContainer(Container) : + """base class for the NectarCAM containers. This contaner cannot berecursive, to be directly written with a HDF5TableWriter """ + +class ArrayDataContainer(NectarCAMContainer): + """ + A container that holds information about waveforms from a specific run. + + Attributes: + run_number (int): The run number associated with the waveforms. + nevents (int): The number of events. + npixels (int): The number of pixels. + camera (str): The name of the camera. + pixels_id (np.ndarray): An array of pixel IDs. + broken_pixels_hg (np.ndarray): An array of high gain broken pixels. + broken_pixels_lg (np.ndarray): An array of low gain broken pixels. + ucts_timestamp (np.ndarray): An array of events' UCTS timestamps. + ucts_busy_counter (np.ndarray): An array of UCTS busy counters. + ucts_event_counter (np.ndarray): An array of UCTS event counters. + event_type (np.ndarray): An array of trigger event types. + event_id (np.ndarray): An array of event IDs. + trig_pattern_all (np.ndarray): An array of trigger patterns. + trig_pattern (np.ndarray): An array of reduced trigger patterns. + multiplicity (np.ndarray): An array of events' multiplicities. + """ + + run_number = Field( + type=int, + description="run number associated to the waveforms", + ) + nevents = Field( + type=int, + description="number of events", + ) + npixels = Field( + type=int, + description="number of effective pixels", + ) + pixels_id = Field(type=np.ndarray, dtype = np.uint16, ndim = 1,description="pixel ids") + broken_pixels_hg = Field(type=np.ndarray, dtype = bool, ndim = 2, description="high gain broken pixels") + broken_pixels_lg = Field(type=np.ndarray, dtype = bool, ndim = 2, description="low gain broken pixels") + camera = Field( + type=str, + description="camera name", + ) + ucts_timestamp = Field(type=np.ndarray, dtype = np.uint64, ndim = 1, description="events ucts timestamp") + ucts_busy_counter = Field(type=np.ndarray, dtype = np.uint32, ndim = 1, description="ucts busy counter") + ucts_event_counter = Field(type=np.ndarray, dtype = np.uint32, ndim = 1, description="ucts event counter") + 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") + trig_pattern_all = Field(type=np.ndarray, dtype = bool, ndim = 3, description="trigger pattern") + trig_pattern = Field(type=np.ndarray, dtype = bool, ndim = 2, description="reduced trigger pattern") + multiplicity = Field(type=np.ndarray, dtype = np.uint16, ndim = 1, description="events multiplicity") + +class TriggerMapContainer(Container) : + containers = Field(default_factory=partial(Map, Container), + description = "trigger mapping of Container" + ) diff --git a/src/nectarchain/data/container/tests/test_charge.py b/src/nectarchain/data/container/tests/test_charge.py new file mode 100644 index 00000000..270d5056 --- /dev/null +++ b/src/nectarchain/data/container/tests/test_charge.py @@ -0,0 +1,165 @@ +import glob + +import numpy as np + +from nectarchain.data.container import ChargesContainer, ChargesContainerIO +from nectarchain.makers import ChargesMaker + + +def create_fake_chargeContainer(): + nevents = TestChargesContainer.nevents + npixels = TestChargesContainer.npixels + rng = np.random.default_rng() + return ChargesContainer( + pixels_id=np.array([2, 4, 3, 8, 6, 9, 7, 1, 5, 10]), + nevents=nevents, + npixels=npixels, + charges_hg=rng.integers(low=0, high=1000, size=(nevents, npixels)), + charges_lg=rng.integers(low=0, high=1000, size=(nevents, npixels)), + peak_hg=rng.integers(low=0, high=60, size=(nevents, npixels)), + peak_lg=rng.integers(low=0, high=60, size=(nevents, npixels)), + run_number=TestChargesContainer.run_number, + camera="TEST", + broken_pixels_hg=rng.integers(low=0, high=1, size=(nevents, npixels)), + broken_pixels_lg=rng.integers(low=0, high=1, size=(nevents, npixels)), + ucts_timestamp=rng.integers(low=0, high=100, size=(nevents)), + ucts_busy_counter=rng.integers(low=0, high=100, size=(nevents)), + ucts_event_counter=rng.integers(low=0, high=100, size=(nevents)), + event_type=rng.integers(low=0, high=1, size=(nevents)), + event_id=rng.integers(low=0, high=1000, size=(nevents)), + trig_pattern_all=rng.integers(low=0, high=1, size=(nevents, npixels, 4)), + trig_pattern=rng.integers(low=0, high=1, size=(nevents, npixels)), + multiplicity=rng.integers(low=0, high=1, size=(nevents)), + ) + + +class TestChargesContainer: + run_number = 1234 + nevents = 140 + npixels = 10 + + # Tests that a ChargeContainer object can be created with valid input parameters. + def test_create_charge_container(self): + pixels_id = np.array([2, 4, 3, 8, 6, 9, 7, 1, 5, 10]) + nevents = TestChargesContainer.nevents + npixels = TestChargesContainer.npixels + run_number = TestChargesContainer.run_number + charges_hg = np.random.randn(nevents, npixels) + charges_lg = np.random.randn(nevents, npixels) + peak_hg = np.random.randn(nevents, npixels) + peak_lg = np.random.randn(nevents, npixels) + method = "FullWaveformSum" + charge_container = ChargesContainer( + charges_hg=charges_hg, + charges_lg=charges_lg, + peak_hg=peak_hg, + peak_lg=peak_lg, + run_number=run_number, + pixels_id=pixels_id, + nevents=nevents, + npixels=npixels, + method=method, + ) + + assert np.allclose(charge_container.charges_hg, charges_hg) + assert np.allclose(charge_container.charges_lg, charges_lg) + assert np.allclose(charge_container.peak_hg, peak_hg) + assert np.allclose(charge_container.peak_lg, peak_lg) + assert charge_container.run_number == run_number + assert charge_container.pixels_id.tolist() == pixels_id.tolist() + assert charge_container.nevents == nevents + assert charge_container.npixels == npixels + assert charge_container.method == method + + # Tests that the from_waveforms method can be called with a valid waveformContainer and method parameter. + # def test_from_waveforms_valid_input(self): + # waveform_container = WaveformsContainer(...) + # method = 'FullWaveformSum' + # + # charge_container = ChargeContainer.from_waveforms(waveform_container, method) + # + # assert isinstance(charge_container, ChargeContainer) + + # Tests that the ChargeContainer object can be written to a file and the file is created. + def test_write_charge_container(self, tmp_path="/tmp"): + charge_container = create_fake_chargeContainer() + tmp_path += f"/{np.random.randn(1)[0]}" + + ChargesContainerIO.write(tmp_path, charge_container) + + assert ( + len( + glob.glob( + f"{tmp_path}/charge_run{TestChargesContainer.run_number}.fits" + ) + ) + == 1 + ) + + # Tests that a ChargeContainer object can be loaded from a file and the object is correctly initialized. + def test_load_charge_container(self, tmp_path="/tmp"): + charge_container = create_fake_chargeContainer() + tmp_path += f"/{np.random.randn(1)[0]}" + + ChargesContainerIO.write(tmp_path, charge_container) + + loaded_charge_container = ChargesContainerIO.load( + tmp_path, TestChargesContainer.run_number + ) + + assert isinstance(loaded_charge_container, ChargesContainer) + assert np.allclose( + loaded_charge_container.charges_hg, charge_container.charges_hg + ) + assert np.allclose( + loaded_charge_container.charges_lg, charge_container.charges_lg + ) + assert np.allclose(loaded_charge_container.peak_hg, charge_container.peak_hg) + assert np.allclose(loaded_charge_container.peak_lg, charge_container.peak_lg) + assert loaded_charge_container.run_number == charge_container.run_number + assert ( + loaded_charge_container.pixels_id.tolist() + == charge_container.pixels_id.tolist() + ) + assert loaded_charge_container.nevents == charge_container.nevents + assert loaded_charge_container.npixels == charge_container.npixels + assert loaded_charge_container.method == charge_container.method + + # Tests that the ChargeContainer object can be sorted by event_id and the object is sorted accordingly. + def test_sort_charge_container(self): + charge_container = create_fake_chargeContainer() + + sorted_charge_container = ChargesMaker.sort(charge_container) + + assert sorted_charge_container.event_id.tolist() == sorted( + charge_container.event_id.tolist() + ) + + # Tests that the run_number, pixels_id, npixels, nevents, method, multiplicity, and trig_pattern properties of the ChargeContainer object can be accessed and the values are correct. + def test_access_properties(self): + pixels_id = np.array([2, 4, 3, 8, 6, 9, 7, 1, 5, 10]) + nevents = 40 + npixels = 10 + charges_hg = np.random.randn(nevents, npixels) + charges_lg = np.random.randn(nevents, npixels) + peak_hg = np.random.randn(nevents, npixels) + peak_lg = np.random.randn(nevents, npixels) + run_number = 1234 + method = "FullWaveformSum" + charge_container = ChargesContainer( + charges_hg=charges_hg, + charges_lg=charges_lg, + peak_hg=peak_hg, + peak_lg=peak_lg, + run_number=run_number, + pixels_id=pixels_id, + nevents=nevents, + npixels=npixels, + method=method, + ) + + assert charge_container.run_number == run_number + assert charge_container.pixels_id.tolist() == pixels_id.tolist() + assert charge_container.npixels == npixels + assert charge_container.nevents == nevents + assert charge_container.method == method diff --git a/src/nectarchain/data/container/tests/test_waveforms.py b/src/nectarchain/data/container/tests/test_waveforms.py new file mode 100644 index 00000000..be38c45f --- /dev/null +++ b/src/nectarchain/data/container/tests/test_waveforms.py @@ -0,0 +1,97 @@ +import glob + +import numpy as np +from ctapipe.instrument import SubarrayDescription + +from nectarchain.data.container import WaveformsContainer, WaveformsContainerIO +from nectarchain.makers import WaveformsMaker + + +def create_fake_waveformsContainer(): + nevents = TestWaveformsContainer.nevents + npixels = TestWaveformsContainer.npixels + nsamples = TestWaveformsContainer.nsamples + rng = np.random.default_rng() + faked_subarray = SubarrayDescription(name="TEST") + + return WaveformsContainer( + pixels_id=np.array([2, 4, 3, 8, 6, 9, 7, 1, 5, 10]), + nevents=nevents, + npixels=npixels, + wfs_hg=rng.integers(low=0, high=1000, size=(nevents, npixels, nsamples)), + wfs_lg=rng.integers(low=0, high=1000, size=(nevents, npixels, nsamples)), + run_number=TestWaveformsContainer.run_number, + camera="TEST", + subarray=faked_subarray, + broken_pixels_hg=rng.integers(low=0, high=1, size=(nevents, npixels)), + broken_pixels_lg=rng.integers(low=0, high=1, size=(nevents, npixels)), + ucts_timestamp=rng.integers(low=0, high=100, size=(nevents)), + ucts_busy_counter=rng.integers(low=0, high=100, size=(nevents)), + ucts_event_counter=rng.integers(low=0, high=100, size=(nevents)), + event_type=rng.integers(low=0, high=1, size=(nevents)), + event_id=rng.integers(low=0, high=1000, size=(nevents)), + trig_pattern_all=rng.integers(low=0, high=1, size=(nevents, npixels, 4)), + trig_pattern=rng.integers(low=0, high=1, size=(nevents, npixels)), + multiplicity=rng.integers(low=0, high=1, size=(nevents)), + ) + + +class TestWaveformsContainer: + run_number = 1234 + nevents = 140 + npixels = 10 + nsamples = 5 + + # Tests that a ChargeContainer object can be created with valid input parameters. + def test_create_waveform_container(self): + waveform_container = create_fake_waveformsContainer() + assert isinstance(waveform_container, WaveformsContainer) + + # Tests that the ChargeContainer object can be written to a file and the file is created. + def test_write_waveform_container(self, tmp_path="/tmp"): + waveform_container = create_fake_waveformsContainer() + tmp_path += f"/{np.random.randn(1)[0]}" + + WaveformsContainerIO.write(tmp_path, waveform_container) + + assert ( + len(glob.glob(f"{tmp_path}/*_run{TestWaveformsContainer.run_number}.fits")) + == 1 + ) + + # Tests that a ChargeContainer object can be loaded from a file and the object is correctly initialized. + def test_load_waveform_container(self, tmp_path="/tmp"): + waveform_container = create_fake_waveformsContainer() + tmp_path += f"/{np.random.randn(1)[0]}" + + WaveformsContainerIO.write(tmp_path, waveform_container) + + loaded_waveform_container = WaveformsContainerIO.load( + tmp_path, TestWaveformsContainer.run_number + ) + + assert isinstance(loaded_waveform_container, WaveformsContainer) + assert np.allclose(loaded_waveform_container.wfs_hg, waveform_container.wfs_hg) + assert np.allclose(loaded_waveform_container.wfs_lg, waveform_container.wfs_lg) + assert loaded_waveform_container.run_number == waveform_container.run_number + assert ( + loaded_waveform_container.pixels_id.tolist() + == waveform_container.pixels_id.tolist() + ) + assert loaded_waveform_container.nevents == waveform_container.nevents + assert loaded_waveform_container.npixels == waveform_container.npixels + assert loaded_waveform_container.nsamples == waveform_container.nsamples + + # Tests that the ChargeContainer object can be sorted by event_id and the object is sorted accordingly. + def test_sort_waveform_container(self): + waveform_container = create_fake_waveformsContainer() + + sorted_waveform_container = WaveformsMaker.sort(waveform_container) + + assert sorted_waveform_container.event_id.tolist() == sorted( + waveform_container.event_id.tolist() + ) + + +if __name__ == "__main__": + TestWaveformsContainer().test_create_waveform_container() diff --git a/src/nectarchain/data/container/waveformsContainer.py b/src/nectarchain/data/container/waveformsContainer.py new file mode 100644 index 00000000..0eda692c --- /dev/null +++ b/src/nectarchain/data/container/waveformsContainer.py @@ -0,0 +1,48 @@ +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 os +from abc import ABC +from enum import Enum +from pathlib import Path + +import astropy.units as u +import numpy as np +from astropy.io import fits +from astropy.table import Column, QTable, Table +from ctapipe.containers import Field,Container,partial,Map +from ctapipe.instrument.subarray import SubarrayDescription +from tqdm import tqdm + +from .core import ArrayDataContainer,TriggerMapContainer + + +class WaveformsContainer(ArrayDataContainer): + """ + A container that holds information about waveforms from a specific run. + + Fields: + nsamples (int): The number of samples in the waveforms. + subarray (SubarrayDescription): The subarray description instance. + wfs_hg (np.ndarray): An array of high gain waveforms. + wfs_lg (np.ndarray): An array of low gain waveforms. + """ + + nsamples = Field( + type=int, + description="number of samples in the waveforms", + ) + #subarray = Field(type=SubarrayDescription, description="The subarray description") + wfs_hg = Field(type=np.ndarray, dtype = np.uint16, ndim = 3, description="high gain waveforms") + wfs_lg = Field(type=np.ndarray, dtype = np.uint16, ndim = 3, description="low gain waveforms") + + +class WaveformsContainers(TriggerMapContainer): + containers = Field(default_factory=partial(Map, WaveformsContainer), + description = "trigger mapping of WaveformContainer" + ) + diff --git a/src/nectarchain/data/management.py b/src/nectarchain/data/management.py new file mode 100644 index 00000000..b524bc40 --- /dev/null +++ b/src/nectarchain/data/management.py @@ -0,0 +1,164 @@ +import logging + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + +import glob +import os +from pathlib import Path +from typing import List, Tuple + +__all__ = ["DataManagement"] + + +class DataManagement: + @staticmethod + def findrun(run_number: int, search_on_GRID=True) -> Tuple[Path, List[Path]]: + """method to find in NECTARCAMDATA the list of *.fits.fz files associated to run_number + + Args: + run_number (int): the run number + + Returns: + (PosixPath,list): the path list of *fits.fz files + """ + basepath = os.environ["NECTARCAMDATA"] + list = glob.glob( + basepath + "**/*" + str(run_number) + "*.fits.fz", recursive=True + ) + list_path = [Path(chemin) for chemin in list] + if len(list_path) == 0: + e = FileNotFoundError(f"run {run_number} is not present in {basepath}") + if search_on_GRID: + log.warning(e, exc_info=True) + log.info("will search files on GRID and fetch them") + lfns = DataManagement.get_GRID_location(run_number) + DataManagement.getRunFromDIRAC(lfns) + list = glob.glob( + basepath + "**/*" + str(run_number) + "*.fits.fz", recursive=True + ) + list_path = [Path(chemin) for chemin in list] + else: + log.error(e, exc_info=True) + raise e + + name = list_path[0].name.split(".") + name[2] = "*" + name = Path(str(list_path[0].parent)) / ( + f"{name[0]}.{name[1]}.{name[2]}.{name[3]}.{name[4]}" + ) + log.info(f"Found {len(list_path)} files matching {name}") + + # to sort list path + _sorted = sorted([[file, int(file.suffixes[1][1:])] for file in list_path]) + list_path = [_sorted[i][0] for i in range(len(_sorted))] + + return name, list_path + + @staticmethod + def getRunFromDIRAC(lfns: list): + """method do get run files from GRID-EGI from input lfns + + Args: + lfns (list): list of lfns path + """ + from DIRAC.Interfaces.API.Dirac import Dirac + + dirac = Dirac() + for lfn in lfns: + if not ( + os.path.exists(f'{os.environ["NECTARCAMDATA"]}/{os.path.basename(lfn)}') + ): + dirac.getFile( + lfn=lfn, destDir=os.environ["NECTARCAMDATA"], printOutput=True + ) + + @staticmethod + def get_GRID_location( + run_number: int, output_lfns=True, username=None, password=None + ): + """method to get run location on GRID from Elog (work in progress!) + + Args: + run_number (int): run number + output_lfns (bool, optional): if True, return lfns path of fits.gz files, else return parent directory of run location. Defaults to True. + username (_type_, optional): username for Elog login. Defaults to None. + password (_type_, optional): password for Elog login. Defaults to None. + + Returns: + _type_: _description_ + """ + import browser_cookie3 + import mechanize + import requests + + url = "http://nectarcam.in2p3.fr/elog/nectarcam-data-qm/?cmd=Find" + + # url_run = f"http://nectarcam.in2p3.fr/elog/nectarcam-data-qm/?mode=full&reverse=0&reverse=1&npp=20&subtext=%23{run_number}" + + if not (username is None or password is None): + log.debug("log to Elog with username and password") + # log to Elog + br = mechanize.Browser() + br.open(url) + # form = br.select_form("form1") + for i in range(4): + log.debug(br.form.find_control(nr=i).name) + br.form["uname"] = username + br.form["upassword"] = password + br.method = "POST" + req = br.submit() + # html_page = req.get_data() + cookies = br._ua_handlers["_cookies"].cookiejar + # get data + req = requests.get( + f"http://nectarcam.in2p3.fr/elog/nectarcam-data-qm/?jcmd=&mode=Raw&attach=1&printable=1&reverse=0&reverse=1&npp=20&ma=&da=&ya=&ha=&na=&ca=&last=&mb=&db=&yb=&hb=&nb=&cb=&Author=&Setup=&Category=&Keyword=&Subject=%23{run_number}&ModuleCount=&subtext=", + cookies=cookies, + ) + + else: + # try to acces data by getting cookies from firefox and Chrome + log.debug("try to get data with cookies from Firefox abnd Chrome") + cookies = browser_cookie3.load() + req = requests.get( + f"http://nectarcam.in2p3.fr/elog/nectarcam-data-qm/?jcmd=&mode=Raw&attach=1&printable=1&reverse=0&reverse=1&npp=20&ma=&da=&ya=&ha=&na=&ca=&last=&mb=&db=&yb=&hb=&nb=&cb=&Author=&Setup=&Category=&Keyword=&Subject=%23{run_number}&ModuleCount=&subtext=", + cookies=cookies, + ) + + # if "" in line: + url_data = line.split("
")[0].split("FC:")[1] + log.debug(f"url_data found {url_data}") + break + + if i == len(lines) - 1: + e = Exception("lfns not found on GRID") + log.error(e, exc_info=True) + log.debug(lines) + raise e + + if output_lfns: + lfns = [] + try: + # Dirac + from DIRAC.Interfaces.API.Dirac import Dirac + + dirac = Dirac() + loc = f"/vo.cta.in2p3.fr/nectarcam/{url_data.split('/')[-2]}/{url_data.split('/')[-1]}" + log.debug(f"searching in Dirac filecatalog at {loc}") + res = dirac.listCatalogDirectory(loc, printOutput=True) + + for key in res["Value"]["Successful"][loc]["Files"].keys(): + if str(run_number) in key and "fits.fz" in key: + lfns.append(key) + except Exception as e: + log.error(e, exc_info=True) + return lfns + else: + return url_data diff --git a/src/nectarchain/display/__init__.py b/src/nectarchain/display/__init__.py new file mode 100644 index 00000000..4941e7d3 --- /dev/null +++ b/src/nectarchain/display/__init__.py @@ -0,0 +1 @@ +from .display import * diff --git a/src/nectarchain/display/display.py b/src/nectarchain/display/display.py new file mode 100644 index 00000000..5bb2e414 --- /dev/null +++ b/src/nectarchain/display/display.py @@ -0,0 +1,70 @@ +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 glob +import os +import sys +from abc import ABC +from argparse import ArgumentError +from enum import Enum +from pathlib import Path + +import astropy.units as u +import numpy as np +from astropy.io import fits +from astropy.table import Column, QTable, Table +from ctapipe.containers import EventType +from ctapipe.coordinates import CameraFrame, EngineeringCameraFrame +from ctapipe.instrument import CameraGeometry, SubarrayDescription, TelescopeDescription +from ctapipe.visualization import CameraDisplay +from ctapipe_io_nectarcam import NectarCAMEventSource +from matplotlib import pyplot as plt +from tqdm import tqdm + +from ..data import DataManagement +from ..data.container import ArrayDataContainer, ChargesContainer, WaveformsContainer + + +class ContainerDisplay(ABC): + @staticmethod + def display(container: ArrayDataContainer, evt, geometry, cmap="gnuplot2"): + """plot camera display for HIGH GAIN channel + + Args: + evt (int): event index + cmap (str, optional): colormap. Defaults to 'gnuplot2'. + + Returns: + CameraDisplay: thoe cameraDisplay plot + """ + if isinstance(container, ChargesContainer): + image = container.charges_hg + elif isinstance(container, WaveformsContainer): + image = container.wfs_hg.sum(axis=2) + else: + log.warning("container can't be displayed") + disp = CameraDisplay(geometry=geometry, image=image[evt], cmap=cmap) + disp.add_colorbar() + return disp + + @staticmethod + def plot_waveform(waveformsContainer: WaveformsContainer, evt, **kwargs): + """plot the waveform of the evt in the HIGH GAIN channel + + Args: + evt (int): the event index + + Returns: + tuple: the figure and axes + """ + if "figure" in kwargs.keys() and "ax" in kwargs.keys(): + fig = kwargs.get("figure") + ax = kwargs.get("ax") + else: + fig, ax = plt.subplots(1, 1) + ax.plot(waveformsContainer.wfs_hg[evt].T) + return fig, ax diff --git a/src/nectarchain/dqm/start_calib.py b/src/nectarchain/dqm/start_calib.py index 6eef50b7..a9b08d11 100644 --- a/src/nectarchain/dqm/start_calib.py +++ b/src/nectarchain/dqm/start_calib.py @@ -52,7 +52,7 @@ if args.runnb is not None: # Grab runs automatically from DIRAC is the -r option is provided - from nectarchain.calibration.container import utils + 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 new file mode 100644 index 00000000..f059c4fa --- /dev/null +++ b/src/nectarchain/makers/__init__.py @@ -0,0 +1,4 @@ +# from .chargesMakers import * +from .core import * +from .waveformsMakers import * +from .chargesMakers import * diff --git a/src/nectarchain/makers/calibration/__init__.py b/src/nectarchain/makers/calibration/__init__.py new file mode 100644 index 00000000..81f96f56 --- /dev/null +++ b/src/nectarchain/makers/calibration/__init__.py @@ -0,0 +1,3 @@ +# 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 new file mode 100644 index 00000000..58f0c991 --- /dev/null +++ b/src/nectarchain/makers/calibration/core.py @@ -0,0 +1,150 @@ +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 +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 + +__all__ = [""] + + +class CalibrationMaker(BaseMaker): + """ + Mother class for all calibration makers that can be defined to compute calibration coefficients from data. + + 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. + """ + + _reduced_name = "Calibration" + PIXELS_ID_COLUMN = "pixels_id" + NP_PIXELS = "npixels" + + 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) + self.__results = QTable() + self.__results.add_column( + Column( + self.__pixels_id, + __class__.PIXELS_ID_COLUMN, + unit=u.dimensionless_unscaled, + ) + ) + self.__results.meta[__class__.NP_PIXELS] = self.npixels + self.__results.meta[ + "comments" + ] = f'Produced with NectarChain, Credit : CTA NectarCam {date.today().strftime("%B %d, %Y")}' + + def save(self, path, **kwargs): + """ + Saves the results to a file in the specified path. + + Args: + path (str): The path to save the results. + **kwargs: Additional keyword arguments. + + Keyword Args: + overwrite (bool): Whether to overwrite an existing file. Defaults to False. + """ + path = Path(path) + path.mkdir(parents=True, exist_ok=True) + log.info(f"data saved in {path}") + self._results.write( + f"{path}/results_{self._reduced_name}.ecsv", + format="ascii.ecsv", + 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): + """ + Get the result table. + + Returns: + QTable: The result table. + """ + return self.__results + + @property + def results(self): + """ + Get a copy of the result table. + + Returns: + QTable: A copy of the result table. + """ + return copy(self.__results) diff --git a/src/nectarchain/makers/calibration/flatfieldMakers.py b/src/nectarchain/makers/calibration/flatfieldMakers.py new file mode 100644 index 00000000..f7c1c5c8 --- /dev/null +++ b/src/nectarchain/makers/calibration/flatfieldMakers.py @@ -0,0 +1,19 @@ +import logging + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + +from .core import CalibrationMaker + +__all__ = ["FlatfieldMaker"] + + +class FlatfieldMaker(CalibrationMaker): + def __init__(self, *args, **kwargs) -> None: + super().__init__(*args, **kwargs) + + def make(self): + raise NotImplementedError( + "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..31834773 --- /dev/null +++ b/src/nectarchain/makers/calibration/gain/FlatFieldSPEMakers.py @@ -0,0 +1,1112 @@ +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 time +from inspect import signature +from multiprocessing import Pool +from typing import Tuple + +import astropy.units as u +import matplotlib.pyplot as plt +import numpy as np +import yaml +from astropy.table import Column, QTable +from iminuit import Minuit +from matplotlib.colors import to_rgba +from matplotlib.patches import Rectangle +from scipy.optimize import curve_fit +from scipy.signal import find_peaks, savgol_filter +from scipy.special import gammainc + +from ....data.container import ChargesContainer +from ...chargesMakers import ChargesMaker +from .gainMakers import GainMaker +from .parameters import Parameter, Parameters +from .utils import ( + MPE2, + MeanValueError, + PedestalValueError, + Statistics, + UtilsMinuit, + weight_gaussian, +) + +__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_lenght`: 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_lenght = 40 + _Order = 2 + + # constructors + def __init__(self, *args, **kwargs) -> None: + """ + Initializes the FlatFieldSPEMaker class. + + Parameters: + *args: Variable length argument list. + **kwargs: Arbitrary keyword arguments. + + Returns: + None + """ + super().__init__(*args, **kwargs) + self.__parameters = Parameters() + + @property + def npixels(self): + """ + Returns the number of pixels. + + Returns: + int: The number of pixels. + """ + return len(self._pixels_id) + + @property + def parameters(self): + """ + Returns a deep copy of the internal parameters. + + Returns: + dict: A deep copy of the internal parameters. + """ + return copy.deepcopy(self.__parameters) + + @property + def _parameters(self): + """ + Returns the internal parameters. + + Returns: + dict: The internal parameters. + """ + return self.__parameters + + # 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. + + 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. + + Returns: + None + """ + with open( + f"{os.path.dirname(os.path.abspath(__file__))}/{parameters_file}" + ) as parameters: + param = yaml.safe_load(parameters) + if only_update: + for i, name in enumerate(self.__parameters.parnames): + dico = param.get(name, False) + if dico: + self._parameters.parameters[i].value = dico.get("value") + self._parameters.parameters[i].min = dico.get("min", np.nan) + self._parameters.parameters[i].max = dico.get("max", np.nan) + else: + for name, dico in param.items(): + setattr( + self, + f"__{name}", + Parameter( + name=name, + value=dico["value"], + min=dico.get("min", np.nan), + max=dico.get("max", np.nan), + unit=dico.get("unit", u.dimensionless_unscaled), + ), + ) + self._parameters.append(eval(f"self.__{name}")) + + @staticmethod + 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. + + Args: + 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. + """ + try: + coeff_ped, coeff_mean = __class__._get_mean_gaussian_fit( + charge, counts, **kwargs + ) + pedestal = parameters["pedestal"] + pedestal.value = coeff_ped[1] + pedestal.min = coeff_ped[1] - coeff_ped[2] + pedestal.max = coeff_ped[1] + coeff_ped[2] + log.debug(f"pedestal updated: {pedestal}") + pedestalWidth = parameters["pedestalWidth"] + pedestalWidth.value = pedestal.max - pedestal.value + pedestalWidth.max = 3 * pedestalWidth.value + log.debug(f"pedestalWidth updated: {pedestalWidth.value}") + + if (coeff_mean[1] - pedestal.value < 0) or ( + (coeff_mean[1] - coeff_mean[2]) - pedestal.max < 0 + ): + raise MeanValueError("mean gaussian fit not good") + mean = parameters["mean"] + mean.value = coeff_mean[1] - pedestal.value + mean.min = (coeff_mean[1] - coeff_mean[2]) - pedestal.max + mean.max = (coeff_mean[1] + coeff_mean[2]) - pedestal.min + log.debug(f"mean updated: {mean}") + except MeanValueError as e: + log.warning(e, exc_info=True) + log.warning("mean parameters limits and starting value not changed") + except Exception as e: + log.warning(e, exc_info=True) + log.warning( + "pedestal and mean parameters limits and starting value not changed" + ) + return parameters + + @staticmethod + def _get_mean_gaussian_fit( + charge: np.ndarray, counts: np.ndarray, extension: str = "", **kwargs + ) -> Tuple[np.ndarray, np.ndarray]: + """ + Perform a Gaussian fit on the data to determine the pedestal and mean values. + + 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 "". + **kwargs: Additional keyword arguments. + + Returns: + Tuple[np.ndarray, np.ndarray]: A tuple of fit coefficients for the pedestal and mean. + + Example Usage: + flat_field_maker = FlatFieldSPEMaker() + charge = np.array([1, 2, 3, 4, 5]) + counts = np.array([10, 20, 30, 40, 50]) + coeff, coeff_mean = flat_field_maker._get_mean_gaussian_fit(charge, counts) + print(coeff) # Output: [norm,peak_value, peak_width] + print(coeff_mean) # Output: [norm,peak_value_mean, peak_width_mean] + """ + windows_lenght = __class__._Windows_lenght + order = __class__._Order + 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]] + coeff, _ = curve_fit( + weight_gaussian, + charge[: peaks[0][peak_max]], + histo_smoothed[: peaks[0][peak_max]], + p0=[peak_value, peak_pos, 1], + ) + if log.getEffectiveLevel() == logging.DEBUG and kwargs.get("display", False): + log.debug("plotting figures with prefit parameters computation") + fig, ax = plt.subplots(1, 1, figsize=(5, 5)) + ax.errorbar( + charge, counts, np.sqrt(counts), zorder=0, fmt=".", label="data" + ) + ax.plot( + charge, + histo_smoothed, + label=f"smoothed data with savgol filter (windows lenght : {windows_lenght}, order : {order})", + ) + ax.plot( + charge, + weight_gaussian(charge, coeff[0], coeff[1], coeff[2]), + label="gaussian fit of the pedestal, left tail only", + ) + ax.set_xlim([peak_pos - 500, None]) + ax.vlines( + coeff[1], + 0, + peak_value, + label=f"pedestal initial value = {coeff[1]:.0f}", + color="red", + ) + ax.add_patch( + Rectangle( + (coeff[1] - coeff[2], 0), + 2 * coeff[2], + peak_value, + fc=to_rgba("red", 0.5), + ) + ) + ax.set_xlabel("Charge (ADC)", size=15) + ax.set_ylabel("Events", size=15) + ax.legend(fontsize=7) + os.makedirs( + f"{os.environ.get('NECTARCHAIN_LOG')}/{os.getpid()}/figures/", + exist_ok=True, + ) + fig.savefig( + f"{os.environ.get('NECTARCHAIN_LOG')}/{os.getpid()}/figures/initialization_pedestal_pixel{extension}_{os.getpid()}.pdf" + ) + fig.clf() + plt.close(fig) + del fig, ax + mask = charge > coeff[1] + 3 * coeff[2] + peaks_mean = find_peaks(histo_smoothed[mask]) + peak_max_mean = np.argmax(histo_smoothed[mask][peaks_mean[0]]) + peak_pos_mean, peak_value_mean = ( + charge[mask][peaks_mean[0][peak_max_mean]], + histo_smoothed[mask][peaks_mean[0][peak_max_mean]], + ) + mask = (charge > ((coeff[1] + peak_pos_mean) / 2)) * ( + charge < (peak_pos_mean + (peak_pos_mean - coeff[1]) / 2) + ) + coeff_mean, _ = curve_fit( + weight_gaussian, + charge[mask], + histo_smoothed[mask], + p0=[peak_value_mean, peak_pos_mean, 1], + ) + if log.getEffectiveLevel() == logging.DEBUG and kwargs.get("display", False): + log.debug("plotting figures with prefit parameters computation") + fig, ax = plt.subplots(1, 1, figsize=(5, 5)) + ax.errorbar( + charge, counts, np.sqrt(counts), zorder=0, fmt=".", label="data" + ) + ax.plot( + charge, + histo_smoothed, + label=f"smoothed data with savgol filter (windows lenght : {windows_lenght}, order : {order})", + ) + ax.plot( + charge, + weight_gaussian(charge, coeff_mean[0], coeff_mean[1], coeff_mean[2]), + label="gaussian fit of the SPE", + ) + ax.vlines( + coeff_mean[1], + 0, + peak_value, + label=f"mean initial value = {coeff_mean[1] - coeff[1]:.0f}", + color="red", + ) + ax.add_patch( + Rectangle( + (coeff_mean[1] - coeff_mean[2], 0), + 2 * coeff_mean[2], + peak_value_mean, + fc=to_rgba("red", 0.5), + ) + ) + ax.set_xlim([peak_pos - 500, None]) + ax.set_xlabel("Charge (ADC)", size=15) + ax.set_ylabel("Events", size=15) + ax.legend(fontsize=7) + os.makedirs( + f"{os.environ.get('NECTARCHAIN_LOG')}/{os.getpid()}/figures/", + exist_ok=True, + ) + fig.savefig( + f"{os.environ.get('NECTARCHAIN_LOG')}/{os.getpid()}/figures/initialization_mean_pixel{extension}_{os.getpid()}.pdf" + ) + fig.clf() + 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. + """ + + for param in self._parameters.parameters: + if not (param.name in self._results.colnames): + self._results.add_column( + Column( + data=np.empty((self.npixels), dtype=np.float64), + name=param.name, + unit=param.unit, + ) + ) + self._results.add_column( + Column( + data=np.empty((self.npixels), dtype=np.float64), + name=f"{param.name}_error", + unit=param.unit, + ) + ) + + +class FlatFieldSingleHHVSPEMaker(FlatFieldSPEMaker): + """ + This class represents a FlatFieldSingleHHVSPEMaker object. + + 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: + """ + Initializes the FlatFieldSingleHHVSPEMaker object. + 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. + """ + super().__init__(*args, **kwargs) + if isinstance(charge, np.ma.masked_array): + self.__charge = charge + else: + self.__charge = np.ma.asarray(charge) + if isinstance(counts, np.ma.masked_array): + 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, + ) + ) + + @classmethod + def create_from_chargesContainer(cls, signal: ChargesContainer, **kwargs): + """ + 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) + 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__" + ) + + # getters and setters + @property + def charge(self): + """ + Returns a deep copy of the __charge attribute. + """ + return copy.deepcopy(self.__charge) + + @property + def _charge(self): + """ + Returns the __charge attribute. + """ + return self.__charge + + @property + def counts(self): + """ + Returns a deep copy of the __counts attribute. + """ + return copy.deepcopy(self.__counts) + + @property + def _counts(self): + """ + Returns the __counts attribute. + """ + return self.__counts + + # 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. + + Args: + dico (dict): A dictionary containing fit values and errors for each pixel. + pixels_id (np.ndarray): An array of pixel IDs. + + Returns: + None + """ + 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]) + 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" + ) + 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] + if key == "mean": + self._high_gain[index] = values[j] + self._results[f"high_gain_error"][index] = [ + errors[j], + errors[j], + ] + self._results[f"high_gain"][index] = values[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( + ndof, __class__.__fit_array[i].fcn(__class__.__fit_array[i].values) + ) + + @staticmethod + def _NG_Likelihood_Chi2( + pp: float, + res: float, + mu2: float, + n: float, + muped: float, + sigped: float, + lum: float, + charge: np.ndarray, + counts: np.ndarray, + **kwargs, + ): + """ + Calculates the chi-square value using the MPE2 function. + Parameters: + pp (float): The pp parameter. + res (float): The res parameter. + mu2 (float): The mu2 parameter. + n (float): The n parameter. + muped (float): The muped parameter. + sigped (float): The sigped parameter. + lum (float): The lum parameter. + charge (np.ndarray): An array of charge values. + counts (np.ndarray): An array of count values. + **kwargs: Additional keyword arguments. + Returns: + float: The chi-square value. + """ + pdf = MPE2(charge, pp, res, mu2, n, muped, sigped, lum, **kwargs) + # log.debug(f"pdf : {np.sum(pdf)}") + Ntot = np.sum(counts) + # log.debug(f'Ntot : {Ntot}') + mask = counts > 0 + Lik = np.sum( + ((pdf * Ntot - counts)[mask]) ** 2 / counts[mask] + ) # 2 times faster + return Lik + + @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. + Parameters: + charge (np.ndarray): An array of charge values. + counts (np.ndarray): An array of count values. + Returns: + function: The Chi2 function. + """ + + def Chi2( + pedestal: float, + pp: float, + luminosity: float, + resolution: float, + mean: float, + 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)) + for i in range(1000): + if gammainc(i + 1, luminosity) < 1e-5: + ntotalPE = i + break + kwargs = {"ntotalPE": ntotalPE} + return __class__._NG_Likelihood_Chi2( + pp, + resolution, + mean, + n, + pedestal, + pedestalWidth, + luminosity, + charge, + counts, + **kwargs, + ) + + return Chi2 + + # @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. + + Args: + 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. + """ + if pixels_id is None: + npix = self.npixels + pixels_id = self.pixels_id + else: + npix = len(pixels_id) + + fit_array = np.empty((npix), dtype=np.object_) + + for i, _id in enumerate(pixels_id): + index = np.where(self.pixels_id == _id)[0][0] + parameters = __class__._update_parameters( + self.parameters, + self._charge[index].data[~self._charge[index].mask], + self._counts[index].data[~self._charge[index].mask], + pixel_id=_id, + **kwargs, + ) + minuitParameters = UtilsMinuit.make_minuit_par_kwargs(parameters) + minuit_kwargs = { + parname: minuitParameters["values"][parname] + for parname in minuitParameters["values"] + } + log.info(f"creation of fit instance for pixel: {_id}") + fit_array[i] = Minuit( + __class__.cost( + self._charge[index].data[~self._charge[index].mask], + self._counts[index].data[~self._charge[index].mask], + ), + **minuit_kwargs, + ) + log.debug("fit created") + fit_array[i].errordef = Minuit.LIKELIHOOD + fit_array[i].strategy = 0 + fit_array[i].tol = 1e40 + 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) + + return fit_array + + @staticmethod + def run_fit(i: int) -> dict: + """ + Perform a fit on a specific pixel using the Minuit package. + + Args: + 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. + """ + log.info("Starting") + __class__.__fit_array[i].migrad() + __class__.__fit_array[i].hesse() + _values = np.array([params.value for params in __class__.__fit_array[i].params]) + _errors = np.array([params.error for params in __class__.__fit_array[i].params]) + log.info("Finished") + return {f"values_{i}": _values, f"errors_{i}": _errors} + + def make( + 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") + if pixels_id is None: + pixels_id = self.pixels_id + npix = self.npixels + else: + 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]}") + pixels_id = pixels_id[mask] + npix = len(pixels_id) + + if npix == 0: + log.warning("The asked pixels id are all out of the data") + return None + else: + 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) + chunksize = kwargs.get( + "chunksize", + max(__class__.__chunksize_default, npix // (nproc * 10)), + ) + log.info(f"pooling with nproc {nproc}, chunksize {chunksize}") + + t = time.time() + with Pool(nproc) as pool: + result = pool.starmap_async( + __class__.run_fit, + [(i,) for i in range(npix)], + chunksize=chunksize, + ) + result.wait() + try: + res = result.get() + except Exception as e: + log.error(e, exc_info=True) + raise e + log.debug(res) + log.info( + f"time for multiproc with starmap_async execution is {time.time() - t:.2e} sec" + ) + else: + 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") + + 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.display(pixels_id, **kwargs) + + return output + + def plot_single( + 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: + """ + Generate a plot of the data and a model fit for a specific pixel. + + 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. + pp (float): The value of the `pp` parameter. + resolution (float): The value of the `resolution` parameter. + gain (float): The value of the `gain` parameter. + gain_error (float): The value of the `gain_error` parameter. + n (float): The value of the `n` parameter. + pedestal (float): The value of the `pedestal` parameter. + pedestalWidth (float): The value of the `pedestalWidth` parameter. + luminosity (float): The value of the `luminosity` parameter. + likelihood (float): The value of the `likelihood` parameter. + + Returns: + 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") + ax.plot( + charge, + np.trapz(counts, charge) + * MPE2( + charge, + pp, + resolution, + gain, + n, + pedestal, + pedestalWidth, + luminosity, + ), + zorder=1, + linewidth=2, + 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.legend(fontsize=18) + return fig, ax + + def display(self, pixels_id: np.ndarray, **kwargs) -> None: + """ + Display and save the plot for each specified pixel ID. + + Args: + pixels_id (np.ndarray): An array of pixel IDs. + **kwargs: Additional keyword arguments. + 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()}") + 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 + + +class FlatFieldSingleHHVStdSPEMaker(FlatFieldSingleHHVSPEMaker): + """class to perform fit of the SPE signal with n and pp fixed""" + + __parameters_file = "parameters_signalStd.yaml" + _reduced_name = "FlatFieldSingleStdSPE" + + def __init__(self, charge: np.ndarray, counts: np.ndarray, *args, **kwargs) -> None: + """ + Initializes a new instance of the FlatFieldSingleHHVStdSPEMaker class. + + Args: + charge (np.ndarray): The charge data. + counts (np.ndarray): The counts data. + *args: Additional positional arguments. + **kwargs: Additional keyword arguments. + """ + super().__init__(charge, counts, *args, **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. + """ + 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" + + def __init__( + self, + charge: np.ndarray, + counts: np.ndarray, + nectarGainSPEresult: str, + same_luminosity: bool = False, + *args, + **kwargs, + ) -> None: + """ + Initializes an instance of FlatFieldSingleNominalSPEMaker. + + 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): + """ + 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. + + 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)] + + def __fix_parameters(self, same_luminosity: bool) -> 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") + 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") + 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. + + Args: + pixels_id (array-like, optional): The pixels to generate the fit array for. Defaults to None. + **kwargs: Arbitrary keyword arguments. + + Returns: + array-like: The fit array. + """ + return super()._make_fit_array_from_parameters( + pixels_id=pixels_id, + nectarGainSPEresult=self.__nectarGainSPEresult, + **kwargs, + ) + + @staticmethod + def _update_parameters( + parameters: Parameters, + charge: np.ndarray, + counts: np.ndarray, + pixel_id, + nectarGainSPEresult: QTable, + **kwargs, + ): + """ + Updates the parameters with the fixed values from the fitted data obtained from a 1400V run. + + Args: + parameters (Parameters): The parameters to update. + charge (np.ndarray): The charge values. + counts (np.ndarray): The counts values. + pixel_id (int): The pixel ID. + nectarGainSPEresult (QTable): The fitted data obtained from a 1400V run. + **kwargs: Arbitrary keyword arguments. + + Returns: + dict: The updated parameters. + """ + param = super()._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] + + resolution.value = nectarGainSPEresult[index]["resolution"].value + pp.value = nectarGainSPEresult[index]["pp"].value + n.value = nectarGainSPEresult[index]["n"].value + + if luminosity.frozen: + luminosity.value = nectarGainSPEresult[index]["luminosity"].value + return param diff --git a/src/nectarchain/makers/calibration/gain/PhotoStatisticMakers.py b/src/nectarchain/makers/calibration/gain/PhotoStatisticMakers.py new file mode 100644 index 00000000..3aa17e54 --- /dev/null +++ b/src/nectarchain/makers/calibration/gain/PhotoStatisticMakers.py @@ -0,0 +1,578 @@ +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 +from datetime import date +from pathlib import Path + +import astropy.units as u +import numpy as np +from astropy.table import Column, 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 ...chargesMakers import ChargesMaker +from .gainMakers import GainMaker + +__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( + f"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 + 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))) diff --git a/src/nectarchain/makers/calibration/gain/WhiteTargetSPEMakers.py b/src/nectarchain/makers/calibration/gain/WhiteTargetSPEMakers.py new file mode 100644 index 00000000..be726219 --- /dev/null +++ b/src/nectarchain/makers/calibration/gain/WhiteTargetSPEMakers.py @@ -0,0 +1,19 @@ +import logging + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + +from .gainMakers import GainMaker + +__all__ = ["FlatfieldMaker"] + + +class WhiteTargetSPEMaker(GainMaker): + def __init__(self, *args, **kwargs) -> None: + super().__init__(*args, **kwargs) + + def make(self): + raise NotImplementedError( + "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 new file mode 100644 index 00000000..26159534 --- /dev/null +++ b/src/nectarchain/makers/calibration/gain/__init__.py @@ -0,0 +1,4 @@ +from .FlatFieldSPEMakers import * + +# from .WhiteTargetSPEMakers import * +from .PhotoStatisticMakers import * diff --git a/src/nectarchain/makers/calibration/gain/gainMakers.py b/src/nectarchain/makers/calibration/gain/gainMakers.py new file mode 100644 index 00000000..ea8c31af --- /dev/null +++ b/src/nectarchain/makers/calibration/gain/gainMakers.py @@ -0,0 +1,135 @@ +import logging + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + +from copy import copy + +import astropy.units as u +import numpy as np +from astropy.table import Column + +from ..core import CalibrationMaker + +__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/parameters.py b/src/nectarchain/makers/calibration/gain/parameters.py new file mode 100644 index 00000000..28d384a3 --- /dev/null +++ b/src/nectarchain/makers/calibration/gain/parameters.py @@ -0,0 +1,145 @@ +import logging + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) + +import copy + +import astropy.units as u +import numpy as np + +__all__ = ["Parameter", "Parameters"] + + +class Parameter: + def __init__( + self, + name: str, + value, + min=np.nan, + max=np.nan, + error=np.nan, + unit=u.dimensionless_unscaled, + frozen: bool = False, + ) -> None: + self.__name = name + self.__value = value + self.__error = error + self.__min = min + self.__max = max + self.__unit = unit + self.__frozen = frozen + + @classmethod + def from_instance(cls, parameter): + return cls( + parameter.name, + parameter.value, + parameter.min, + parameter.max, + parameter.unit, + parameter.frozen, + ) + + def __str__(self): + return f"name : {self.__name}, value : {self.__value}, error : {self.__error}, unit : {self.__unit}, min : {self.__min}, max : {self.__max},frozen : {self.__frozen}" + + @property + def name(self): + return self.__name + + @name.setter + def name(self, value: str): + self.__name = value + + @property + def value(self): + return self.__value + + @value.setter + def value(self, value): + self.__value = value + + @property + def min(self): + return self.__min + + @min.setter + def min(self, value): + self.__min = value + + @property + def max(self): + return self.__max + + @max.setter + def max(self, value): + self.__max = value + + @property + def unit(self): + return self.__unit + + @unit.setter + def unit(self, value: u.Unit): + self.__unit = value + + @property + def error(self): + return self.__error + + @error.setter + def error(self, value): + self.__error = value + + @property + def frozen(self): + return self.__frozen + + @frozen.setter + def frozen(self, value: bool): + self.__frozen = value + + +class Parameters: + def __init__(self, parameters_liste: list = []) -> None: + self.__parameters = copy.deepcopy(parameters_liste) + + def append(self, parameter: Parameter) -> None: + self.__parameters.append(parameter) + + def __getitem__(self, key: str): + for parameter in self.__parameters: + if parameter.name == key: + return parameter + return [] + + def __str__(self): + string = "" + for parameter in self.__parameters: + string += str(parameter) + "\n" + return string + + @property + def parameters(self): + return self.__parameters + + @property + def size(self): + return len(self.__parameters) + + @property + def parnames(self): + return [parameter.name for parameter in self.__parameters] + + @property + def parvalues(self): + return [parameter.value for parameter in self.__parameters] + + @property + def unfrozen(self): + parameters = Parameters() + for parameter in self.__parameters: + if not (parameter.frozen): + parameters.append(parameter) + return parameters diff --git a/src/nectarchain/calibration/NectarGain/SPEfit/parameters_signal.yaml b/src/nectarchain/makers/calibration/gain/parameters_signal.yaml similarity index 90% rename from src/nectarchain/calibration/NectarGain/SPEfit/parameters_signal.yaml rename to src/nectarchain/makers/calibration/gain/parameters_signal.yaml index df2ec0f0..89ece153 100644 --- a/src/nectarchain/calibration/NectarGain/SPEfit/parameters_signal.yaml +++ b/src/nectarchain/makers/calibration/gain/parameters_signal.yaml @@ -10,7 +10,7 @@ max : 5.0 }, pp : { - value : 0.3735, + value : 0.45, min : 0.2, max : 0.8 }, @@ -20,7 +20,7 @@ max : 0.7 }, n : { - value : 0.708, + value : 0.697, min : 0.5, max : 0.9 }, diff --git a/src/nectarchain/calibration/NectarGain/SPEfit/parameters_signalStd.yaml b/src/nectarchain/makers/calibration/gain/parameters_signalStd.yaml similarity index 90% rename from src/nectarchain/calibration/NectarGain/SPEfit/parameters_signalStd.yaml rename to src/nectarchain/makers/calibration/gain/parameters_signalStd.yaml index b887be64..12746d7f 100644 --- a/src/nectarchain/calibration/NectarGain/SPEfit/parameters_signalStd.yaml +++ b/src/nectarchain/makers/calibration/gain/parameters_signalStd.yaml @@ -10,7 +10,7 @@ max : 5.0 }, pp : { - value : 0.463, + value : 0.45, min : .NAN, max : .NAN }, @@ -20,7 +20,7 @@ max : 0.7 }, n : { - value : 0.698, + value : 0.697, min : .NAN, max : .NAN }, diff --git a/src/nectarchain/calibration/NectarGain/SPEfit/parameters_signal_combined.yaml b/src/nectarchain/makers/calibration/gain/parameters_signal_combined.yaml similarity index 92% rename from src/nectarchain/calibration/NectarGain/SPEfit/parameters_signal_combined.yaml rename to src/nectarchain/makers/calibration/gain/parameters_signal_combined.yaml index 93dd3124..a04fc05a 100644 --- a/src/nectarchain/calibration/NectarGain/SPEfit/parameters_signal_combined.yaml +++ b/src/nectarchain/makers/calibration/gain/parameters_signal_combined.yaml @@ -15,7 +15,7 @@ max : 5.0 }, pp : { - value : 0.463, + value : 0.45, min : .NAN, max : .NAN }, @@ -25,7 +25,7 @@ max : 0.7 }, n : { - value : 0.698, + value : 0.697, min : .NAN, max : .NAN }, diff --git a/src/nectarchain/calibration/NectarGain/SPEfit/parameters_signal_fromHHVFit.yaml b/src/nectarchain/makers/calibration/gain/parameters_signal_fromHHVFit.yaml similarity index 90% rename from src/nectarchain/calibration/NectarGain/SPEfit/parameters_signal_fromHHVFit.yaml rename to src/nectarchain/makers/calibration/gain/parameters_signal_fromHHVFit.yaml index ea6a7a0e..62515304 100644 --- a/src/nectarchain/calibration/NectarGain/SPEfit/parameters_signal_fromHHVFit.yaml +++ b/src/nectarchain/makers/calibration/gain/parameters_signal_fromHHVFit.yaml @@ -10,7 +10,7 @@ max : 5.0 }, pp : { - value : 0.463, + value : 0.45, min : .NAN, max : .NAN }, @@ -20,7 +20,7 @@ max : .NAN }, n : { - value : 0.698, + value : 0.697, min : .NAN, max : .NAN }, diff --git a/src/nectarchain/makers/calibration/gain/tests/test_FlatFieldSPEMakers.py b/src/nectarchain/makers/calibration/gain/tests/test_FlatFieldSPEMakers.py new file mode 100644 index 00000000..2bd60e2b --- /dev/null +++ b/src/nectarchain/makers/calibration/gain/tests/test_FlatFieldSPEMakers.py @@ -0,0 +1,160 @@ +import astropy.units as u +import numpy as np +import pytest + +from nectarchain.data.container import ChargesContainer +from nectarchain.makers.calibration.gain import ( + FlatFieldSingleHHVSPEMaker, + FlatFieldSingleHHVStdSPEMaker, +) +from nectarchain.makers.calibration.gain.FlatFieldSPEMakers import FlatFieldSPEMaker +from nectarchain.makers.calibration.gain.parameters import Parameter, Parameters + + +class FlatFieldSPEMakerforTest(FlatFieldSPEMaker): + def make(): + pass + + +def create_fake_chargeContainer(): + pixels_id = np.array([2, 4, 3, 8, 6, 9, 7, 1, 5, 10]) + nevents = 40 + npixels = 10 + rng = np.random.default_rng() + charges_hg = rng.integers(low=0, high=1000, size=(nevents, npixels)) + charges_lg = rng.integers(low=0, high=1000, size=(nevents, npixels)) + peak_hg = rng.integers(low=0, high=60, size=(nevents, npixels)) + peak_lg = rng.integers(low=0, high=60, size=(nevents, npixels)) + run_number = 1234 + return ChargesContainer( + charges_hg=charges_hg, + charges_lg=charges_lg, + peak_hg=peak_hg, + peak_lg=peak_lg, + run_number=run_number, + pixels_id=pixels_id, + nevents=nevents, + npixels=npixels, + ) + + +class TestFlatFieldSPEMaker: + # Tests that the object can be initialized without errors + def test_initialize_object(self): + pixels_id = [2, 3, 5] + flat_field_spe_maker = FlatFieldSPEMakerforTest(pixels_id) + assert isinstance(flat_field_spe_maker, FlatFieldSPEMakerforTest) + + # Tests that parameters can be read from a YAML file + def test_read_parameters_from_yaml(self): + pixels_id = [2, 3, 5] + flat_field_spe_maker = FlatFieldSPEMakerforTest(pixels_id) + flat_field_spe_maker.read_param_from_yaml("parameters_signal.yaml") + assert flat_field_spe_maker.parameters.size == 6 + assert isinstance(flat_field_spe_maker.parameters, Parameters) + + # Tests that parameters can be updated from a YAML file + def test_update_parameters_from_yaml(self): + pixels_id = [2, 3, 5] + flat_field_spe_maker = FlatFieldSPEMakerforTest(pixels_id) + flat_field_spe_maker.read_param_from_yaml("parameters_signal.yaml") + flat_field_spe_maker.read_param_from_yaml( + "parameters_signalStd.yaml", only_update=True + ) + assert flat_field_spe_maker.parameters.parameters[-2].value == 0.697 + + # Tests that parameters can be updated from a fit + def test_update_parameters_from_fit(self): + pixels_id = [2] + flat_field_spe_maker = FlatFieldSPEMakerforTest(pixels_id) + flat_field_spe_maker.read_param_from_yaml("parameters_signal.yaml") + updated_parameters = flat_field_spe_maker._update_parameters( + flat_field_spe_maker.parameters, + charge=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10], + counts=[1, 3, 9, 5, 3, 5, 6, 3, 2, 1], + ) + + # Tests that the table can be updated from parameters + def test_update_table_from_parameters(self): + pixels_id = [2, 3, 5] + flat_field_spe_maker = FlatFieldSPEMakerforTest(pixels_id) + flat_field_spe_maker._parameters.append( + Parameter(name="param1", value=1, unit=u.dimensionless_unscaled) + ) + flat_field_spe_maker._parameters.append( + Parameter(name="param2", value=2, unit=u.dimensionless_unscaled) + ) + + flat_field_spe_maker._update_table_from_parameters() + + assert "param1" in flat_field_spe_maker._results.colnames + assert "param1_error" in flat_field_spe_maker._results.colnames + assert "param2" in flat_field_spe_maker._results.colnames + assert "param2_error" in flat_field_spe_maker._results.colnames + + +class TestFlatFieldSingleHHVSPEMaker: + # Tests that creating an instance of FlatFieldSingleHHVSPEMaker with valid input parameters is successful + def test_create_instance_valid_input(self): + charge = [1, 2, 3] + counts = [10, 20, 30] + pixels_id = [2, 3, 5] + maker = FlatFieldSingleHHVSPEMaker(charge, counts, pixels_id) + assert isinstance(maker, FlatFieldSingleHHVSPEMaker) + + # Tests that creating an instance of FlatFieldSingleHHVSPEMaker with invalid input parameters raises an error + def test_create_instance_invalid_input(self): + charge = [1, 2, 3] + counts = [10, 20] # Invalid input, counts and charge must have the same length + pixels_id = [2, 3, 5] + + with pytest.raises(Exception): + FlatFieldSingleHHVSPEMaker(charge, counts, pixels_id) + + # Tests that calling create_from_chargeContainer method with valid input parameters is successful + def test_create_from_ChargeContainer_valid_input(self): + chargeContainer = create_fake_chargeContainer() + maker = FlatFieldSingleHHVSPEMaker.create_from_chargesContainer(chargeContainer) + assert isinstance(maker, FlatFieldSingleHHVSPEMaker) + + def test_fill_results_table_from_dict(self): + pass + + def test_NG_Likelihood_Chi2(self): + pass + + def test_cost(self): + pass + + def test_make_fit_array_from_parameters(self): + pass + + def test_run_fit(self): + pass + + def test_make(self): + pass + + def test_plot_single(self): + pass + + def test_display(self): + pass + + +class TestFlatFieldSingleHHVStdSPEMaker: + def test_create_instance(self): + charge = [1, 2, 3] + counts = [ + 10, + 20, + 30, + ] # Invalid input, counts and charge must have the same length + pixels_id = [2, 3, 5] + instance = FlatFieldSingleHHVStdSPEMaker(charge, counts, pixels_id) + assert isinstance(instance, FlatFieldSingleHHVStdSPEMaker) + + +class TestFlatFieldSingleNominalSPEMaker: + def test_create_instance(self): + pass diff --git a/src/nectarchain/makers/calibration/gain/tests/test_gainMakers.py b/src/nectarchain/makers/calibration/gain/tests/test_gainMakers.py new file mode 100644 index 00000000..16be40b0 --- /dev/null +++ b/src/nectarchain/makers/calibration/gain/tests/test_gainMakers.py @@ -0,0 +1,60 @@ +# Generated by CodiumAI +from pathlib import Path + +import numpy as np +from astropy.table import QTable + +from nectarchain.makers.calibration.gain.gainMakers import GainMaker + + +class GainMakerforTest(GainMaker): + _reduced_name = "test" + + def make(): + pass + + +class TestGainMaker: + # Tests that an instance of GainMaker can be created with a list of pixel ids as input. + def test_create_instance_with_pixel_ids(self): + pixel_ids = [1, 2, 3, 4, 5] + gain_maker = GainMakerforTest(pixel_ids) + assert isinstance(gain_maker, GainMakerforTest) + assert np.array_equal(gain_maker.pixels_id, np.array(pixel_ids)) + + # Tests that high gain values can be set and retrieved for all pixels. + def test_set_and_get_high_gain_values(self): + pixel_ids = [1, 2, 3, 4, 5] + gain_maker = GainMakerforTest(pixel_ids) + high_gain_values = np.array([0.5, 0.6, 0.7, 0.8, 0.9]) + gain_maker._high_gain = high_gain_values + assert np.array_equal(gain_maker.high_gain, high_gain_values) + + # Tests that low gain values can be set and retrieved for all pixels. + def test_set_and_get_low_gain_values(self): + pixel_ids = [1, 2, 3, 4, 5] + gain_maker = GainMakerforTest(pixel_ids) + low_gain_values = np.array([0.1, 0.2, 0.3, 0.4, 0.5]) + gain_maker._low_gain = low_gain_values + assert np.array_equal(gain_maker.low_gain, low_gain_values) + + # Tests that the results can be saved to a file. + def test_save_results_to_file(self, tmp_path=Path(f"/tmp/{np.random.rand()}")): + pixel_ids = [1, 2, 3, 4, 5] + gain_maker = GainMakerforTest(pixel_ids) + high_gain_values = np.array([0.5, 0.6, 0.7, 0.8, 0.9]) + gain_maker._high_gain = high_gain_values + low_gain_values = np.array([0.1, 0.2, 0.3, 0.4, 0.5]) + gain_maker._low_gain = low_gain_values + gain_maker.save(tmp_path) + assert (tmp_path / "results_test.ecsv").exists() + + # Tests that a copy of the result table can be obtained. + def test_get_copy_of_result_table(self): + pixel_ids = [1, 2, 3, 4, 5] + gain_maker = GainMakerforTest(pixel_ids) + result_table_copy = gain_maker.results + assert isinstance(result_table_copy, QTable) + assert np.array_equal( + result_table_copy[GainMakerforTest.PIXELS_ID_COLUMN], np.array(pixel_ids) + ) diff --git a/src/nectarchain/makers/calibration/gain/utils/__init__.py b/src/nectarchain/makers/calibration/gain/utils/__init__.py new file mode 100644 index 00000000..6f3aa3a4 --- /dev/null +++ b/src/nectarchain/makers/calibration/gain/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/makers/calibration/gain/utils/error.py new file mode 100644 index 00000000..9e1c9e75 --- /dev/null +++ b/src/nectarchain/makers/calibration/gain/utils/error.py @@ -0,0 +1,32 @@ +import logging + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + + +class DifferentPixelsID(Exception): + def __init__(self, message): + self.__message = message + + @property + def message(self): + return self.__message + + +class PedestalValueError(ValueError): + def __init__(self, message): + self.__message = message + + @property + def message(self): + return self.__message + + +class MeanValueError(ValueError): + def __init__(self, message): + self.__message = message + + @property + def message(self): + return self.__message diff --git a/src/nectarchain/calibration/NectarGain/SPEfit/utils.py b/src/nectarchain/makers/calibration/gain/utils/utils.py similarity index 54% rename from src/nectarchain/calibration/NectarGain/SPEfit/utils.py rename to src/nectarchain/makers/calibration/gain/utils/utils.py index 0dbb91f3..d5577bac 100644 --- a/src/nectarchain/calibration/NectarGain/SPEfit/utils.py +++ b/src/nectarchain/makers/calibration/gain/utils/utils.py @@ -5,35 +5,39 @@ from iminuit import Minuit from scipy import interpolate, signal from scipy.special import gammainc -from scipy.stats import norm,chi2 +from scipy.stats import chi2, norm -logging.basicConfig(format='%(asctime)s %(name)s %(levelname)s %(message)s') +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") log = logging.getLogger(__name__) -log.handlers = logging.getLogger('__main__').handlers +log.handlers = logging.getLogger("__main__").handlers -__all__ = ['UtilsMinuit','Multiprocessing'] -from .parameters import Parameters +# __all__ = ['UtilsMinuit','Multiprocessing'] -class Multiprocessing() : +from ..parameters import Parameters + + +class Multiprocessing: @staticmethod def custom_error_callback(error): - log.error(f'Got an error: {error}') - log.error(error,exc_info=True) + log.error(f"Got an error: {error}") + log.error(error, exc_info=True) -class Statistics() : + +class Statistics: @staticmethod - def chi2_pvalue(ndof : int, likelihood : float) : - return 1 - chi2(df = ndof).cdf(likelihood) + def chi2_pvalue(ndof: int, likelihood: float): + return 1 - chi2(df=ndof).cdf(likelihood) + -class UtilsMinuit() : +class UtilsMinuit: @staticmethod - def make_minuit_par_kwargs(parameters : Parameters): + def make_minuit_par_kwargs(parameters: Parameters): """Create *Parameter Keyword Arguments* for the `Minuit` constructor. updated for Minuit >2.0 """ names = parameters.parnames - kwargs = {"names": names,"values" : {}} + kwargs = {"names": names, "values": {}} for parameter in parameters.parameters: kwargs["values"][parameter.name] = parameter.value @@ -42,31 +46,38 @@ def make_minuit_par_kwargs(parameters : Parameters): error = 0.1 if np.isnan(parameter.error) else parameter.error kwargs[f"limit_{parameter.name}"] = (min_, max_) kwargs[f"error_{parameter.name}"] = error + if parameter.frozen: + kwargs[f"fix_{parameter.name}"] = True return kwargs - + @staticmethod - def set_minuit_parameters_limits_and_errors(m : Minuit,parameters : dict) : + def set_minuit_parameters_limits_and_errors(m: Minuit, parameters: dict): """function to set minuit parameter limits and errors with Minuit >2.0 Args: m (Minuit): a Minuit instance parameters (dict): dict containing parameters names, limits errors and values """ - for name in parameters["names"] : + for name in parameters["names"]: m.limits[name] = parameters[f"limit_{name}"] m.errors[name] = parameters[f"error_{name}"] + if parameters.get(f"fix_{name}", False): + m.fixed[name] = True # Usefull fucntions for the fit def gaussian(x, mu, sig): - #return (1./(sig*np.sqrt(2*math.pi)))*np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.))) - return norm.pdf(x,loc = mu,scale = sig) + # return (1./(sig*np.sqrt(2*math.pi)))*np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.))) + return norm.pdf(x, loc=mu, scale=sig) + -def weight_gaussian(x,N, mu, sig) : +def weight_gaussian(x, N, mu, sig): return N * gaussian(x, mu, sig) -def doubleGauss(x,sig1,mu2,sig2,p): - return p *2 *gaussian(x, 0, sig1) + (1-p) * gaussian(x, mu2, sig2) + +def doubleGauss(x, sig1, mu2, sig2, p): + return p * 2 * gaussian(x, 0, sig1) + (1 - p) * gaussian(x, mu2, sig2) + def PMax(r): """p_{max} in equation 6 in Caroff et al. (2019) @@ -77,13 +88,14 @@ def PMax(r): Returns: float : p_{max} """ - if r > np.sqrt((np.pi -2 )/ 2) : - pmax = np.pi/(2 * (r**2 + 1)) - else : - pmax = np.pi*r**2/(np.pi*r**2 + np.pi - 2*r**2 - 2) + if r > np.sqrt((np.pi - 2) / 2): + pmax = np.pi / (2 * (r**2 + 1)) + else: + pmax = np.pi * r**2 / (np.pi * r**2 + np.pi - 2 * r**2 - 2) return pmax -def ax(p,res): + +def ax(p, res): """a in equation 4 in Caroff et al. (2019) Args: @@ -93,9 +105,10 @@ def ax(p,res): Returns: float : a """ - return ((2/np.pi)*p**2-p/(res**2+1)) + return (2 / np.pi) * p**2 - p / (res**2 + 1) + -def bx(p,mu2): +def bx(p, mu2): """b in equation 4 in Caroff et al. (2019) Args: @@ -105,11 +118,12 @@ def bx(p,mu2): Returns: float : b """ - return (np.sqrt(2/np.pi)*2*p*(1-p)*mu2) + return np.sqrt(2 / np.pi) * 2 * p * (1 - p) * mu2 -def cx(sig2,mu2,res,p): + +def cx(sig2, mu2, res, p): """c in equation 4 in Caroff et al. (2019) - Note : There is a typo in the article 1-p**2 -> (1-p)**2 + Note : There is a typo in the article 1-p**2 -> (1-p)**2 Args: sig2 (float): width of the high charge Gaussian @@ -120,9 +134,10 @@ def cx(sig2,mu2,res,p): Returns: float : c """ - return (1-p)**2*mu2**2 - (1-p)*(sig2**2+mu2**2)/(res**2+1) + return (1 - p) ** 2 * mu2**2 - (1 - p) * (sig2**2 + mu2**2) / (res**2 + 1) + -def delta(p,res,sig2,mu2): +def delta(p, res, sig2, mu2): """well known delta in 2nd order polynom Args: @@ -134,9 +149,10 @@ def delta(p,res,sig2,mu2): Returns: float : b**2 - 4*a*c """ - return bx(p,mu2)*bx(p,mu2) - 4*ax(p,res)*cx(sig2,mu2,res,p) + return bx(p, mu2) * bx(p, mu2) - 4 * ax(p, res) * cx(sig2, mu2, res, p) + -def ParamU(p,r): +def ParamU(p, r): """d in equation 6 in Caroff et al. (2019) Args: @@ -146,9 +162,12 @@ def ParamU(p,r): Returns: float : d """ - return ((8*(1-p)**2*p**2)/np.pi - 4*(2*p**2/np.pi - p/(r**2+1))*((1-p)**2-(1-p)/(r**2+1))) + return (8 * (1 - p) ** 2 * p**2) / np.pi - 4 * ( + 2 * p**2 / np.pi - p / (r**2 + 1) + ) * ((1 - p) ** 2 - (1 - p) / (r**2 + 1)) -def ParamS(p,r): + +def ParamS(p, r): """e in equation 6 in Caroff et al. (2019) Args: @@ -158,10 +177,11 @@ def ParamS(p,r): Returns: float : e """ - e = (4*(2*p**2/np.pi - p/(r**2+1))*(1-p))/(r**2+1) + e = (4 * (2 * p**2 / np.pi - p / (r**2 + 1)) * (1 - p)) / (r**2 + 1) return e -def SigMin(p,res,mu2): + +def SigMin(p, res, mu2): """sigma_{high,min} in equation 6 in Caroff et al. (2019) Args: @@ -172,9 +192,12 @@ def SigMin(p,res,mu2): Returns: float : sigma_{high,min} """ - return mu2*np.sqrt((-ParamU(p,res)+(bx(p,mu2)**2/mu2**2))/(ParamS(p,res))) + return mu2 * np.sqrt( + (-ParamU(p, res) + (bx(p, mu2) ** 2 / mu2**2)) / (ParamS(p, res)) + ) -def SigMax(p,res,mu2): + +def SigMax(p, res, mu2): """sigma_{high,min} in equation 6 in Caroff et al. (2019) Args: @@ -185,15 +208,16 @@ def SigMax(p,res,mu2): Returns: float : sigma_{high,min} """ - temp = (-ParamU(p,res))/(ParamS(p,res)) - if temp < 0 : + temp = (-ParamU(p, res)) / (ParamS(p, res)) + if temp < 0: err = ValueError("-d/e must be < 0") - log.error(err,exc_info=True) + log.error(err, exc_info=True) raise err - else : - return mu2*np.sqrt(temp) + else: + return mu2 * np.sqrt(temp) -def sigma1(p,res,sig2,mu2): + +def sigma1(p, res, sig2, mu2): """sigma_{low} in equation 5 in Caroff et al. (2019) Args: @@ -205,9 +229,10 @@ def sigma1(p,res,sig2,mu2): Returns: float : sigma_{low} """ - return (-bx(p,mu2)+np.sqrt(delta(p,res,sig2,mu2)))/(2*ax(p,res)) + return (-bx(p, mu2) + np.sqrt(delta(p, res, sig2, mu2))) / (2 * ax(p, res)) + -def sigma2(n,p,res,mu2): +def sigma2(n, p, res, mu2): """sigma_{high} in equation 7 in Caroff et al. (2019) Args: @@ -219,39 +244,56 @@ def sigma2(n,p,res,mu2): Returns: float : sigma_{high} """ - if ((-ParamU(p,res)+(bx(p,mu2)**2/mu2**2))/(ParamS(p,res)) > 0): - return SigMin(p,res,mu2)+n*(SigMax(p,res,mu2)-SigMin(p,res,mu2)) + if (-ParamU(p, res) + (bx(p, mu2) ** 2 / mu2**2)) / (ParamS(p, res)) > 0: + return SigMin(p, res, mu2) + n * (SigMax(p, res, mu2) - SigMin(p, res, mu2)) else: - return n*SigMax(p,res,mu2) + return n * SigMax(p, res, mu2) # The real final model callign all the above for luminosity (lum) + PED, wil return probability of number of Spe -def MPE2(x,pp,res,mu2,n,muped,sigped,lum,**kwargs): - log.debug(f"pp = {pp}, res = {res}, mu2 = {mu2}, n = {n}, muped = {muped}, sigped = {sigped}, lum = {lum}") + + +def MPE2(x, pp, res, mu2, n, muped, sigped, lum, **kwargs): + log.debug( + f"pp = {pp}, res = {res}, mu2 = {mu2}, n = {n}, muped = {muped}, sigped = {sigped}, lum = {lum}" + ) f = 0 - ntotalPE = kwargs.get("ntotalPE",0) - if ntotalPE == 0 : - #about 1sec + ntotalPE = kwargs.get("ntotalPE", 0) + if ntotalPE == 0: + # about 1sec for i in range(1000): - if (gammainc(i+1,lum) < 1e-5): + if gammainc(i + 1, lum) < 1e-5: ntotalPE = i break - #print(ntotalPE) - #about 8 sec, 1 sec by nPEPDF call - #for i in range(ntotalPE): + # print(ntotalPE) + # about 8 sec, 1 sec by nPEPDF call + # for i in range(ntotalPE): # f = f + ((lum**i)/math.factorial(i)) * np.exp(-lum) * nPEPDF(x,pp,res,mu2,n,muped,sigped,i,int(mu2*ntotalPE+10*mu2)) - f = np.sum([(lum**i)/math.factorial(i) * np.exp(-lum) * nPEPDF(x,pp,res,mu2,n,muped,sigped,i,int(mu2*ntotalPE+10*mu2)) for i in range(ntotalPE)],axis = 0) # 10 % faster + f = np.sum( + [ + (lum**i) + / math.factorial(i) + * np.exp(-lum) + * nPEPDF( + x, pp, res, mu2, n, muped, sigped, i, int(mu2 * ntotalPE + 10 * mu2) + ) + for i in range(ntotalPE) + ], + axis=0, + ) # 10 % faster return f - + + # Fnal model shape/function (for one SPE) -def doubleGaussConstrained(x,pp,res,mu2,n): - p = pp*PMax(res) - sig2 = sigma2(n,p,res,mu2) - sig1 = sigma1(p,res,sig2,mu2) - return doubleGauss(x,sig1,mu2,sig2,p) +def doubleGaussConstrained(x, pp, res, mu2, n): + p = pp * PMax(res) + sig2 = sigma2(n, p, res, mu2) + sig1 = sigma1(p, res, sig2, mu2) + return doubleGauss(x, sig1, mu2, sig2, p) + # Get the gain from the parameters model -def Gain(pp,res,mu2,n): +def Gain(pp, res, mu2, n): """analytic gain computatuon Args: @@ -263,32 +305,35 @@ def Gain(pp,res,mu2,n): Returns: float : gain """ - p = pp*PMax(res) - sig2 = sigma2(n,p,res,mu2) - return (1-p)*mu2 + 2*p*sigma1(p,res,sig2,mu2)/np.sqrt(2*np.pi) + p = pp * PMax(res) + sig2 = sigma2(n, p, res, mu2) + return (1 - p) * mu2 + 2 * p * sigma1(p, res, sig2, mu2) / np.sqrt(2 * np.pi) + -def nPEPDF(x,pp,res,mu2,n,muped,sigped,nph,size_charge): - allrange = np.linspace(-1 * size_charge,size_charge,size_charge*2) +def nPEPDF(x, pp, res, mu2, n, muped, sigped, nph, size_charge): + allrange = np.linspace(-1 * size_charge, size_charge, size_charge * 2) spe = [] - #about 2 sec this is the main pb - #for i in range(len(allrange)): + # about 2 sec this is the main pb + # for i in range(len(allrange)): # if (allrange[i]>=0): # spe.append(doubleGaussConstrained(allrange[i],pp,res,mu2,n)) # else: # spe.append(0) - - spe = doubleGaussConstrained(allrange,pp,res,mu2,n) * (allrange>=0 * np.ones(allrange.shape)) #100 times faster + + spe = doubleGaussConstrained(allrange, pp, res, mu2, n) * ( + allrange >= 0 * np.ones(allrange.shape) + ) # 100 times faster # ~ plt.plot(allrange,spe) - #npe = semi_gaussian(allrange, muped, sigped) - npe = gaussian(allrange, 0, sigped) + # npe = semi_gaussian(allrange, muped, sigped) + npe = gaussian(allrange, 0, sigped) # ~ plt.plot(allrange,npe) # ~ plt.show() for i in range(nph): - #npe = np.convolve(npe,spe,"same") - npe = signal.fftconvolve(npe,spe,"same") + # npe = np.convolve(npe,spe,"same") + npe = signal.fftconvolve(npe, spe, "same") # ~ plt.plot(allrange,npe) # ~ plt.show() - fff = interpolate.UnivariateSpline(allrange,npe,ext=1,k=3,s=0) - norm = np.trapz(fff(allrange),allrange) - return fff(x-muped)/norm \ No newline at end of file + fff = interpolate.UnivariateSpline(allrange, npe, ext=1, k=3, s=0) + norm = np.trapz(fff(allrange), allrange) + return fff(x - muped) / norm diff --git a/src/nectarchain/makers/calibration/pedestalMakers.py b/src/nectarchain/makers/calibration/pedestalMakers.py new file mode 100644 index 00000000..15a8af1b --- /dev/null +++ b/src/nectarchain/makers/calibration/pedestalMakers.py @@ -0,0 +1,19 @@ +import logging + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + +from .core import CalibrationMaker + +__all__ = ["PedestalMaker"] + + +class PedestalMaker(CalibrationMaker): + def __init__(self, *args, **kwargs) -> None: + super().__init__(*args, **kwargs) + + def make(self): + raise NotImplementedError( + "The computation of the pedestal calibration is not yet implemented, feel free to contribute !:)" + ) diff --git a/src/nectarchain/makers/calibration/tests/test_core.py b/src/nectarchain/makers/calibration/tests/test_core.py new file mode 100644 index 00000000..d2a27fab --- /dev/null +++ b/src/nectarchain/makers/calibration/tests/test_core.py @@ -0,0 +1,58 @@ +from pathlib import Path + +import numpy as np +import pytest + +from nectarchain.makers.calibration.core import CalibrationMaker + + +class CalibrationMakerforTest(CalibrationMaker): + _reduced_name = "test" + + def make(): + pass + + +class TestflatfieldMaker: + # Tests that the constructor initializes the object with the correct attributes and metadata when valid input is provided + def test_constructor_with_valid_input(self): + pixels_id = [1, 2, 3] + calibration_maker = CalibrationMakerforTest(pixels_id) + + assert np.equal(calibration_maker._pixels_id, pixels_id).all() + assert np.equal( + calibration_maker._results[calibration_maker.PIXELS_ID_COLUMN], + np.array(pixels_id), + ).all() + assert calibration_maker._results.meta[calibration_maker.NP_PIXELS] == len( + pixels_id + ) + assert isinstance(calibration_maker._results.meta["comments"], str) + + # Tests that the constructor raises an error when a non-iterable pixels_id is provided + def test_constructor_with_non_iterable_pixels_id(self): + pixels_id = 123 + with pytest.raises(TypeError): + CalibrationMakerforTest(pixels_id) + + # Tests that saving the results to an existing file with overwrite=False raises an error + def test_save_to_existing_file_with_overwrite_false(self, tmp_path=Path("/tmp")): + pixels_id = [1, 2, 3] + calibration_maker = CalibrationMakerforTest(pixels_id) + + # Create a temporary file + file_path = tmp_path / "results_Calibration.ecsv" + file_path.touch() + + with pytest.raises(FileExistsError): + calibration_maker.save(file_path, overwrite=False) + + # Tests that changing the pixels_id attribute updates the results table with the expected values + def test_change_pixels_id_attribute(self): + pixels_id = [1, 2, 3] + calibration_maker = CalibrationMakerforTest(pixels_id) + + new_pixels_id = [4, 5, 6] + calibration_maker._pixels_id = np.array(new_pixels_id) + + assert np.equal(calibration_maker._pixels_id, new_pixels_id).all() diff --git a/src/nectarchain/makers/chargesMakers.py b/src/nectarchain/makers/chargesMakers.py new file mode 100644 index 00000000..f6d82d4c --- /dev/null +++ b/src/nectarchain/makers/chargesMakers.py @@ -0,0 +1,56 @@ +import logging + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + +import time +from argparse import ArgumentError + +import numpy as np +import numpy.ma as ma +from ctapipe.containers import EventType +from ctapipe.image.extractor import ( + BaselineSubtractedNeighborPeakWindowSum, + FixedWindowSum, + FullWaveformSum, + GlobalPeakWindowSum, + LocalPeakWindowSum, + NeighborPeakWindowSum, + SlidingWindowMaxSum, + TwoPassWindowSum, +) +from ctapipe.instrument.subarray 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, WaveformsContainer +from .extractor.utils import CtapipeExtractor + +import numpy as np + +from ctapipe.containers import EventType +from ctapipe.instrument import SubarrayDescription +from ctapipe_io_nectarcam import constants +from ctapipe_io_nectarcam.containers import NectarCAMDataContainer +from ctapipe.core.traits import ComponentNameList,Unicode +from tqdm import tqdm + +from ..data.container import WaveformsContainer +from .core import EventsLoopNectarCAMCalibrationTool +from .component import NectarCAMComponent,ChargesComponent,get_specific_traits + +__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/__init__.py b/src/nectarchain/makers/component/__init__.py new file mode 100644 index 00000000..dee30c9d --- /dev/null +++ b/src/nectarchain/makers/component/__init__.py @@ -0,0 +1,3 @@ +from .core import * +from .waveformsComponent import * +from .chargesComponent import * \ No newline at end of file diff --git a/src/nectarchain/makers/component/chargesComponent.py b/src/nectarchain/makers/component/chargesComponent.py new file mode 100644 index 00000000..7a5947e4 --- /dev/null +++ b/src/nectarchain/makers/component/chargesComponent.py @@ -0,0 +1,613 @@ +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 +import numpy.ma as ma + +from abc import ABC, abstractmethod +import copy +import tqdm +from argparse import ArgumentError +import time +from numba import bool_, float64, guvectorize, int64 + +from ctapipe.core import Component,TelescopeComponent +from ctapipe.instrument import CameraGeometry +from ctapipe.containers import EventType +from ctapipe.core.traits import List,Unicode,Dict +from ctapipe.instrument import SubarrayDescription +from ctapipe.image.extractor import ( + BaselineSubtractedNeighborPeakWindowSum, + FixedWindowSum, + FullWaveformSum, + GlobalPeakWindowSum, + LocalPeakWindowSum, + NeighborPeakWindowSum, + SlidingWindowMaxSum, + TwoPassWindowSum, +) +from ctapipe_io_nectarcam.containers import NectarCAMDataContainer +from ctapipe_io_nectarcam import NectarCAMEventSource, constants + +from .core import ArrayDataComponent +from ..extractor.utils import CtapipeExtractor +from ...data.container import ChargesContainer,ChargesContainers,WaveformsContainer + +__all__ = ["ChargesComponent"] + +list_ctapipe_charge_extractor = [ + "FullWaveformSum", + "FixedWindowSum", + "GlobalPeakWindowSum", + "LocalPeakWindowSum", + "SlidingWindowMaxSum", + "NeighborPeakWindowSum", + "BaselineSubtractedNeighborPeakWindowSum", + "TwoPassWindowSum", +] + + +list_nectarchain_charge_extractor = ["gradient_extractor"] + + +@guvectorize( + [ + (int64[:], float64[:], bool_, bool_[:], int64[:]), + ], + "(s),(n),()->(n),(n)", + nopython=True, + cache=True, +) +def make_histo(charge, all_range, mask_broken_pix, _mask, hist_ma_data): + """compute histogram of charge with numba + + Args: + charge (np.ndarray(pixels,nevents)): charge + all_range (np.ndarray(nbins)): charge range + mask_broken_pix (np.ndarray(pixels)): mask on broxen pixels + _mask (np.ndarray(pixels,nbins)): mask + hist_ma_data (np.ndarray(pixels,nbins)): histogram + """ + # print(f"charge.shape = {charge.shape[0]}") + # print(f"_mask.shape = {_mask.shape[0]}") + # print(f"_mask[0] = {_mask[0]}") + # print(f"hist_ma_data[0] = {hist_ma_data[0]}") + # print(f"mask_broken_pix = {mask_broken_pix}") + + if not (mask_broken_pix): + # print("this pixel is not broken, let's continue computation") + hist, _charge = np.histogram( + charge, + bins=np.arange( + np.uint16(np.min(charge)) - 1, np.uint16(np.max(charge)) + 2, 1 + ), + ) + # print(f"hist.shape[0] = {hist.shape[0]}") + # print(f"charge.shape[0] = {_charge.shape[0]}") + charge_edges = np.array( + [np.mean(_charge[i : i + 2]) for i in range(_charge.shape[0] - 1)] + ) + # print(f"charge_edges.shape[0] = {charge_edges.shape[0]}") + mask = (all_range >= charge_edges[0]) * (all_range <= charge_edges[-1]) + # print(f"all_range = {int(all_range[0])}-{int(all_range[-1])}") + # print(f"charge_edges[0] = {int(charge_edges[0])}") + # print(f"charge_edges[-1] = {int(charge_edges[-1])}") + # print(f"mask[0] = {mask[0]}") + # print(f"mask[-1] = {mask[-1]}") + + # MASK THE DATA + # print(f"mask.shape = {mask.shape[0]}") + _mask[:] = ~mask + # print(f"_mask[0] = {_mask[0]}") + # print(f"_mask[-1] = {_mask[-1]}") + # FILL THE DATA + hist_ma_data[mask] = hist + # print("work done") + else: + # print("this pixel is broken, skipped") + pass + + + + + + + +class ChargesComponent(ArrayDataComponent) : + + 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) + + def __init__(self, subarray, config=None, parent=None,*args, **kwargs): + super().__init__(subarray = subarray, config = config, parent = parent, *args,**kwargs) + + self.__charges_hg = {} + self.__charges_lg = {} + self.__peak_hg = {} + self.__peak_lg = {} + + + def _init_trigger_type(self, trigger_type: EventType, **kwargs): + """ + Initializes the ChargesMaker based on the trigger type. + + Args: + trigger_type (EventType): The type of trigger. + **kwargs: Additional keyword arguments. + + Returns: + None + """ + super()._init_trigger_type(trigger_type, **kwargs) + name = __class__._get_name_trigger(trigger_type) + log.info(f"initialization of the ChargesMaker following trigger type : {name}") + self.__charges_hg[f"{name}"] = [] + self.__charges_lg[f"{name}"] = [] + self.__peak_hg[f"{name}"] = [] + self.__peak_lg[f"{name}"] = [] + + def __call__( + self, + event: NectarCAMDataContainer, + *args, + **kwargs, + ): + + wfs_hg_tmp, wfs_lg_tmp = super(ChargesComponent,self).__call__( + event=event, return_wfs=True, *args, **kwargs + ) + name = __class__._get_name_trigger(event.trigger.event_type) + + broken_pixels_hg, broken_pixels_lg = __class__._compute_broken_pixels_event( + event, self._pixels_id + ) + self._broken_pixels_hg[f"{name}"].append(broken_pixels_hg.tolist()) + self._broken_pixels_lg[f"{name}"].append(broken_pixels_lg.tolist()) + + imageExtractor = __class__._get_imageExtractor( + self.method, self.subarray, **self.extractor_kwargs + ) + + __image = CtapipeExtractor.get_image_peak_time( + imageExtractor( + wfs_hg_tmp, self.TEL_ID, constants.HIGH_GAIN, broken_pixels_hg + ) + ) + self.__charges_hg[f"{name}"].append(__image[0].tolist()) + self.__peak_hg[f"{name}"].append(__image[1].tolist()) + + __image = CtapipeExtractor.get_image_peak_time( + imageExtractor( + wfs_lg_tmp, self.TEL_ID, constants.LOW_GAIN, broken_pixels_lg + ) + ) + self.__charges_lg[f"{name}"].append(__image[0].tolist()) + self.__peak_lg[f"{name}"].append(__image[1].tolist()) + + @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. + 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. + Returns: + 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 {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 + ) + log.debug( + f"Extracting charges with method {method} and extractor_kwargs {extractor_kwargs}" + ) + imageExtractor = eval(method)(subarray, **extractor_kwargs) + return imageExtractor + + def finish( + self, *args, **kwargs + ): + """ + Create an output container for the specified trigger type and method. + Args: + trigger_type (EventType): The type of trigger. + method (str): The name of the charge extraction method. + *args: Additional positional arguments. + **kwargs: Additional keyword arguments. + Returns: + list: A list of ChargesContainer objects. + """ + output = ChargesContainers() + for i,trigger in enumerate(self.trigger_list): + chargesContainer = ChargesContainer( + run_number=self._run_number, + npixels=self._npixels, + camera=self.CAMERA_NAME, + pixels_id=self._pixels_id, + method=self.method, + nevents=self.nevents(trigger), + charges_hg=self.charges_hg(trigger), + charges_lg=self.charges_lg(trigger), + peak_hg=self.peak_hg(trigger), + peak_lg=self.peak_lg(trigger), + broken_pixels_hg=self.broken_pixels_hg(trigger), + broken_pixels_lg=self.broken_pixels_lg(trigger), + ucts_timestamp=self.ucts_timestamp(trigger), + ucts_busy_counter=self.ucts_busy_counter(trigger), + ucts_event_counter=self.ucts_event_counter(trigger), + event_type=self.event_type(trigger), + event_id=self.event_id(trigger), + trig_pattern_all=self.trig_pattern_all(trigger), + trig_pattern=self.trig_pattern(trigger), + multiplicity=self.multiplicity(trigger), + ) + output.containers[trigger] = chargesContainer + return output + + @staticmethod + 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. + 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. + + Raises: + ArgumentError: If the specified method is not valid. + """ + output = ChargesContainer( + run_number=chargesContainer.run_number, + npixels=chargesContainer.npixels, + camera=chargesContainer.camera, + pixels_id=chargesContainer.pixels_id, + nevents=chargesContainer.nevents, + method=chargesContainer.method, + ) + if method == "event_id": + index = np.argsort(chargesContainer.event_id) + for field in chargesContainer.keys(): + if not ( + field + in [ + "run_number", + "npixels", + "camera", + "pixels_id", + "nevents", + "method", + ] + ): + output[field] = chargesContainer[field][index] + else: + raise ArgumentError(f"{method} is not a valid method for sorting") + return output + + @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. + 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. + """ + res = __class__.select_container_array_field( + container=chargesContainer, pixel_id=pixel_id, field="charges_hg" + ) + res = res.transpose(1, 0) + return res + + @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. + 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. + """ + res = __class__.select_container_array_field( + container=chargesContainer, pixel_id=pixel_id, field="charges_lg" + ) + res = res.transpose(1, 0) + return res + + def charges_hg(self, trigger: EventType): + """ + 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 + ) + + def charges_lg(self, trigger: EventType): + """ + 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 + ) + + def peak_hg(self, trigger: EventType): + """ + 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 + ) + + def peak_lg(self, trigger: EventType): + """ + 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 + ) + + @staticmethod + def create_from_waveforms( + waveformsContainer: WaveformsContainer, + method: str = "FullWaveformSum", + **kwargs, + ) -> ChargesContainer: + """ + 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. + Returns: + 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 = 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 = np.array(charges_lg, dtype=np.uint16) + chargesContainer.charges_hg = charges_hg + chargesContainer.charges_lg = charges_lg + chargesContainer.peak_hg = peak_hg + chargesContainer.peak_lg = peak_lg + chargesContainer.method = method + return chargesContainer + + @staticmethod + def compute_charge( + waveformContainer: WaveformsContainer, + channel: int, + subarray : SubarrayDescription, + method: str = "FullWaveformSum", + tel_id : int = None, + **kwargs, + ): + """ + Compute charge from waveforms. + 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. + 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...)) + from ..extractor.utils import CtapipeExtractor + + if tel_id is None : + tel_id = __class__.TEL_ID.default_value + + imageExtractor = __class__._get_imageExtractor( + method=method, subarray=subarray, **kwargs + ) + if channel == constants.HIGH_GAIN: + out = np.array( + [ + CtapipeExtractor.get_image_peak_time( + imageExtractor( + waveformContainer.wfs_hg[i], + tel_id, + channel, + waveformContainer.broken_pixels_hg, + ) + ) + for i in range(len(waveformContainer.wfs_hg)) + ] + ).transpose(1, 0, 2) + return out[0], out[1] + elif channel == constants.LOW_GAIN: + out = np.array( + [ + CtapipeExtractor.get_image_peak_time( + imageExtractor( + waveformContainer.wfs_lg[i], + tel_id, + channel, + waveformContainer.broken_pixels_lg, + ) + ) + for i in range(len(waveformContainer.wfs_lg)) + ] + ).transpose(1, 0, 2) + return out[0], out[1] + else: + raise ArgumentError( + f"channel must be {constants.LOW_GAIN} or {constants.HIGH_GAIN}" + ) + + @staticmethod + def histo_hg( + chargesContainer: ChargesContainer, n_bins: int = 1000, autoscale: bool = True + ) -> ma.masked_array: + """ + 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. + + 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. + """ + return __class__._histo( + chargesContainer=chargesContainer, + field="charges_hg", + n_bins=n_bins, + autoscale=autoscale, + ) + + @staticmethod + def histo_lg( + chargesContainer: ChargesContainer, n_bins: int = 1000, autoscale: bool = True + ) -> ma.masked_array: + """ + 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. + + 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. + """ + return __class__._histo( + chargesContainer=chargesContainer, + field="charges_lg", + n_bins=n_bins, + autoscale=autoscale, + ) + + @staticmethod + def _histo( + chargesContainer: ChargesContainer, + field: str, + n_bins: int = 1000, + autoscale: bool = True, + ) -> ma.masked_array: + """ + Computes histogram of charges for a given field from a ChargesContainer object. + 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. + 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. + + 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. + """ + mask_broken_pix = np.array( + (chargesContainer[field] == chargesContainer[field].mean(axis=0)).mean( + axis=0 + ), + dtype=bool, + ) + log.debug( + f"there are {mask_broken_pix.sum()} broken pixels (charge stays at same level for each events)" + ) + + if autoscale: + all_range = np.arange( + np.uint16(np.min(chargesContainer[field].T[~mask_broken_pix].T)) - 0.5, + np.uint16(np.max(chargesContainer[field].T[~mask_broken_pix].T)) + 1.5, + 1, + ) + charge_ma = ma.masked_array( + ( + all_range.reshape(all_range.shape[0], 1) + @ np.ones((1, chargesContainer[field].shape[1])) + ).T, + mask=np.zeros( + (chargesContainer[field].shape[1], all_range.shape[0]), dtype=bool + ), + ) + broxen_pixels_mask = np.array( + [mask_broken_pix for i in range(charge_ma.shape[1])] + ).T + start = time.time() + _mask, hist_ma_data = make_histo( + chargesContainer[field].T, all_range, mask_broken_pix + ) + charge_ma.mask = np.logical_or(_mask, broxen_pixels_mask) + hist_ma = ma.masked_array(hist_ma_data, mask=charge_ma.mask) + log.debug(f"histogram hg computation time : {time.time() - start} sec") + + return ma.masked_array((hist_ma, charge_ma)) + + else: + hist = np.array( + [ + np.histogram(chargesContainer[field].T[i], bins=n_bins)[0] + for i in range(chargesContainer[field].shape[1]) + ] + ) + charge = np.array( + [ + np.histogram(chargesContainer[field].T[i], bins=n_bins)[1] + for i in range(chargesContainer[field].shape[1]) + ] + ) + charge_edges = np.array( + [ + np.mean(charge.T[i : i + 2], axis=0) + for i in range(charge.shape[1] - 1) + ] + ).T + + return np.array((hist, charge_edges)) diff --git a/src/nectarchain/makers/component/core.py b/src/nectarchain/makers/component/core.py new file mode 100644 index 00000000..c6fb5eb2 --- /dev/null +++ b/src/nectarchain/makers/component/core.py @@ -0,0 +1,505 @@ +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 abc import ABC, abstractmethod +import copy + +from ctapipe.core import Component,TelescopeComponent +from ctapipe.instrument import CameraGeometry +from ctapipe.containers import EventType +from ctapipe.core.traits import Unicode,Integer + + +from ctapipe_io_nectarcam.containers import NectarCAMDataContainer +from ctapipe_io_nectarcam import NectarCAMEventSource, constants + +from ...data.container.core import ArrayDataContainer + + +__all__ = ["ArrayDataComponent", + "NectarCAMComponent", + "get_valid_component", + "get_specific_traits", + "get_configurable_traits", + ] + +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""" + def __init__(self, subarray, config=None, parent=None,*args, **kwargs): + super().__init__(subarray = subarray, config = config, parent = parent, *args,**kwargs) + self.__pixels_id = parent._event_source.camera_config.expected_pixels_id + self.__run_number = parent.run_number + self.__npixels=parent.npixels + @abstractmethod + def __call__( + self, event: NectarCAMDataContainer, *args, **kwargs + ): + pass + + @property + def _pixels_id(self) : return self.__pixels_id + @property + def pixels_id(self) : return copy.deepcopy(self.__pixels_id) + @property + def _run_number(self) : + return self.__run_number + @property + def run_number(self) : + return copy.deepcopy(self.__run_number) + @property + def _npixels(self) : + return self.__npixels + @property + 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, + help = "The telescope ID", + read_only = True, + ).tag(config = True) + + CAMERA_NAME = Unicode(default_value = "NectarCam-003", + help = "The camera name", + read_only = True, + ).tag(config = True) + + CAMERA = CameraGeometry.from_name(CAMERA_NAME.default_value) + + #trigger_list = List( + # help="List of trigger(EventType) inside the instance", + # default_value=[], + #).tag(config=True) + + + + def __init__(self, subarray, config=None, parent=None,*args, **kwargs): + super().__init__(subarray = subarray,config = config, parent = parent,*args, **kwargs) + self.__nsamples = parent._event_source.camera_config.num_samples + + + self.trigger_list = [] + + # data we want to compute + self.__ucts_timestamp = {} + self.__ucts_busy_counter = {} + self.__ucts_event_counter = {} + self.__event_type = {} + self.__event_id = {} + self.__trig_patter_all = {} + self.__broken_pixels_hg = {} + self.__broken_pixels_lg = {} + + 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. + + Returns: + None. The method only initializes the empty lists for the trigger type. + """ + name = __class__._get_name_trigger(trigger) + self.__ucts_timestamp[f"{name}"] = [] + self.__ucts_busy_counter[f"{name}"] = [] + self.__ucts_event_counter[f"{name}"] = [] + self.__event_type[f"{name}"] = [] + self.__event_id[f"{name}"] = [] + self.__trig_patter_all[f"{name}"] = [] + self.__broken_pixels_hg[f"{name}"] = [] + self.__broken_pixels_lg[f"{name}"] = [] + self.trigger_list.append(trigger) + + @staticmethod + def _compute_broken_pixels(wfs_hg, wfs_lg, **kwargs): + """ + Computes broken pixels for high and low gain waveforms. + Args: + wfs_hg (ndarray): High gain waveforms. + 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. + """ + log.warning("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 + ) + + @staticmethod + def _compute_broken_pixels_event( + event: NectarCAMDataContainer, pixels_id: np.ndarray, **kwargs + ): + """ + Computes broken pixels for a specific event and pixel IDs. + Args: + event (NectarCAMDataContainer): An event. + pixels_id (list or np.ndarray): IDs of pixels. + **kwargs: Additional keyword arguments. + Returns: + tuple: Two arrays of zeros with the length of `pixels_id`. + """ + log.warning("computation of broken pixels is not yet implemented") + return np.zeros((len(pixels_id)), dtype=bool), np.zeros( + (len(pixels_id)), dtype=bool + ) + + @staticmethod + def _get_name_trigger(trigger: EventType): + """ + Gets the name of a trigger event. + Args: + trigger (EventType): A trigger event. + Returns: + str: The name of the trigger event. + """ + if trigger is None: + name = "None" + else: + name = trigger.name + return name + + def __call__( + self, event: NectarCAMDataContainer, *args, **kwargs + ): + """ + Method to extract data from the event. + + Args: + event (NectarCAMDataContainer): The event object. + trigger (EventType): The trigger type. + *args: Additional arguments that can be passed to the method. + **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. + """ + name = __class__._get_name_trigger(event.trigger.event_type) + + if not(name in self.__event_id.keys()) : + self._init_trigger_type(event.trigger.event_type) + + self.__event_id[f"{name}"].append(np.uint32(event.index.event_id)) + self.__ucts_timestamp[f"{name}"].append( + event.nectarcam.tel[__class__.TEL_ID.default_value].evt.ucts_timestamp + ) + self.__event_type[f"{name}"].append(event.trigger.event_type.value) + self.__ucts_busy_counter[f"{name}"].append( + event.nectarcam.tel[__class__.TEL_ID.default_value].evt.ucts_busy_counter + ) + self.__ucts_event_counter[f"{name}"].append( + event.nectarcam.tel[__class__.TEL_ID.default_value].evt.ucts_event_counter + ) + self.__trig_patter_all[f"{name}"].append( + event.nectarcam.tel[__class__.TEL_ID.default_value].evt.trigger_pattern.T + ) + + if kwargs.get("return_wfs", False): + get_wfs_hg = event.r0.tel[0].waveform[constants.HIGH_GAIN][self.pixels_id] + get_wfs_lg = event.r0.tel[0].waveform[constants.LOW_GAIN][self.pixels_id] + return get_wfs_hg, get_wfs_lg + + + + @abstractmethod + def finish(self): + pass + + @staticmethod + 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. + + 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. + field (str): The name of the field to be selected from the container. + + Returns: + ndarray: An array containing the selected data for the given pixel IDs. + """ + mask_contain_pixels_id = np.array( + [pixel in container.pixels_id for pixel in pixel_id], dtype=bool + ) + 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 container, skip this one" + ) + res = np.array( + [ + np.take( + container[field], + np.where(container.pixels_id == pixel)[0][0], + axis=1, + ) + 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])) + return res + + @staticmethod + def merge( + container_a: ArrayDataContainer, container_b: ArrayDataContainer + ) -> ArrayDataContainer: + """method to merge 2 ArrayDataContainer into one single ArrayDataContainer + + Returns: + ArrayDataContainer: the merged object + """ + 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): + raise Exception("The containers have not the same pixels ids") + + merged_container = container_a.__class__.__new__() + + for field in container_a.keys(): + if ~isinstance(container_a[field], np.ndarray): + if container_a[field] != container_b[field]: + raise Exception( + 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 + ) + else: + merged_container[field] = container_a[field] + + return merged_container + + @property + def nsamples(self): + """ + Returns a deep copy of the nsamples attribute. + + Returns: + np.ndarray: A deep copy of the nsamples attribute. + """ + return copy.deepcopy(self.__nsamples) + + @property + def _nsamples(self): + """ + Returns the nsamples attribute. + + Returns: + np.ndarray: The nsamples attribute. + """ + return self.__nsamples + + + 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. + + Returns: + int: The number of events for the specified trigger type. + """ + return len(self.__event_id[__class__._get_name_trigger(trigger)]) + + @property + def _broken_pixels_hg(self): + """ + Returns the broken_pixels_hg attribute. + + Returns: + np.ndarray: The broken_pixels_hg attribute. + """ + return self.__broken_pixels_hg + + 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. + + Returns: + 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 + ) + + @property + def _broken_pixels_lg(self): + """ + Returns the broken_pixels_lg attribute. + + Returns: + np.ndarray: The broken_pixels_lg attribute. + """ + return self.__broken_pixels_lg + + 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. + + Returns: + 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 + ) + + 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. + + 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 + ) + + 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. + + 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, + ) + + 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. + + 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, + ) + + 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. + + 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 + ) + + def event_id(self, trigger: EventType): + """ + Returns an array of event IDs for the specified trigger type. + + Args: + trigger (EventType): The trigger type for which the event IDs are requested. + + Returns: + 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 + ) + + 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. + + Returns: + np.ndarray: An array of multiplicities for the specified trigger type. + """ + tmp = self.trig_pattern(trigger) + if len(tmp) == 0: + return np.array([]) + else: + return np.uint16(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. + + Returns: + np.ndarray: An array of trigger patterns for the specified trigger type. + """ + tmp = self.trig_pattern_all(trigger) + if len(tmp) == 0: + return np.array([]) + else: + return tmp.any(axis=2) + + def trig_pattern_all(self, trigger: EventType): + """ + 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. + + Returns: + 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, + ) + + \ No newline at end of file diff --git a/src/nectarchain/makers/component/waveformsComponent.py b/src/nectarchain/makers/component/waveformsComponent.py new file mode 100644 index 00000000..6874ccfb --- /dev/null +++ b/src/nectarchain/makers/component/waveformsComponent.py @@ -0,0 +1,326 @@ +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 abc import ABC, abstractmethod +import copy +import tqdm +from argparse import ArgumentError + + +from ctapipe.core import Component,TelescopeComponent +from ctapipe.instrument import CameraGeometry +from ctapipe.containers import EventType +from ctapipe.core.traits import List +from ctapipe.instrument import SubarrayDescription + + + + +from ctapipe_io_nectarcam.containers import NectarCAMDataContainer +from ctapipe_io_nectarcam import NectarCAMEventSource, constants + +from .core import ArrayDataComponent +from ...data.container import WaveformsContainer,WaveformsContainers + +__all__ = ["WaveformsComponent"] + +class WaveformsComponent(ArrayDataComponent) : + + def __init__(self, subarray, config=None, parent=None,*args, **kwargs): + super().__init__(subarray = subarray, config = config, parent = parent, *args,**kwargs) + + self.__geometry = subarray.tel[self.TEL_ID].camera + self.__wfs_hg = {} + self.__wfs_lg = {} + + @staticmethod + def create_from_events_list( + events_list: list, + run_number: int, + npixels: int, + nsamples: int, + subarray: SubarrayDescription, + pixels_id: int, + tel_id : int = None, + ) -> WaveformsContainer: + """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. + 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. + subarray (SubarrayDescription): The subarray description instance. + 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. + """ + if tel_id is None : + tel_id = __class__.TEL_ID.default_value + + container = WaveformsContainer( + run_number=run_number, + npixels=npixels, + nsamples=nsamples, + subarray=subarray, + camera=__class__.CAMERA_NAME, + pixels_id=pixels_id, + ) + + ucts_timestamp = [] + ucts_busy_counter = [] + ucts_event_counter = [] + event_type = [] + event_id = [] + trig_pattern_all = [] + wfs_hg = [] + wfs_lg = [] + + for event in tqdm(events_list): + ucts_timestamp.append( + event.nectarcam.tel[tel_id].evt.ucts_timestamp + ) + ucts_busy_counter.append( + event.nectarcam.tel[tel_id].evt.ucts_busy_counter + ) + ucts_event_counter.append( + event.nectarcam.tel[tel_id].evt.ucts_event_counter + ) + event_type.append(event.trigger.event_type.value) + event_id.append(event.index.event_id) + trig_pattern_all.append( + event.nectarcam.tel[tel_id].evt.trigger_pattern.T + ) + wfs_hg.append(event.r0.tel[0].waveform[constants.HIGH_GAIN][pixels_id]) + wfs_lg.append(event.r0.tel[0].waveform[constants.HIGH_GAIN][pixels_id]) + + container.wfs_hg = np.array(wfs_hg, dtype=np.uint16) + container.wfs_lg = np.array(wfs_lg, dtype=np.uint16) + + container.ucts_timestamp = np.array(ucts_timestamp, dtype=np.uint64) + container.ucts_busy_counter = np.array(ucts_busy_counter, dtype=np.uint32) + container.ucts_event_counter = np.array(ucts_event_counter, dtype=np.uint32) + container.event_type = np.array(event_type, dtype=np.uint8) + container.event_id = np.array(event_id, dtype=np.uint32) + container.trig_pattern_all = np.array(trig_pattern_all, dtype=bool) + container.trig_pattern = container.trig_pattern_all.any(axis=2) + container.multiplicity = np.uint16( + np.count_nonzero(container.trig_pattern, axis=1) + ) + + broken_pixels = __class__._compute_broken_pixels( + container.wfs_hg, container.wfs_lg + ) + container.broken_pixels_hg = broken_pixels[0] + container.broken_pixels_lg = broken_pixels[1] + return container + + def _init_trigger_type(self, trigger_type: EventType, **kwargs): + """Initialize the waveformsMaker following the trigger type. + + Args: + trigger_type: The type of trigger. + + """ + super()._init_trigger_type(trigger_type, **kwargs) + name = __class__._get_name_trigger(trigger_type) + log.info( + f"initialization of the waveformsMaker following trigger type : {name}" + ) + self.__wfs_hg[f"{name}"] = [] + self.__wfs_lg[f"{name}"] = [] + + def __call__( + self, event: NectarCAMDataContainer, *args, **kwargs + ): + """Process an event and extract waveforms. + + Args: + event (NectarCAMDataContainer): The event to process and extract waveforms from. + trigger (EventType): The type of trigger for the event. + + """ + wfs_hg_tmp = np.zeros((self.npixels, self.nsamples), dtype=np.uint16) + wfs_lg_tmp = np.zeros((self.npixels, self.nsamples), dtype=np.uint16) + + wfs_hg_tmp, wfs_lg_tmp = super(WaveformsComponent,self).__call__( + event=event, return_wfs=True, *args, **kwargs + ) + name = __class__._get_name_trigger(event.trigger.event_type) + + self.__wfs_hg[f"{name}"].append(wfs_hg_tmp.tolist()) + self.__wfs_lg[f"{name}"].append(wfs_lg_tmp.tolist()) + + broken_pixels_hg, broken_pixels_lg = __class__._compute_broken_pixels_event( + event, self._pixels_id + ) + self._broken_pixels_hg[f"{name}"].append(broken_pixels_hg.tolist()) + self._broken_pixels_lg[f"{name}"].append(broken_pixels_lg.tolist()) + + def finish(self, *args, **kwargs): + """Make the output container for the selected trigger types. + + Args: + trigger_type (EventType): The selected trigger types. + + Returns: + 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, + #subarray=self.subarray, + camera=self.CAMERA_NAME, + pixels_id=self._pixels_id, + nevents=self.nevents(trigger), + wfs_hg=self.wfs_hg(trigger), + wfs_lg=self.wfs_lg(trigger), + broken_pixels_hg=self.broken_pixels_hg(trigger), + broken_pixels_lg=self.broken_pixels_lg(trigger), + ucts_timestamp=self.ucts_timestamp(trigger), + ucts_busy_counter=self.ucts_busy_counter(trigger), + ucts_event_counter=self.ucts_event_counter(trigger), + event_type=self.event_type(trigger), + event_id=self.event_id(trigger), + trig_pattern_all=self.trig_pattern_all(trigger), + trig_pattern=self.trig_pattern(trigger), + multiplicity=self.multiplicity(trigger), + ) + output.containers[trigger] = waveformsContainer + return output + + @staticmethod + def sort(waveformsContainer: WaveformsContainer, method: str = "event_id"): + """Sort the waveformsContainer based on a specified method. + + Args: + waveformsContainer (WaveformsContainer): The waveformsContainer to be sorted. + method (str, optional): The sorting method. Defaults to 'event_id'. + + Returns: + WaveformsContainer: The sorted waveformsContainer. + """ + output = WaveformsContainer( + 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, + ) + if method == "event_id": + index = np.argsort(waveformsContainer.event_id) + for field in waveformsContainer.keys(): + if not ( + field + in [ + "run_number", + "npixels", + "nsamples", + "subarray", + "camera", + "pixels_id", + "nevents", + ] + ): + output[field] = waveformsContainer[field][index] + else: + raise ArgumentError(f"{method} is not a valid method for sorting") + return output + + @staticmethod + def select_waveforms_hg( + waveformsContainer: WaveformsContainer, pixel_id: np.ndarray + ): + """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. + + Returns: + np.ndarray: An array of selected waveforms from the container. + """ + res = __class__.select_container_array_field( + container=waveformsContainer, pixel_id=pixel_id, field="wfs_lg" + ) + res = res.transpose(1, 0, 2) + return res + + @staticmethod + def select_waveforms_lg( + waveformsContainer: WaveformsContainer, pixel_id: np.ndarray + ): + """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. + + Returns: + np.ndarray: An array of selected waveforms from the container. + """ + res = __class__.select_container_array_field( + container=waveformsContainer, pixel_id=pixel_id, field="wfs_hg" + ) + res = res.transpose(1, 0, 2) + return res + + @property + def _geometry(self): + """ + Returns the private __geometry attribute of the WaveformsMaker class. + + :return: The value of the private __geometry attribute. + """ + return self.__geometry + + + @property + def geometry(self): + """ + Returns a deep copy of the geometry attribute. + + Returns: + A deep copy of the geometry attribute. + """ + return copy.deepcopy(self.__geometry) + + + 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. + + 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 + ) + + def wfs_lg(self, trigger: EventType): + """ + 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. + + Returns: + 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 + ) diff --git a/src/nectarchain/makers/core.py b/src/nectarchain/makers/core.py new file mode 100644 index 00000000..62a53a25 --- /dev/null +++ b/src/nectarchain/makers/core.py @@ -0,0 +1,345 @@ +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 pathlib +from abc import ABC, abstractmethod + +import numpy as np +from ctapipe.containers import EventType +from ctapipe.core import Tool,Component +from ctapipe.core.traits import Bool, Integer, Path, classes_with_traits, flag, ComponentNameList +from ctapipe.instrument import CameraGeometry +from ctapipe.io import HDF5TableWriter +from ctapipe.io.datawriter import DATA_MODEL_VERSION +from ctapipe_io_nectarcam import NectarCAMEventSource, constants +from ctapipe_io_nectarcam.containers import NectarCAMDataContainer + + + +from tqdm.auto import tqdm + +from ..data import DataManagement +from ..data.container.core import ArrayDataContainer,NectarCAMContainer,TriggerMapContainer +from .component import * + +__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. +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.""" + + name = "BaseNectarCAMCalibration" + + progress_bar = Bool( + help="show progress bar during processing", default_value=False + ).tag(config=True) + + @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 + + 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 everythings. + run_file (optional) : if provided, will load this run file + Returns: + List[ctapipe_io_nectarcam.NectarCAMEventSource]: 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( + input_url=generic_filename, max_events=max_events + ) + else: + log.info(f"{run_file} will be loaded") + eventsource = NectarCAMEventSource( + input_url=run_file, max_events=max_events + ) + return eventsource + + +class EventsLoopNectarCAMCalibrationTool(BaseNectarCAMCalibrationTool): + """ + A class for data processing and computation on events from a specific run. + + 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. + run_file (optional): The specific run file to be loaded. + + Example Usage: + maker = EventsLoopMaker(run_number=1234, max_events=1000) + maker.make(n_events=500) + + + """ + + name = "EventsLoopNectarCAMCalibration" + + description = ( + __doc__ + f" This currently uses data model version {DATA_MODEL_VERSION}" + ) + examples = """To be implemented""" + + + aliases = { + ("i", "input"): "EventsLoopNectarCAMCalibrationTool.run_file", + ("r", "run-number"): "EventsLoopNectarCAMCalibrationTool.run_number", + ("m", "max-events"): "EventsLoopNectarCAMCalibrationTool.max_events", + ("o", "output"): "EventsLoopNectarCAMCalibrationTool.output_path", + } + + flags = { + "overwrite": ( + {"HDF5TableWriter": {"overwrite": True}}, + "Overwrite output file if it exists", + ), + **flag( + "progress", + "ProcessorTool.progress_bar", + "show a progress bar during event processing", + "don't show a progress bar during event processing", + ), + } + + classes = ( + [ + HDF5TableWriter, + ] + + classes_with_traits(NectarCAMEventSource) + + 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 + ) + + max_events = Integer( + help="maximum number of events to be loaded", + default_value=None, + allow_none=True, + ).tag(config=True) + + run_file = Path( + help="file name to be loaded", + default_value=None, + allow_none=True, + ).tag(config=True) + + componentsList = ComponentNameList(NectarCAMComponent, + help="List of Component names to be apply, the order will be respected" + ).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 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") + for componentName in _cls.componentsList : + configurable_traits = 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 _load_eventsource(self): + 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)) + 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 setup(self, *args, **kwargs): + if self.run_number == -1: + raise Exception("run_number need to be set up") + self._load_eventsource() + self.__npixels = self._event_source.camera_config.num_pixels + self.__pixels_id = self._event_source.camera_config.expected_pixels_id + + self.components = [] + for componentName in self.componentsList : + if componentName in get_valid_component(): + component_kwargs = self._get_provided_component_kwargs(componentName) + self.components.append( + # self.add_component( + Component.from_name( + componentName, + subarray = self.event_source.subarray, + parent=self, + **component_kwargs, + ) + # ) + ) + + self.writer = self.enter_context( + HDF5TableWriter( + filename = pathlib.Path(f"{self.output_path.parent}/{self.output_path.stem}_{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, + #trigger_type: list = None, + restart_from_begining: bool = False, + *args, + **kwargs, + ): + """ + 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. + *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): + self.log.warning( + "no needed events number specified, it may cause a memory error" + ) + #if isinstance(trigger_type, EventType) or trigger_type is None: + # trigger_type = [trigger_type] + #for _trigger_type in trigger_type: + # self._init_trigger_type(_trigger_type) + + if restart_from_begining: + self.log.debug("restart from begining : creation of the EventSource reader") + self._load_eventsource() + + n_traited_events = 0 + for i, event in enumerate( + tqdm( + self._event_source, + desc=self._event_source.__class__.__name__, + total=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}") + for component in self.components : + component(event,*args,**kwargs) + n_traited_events += 1 + if n_traited_events >= n_events: + break + + def finish(self, *args, **kwargs): + # self.write = self.enter_context( + # HDF5TableWriter(filename=filename, parent=self) + # ) + output = [] + for component in self.components : + output.append(component.finish(*args,**kwargs)) + log.info(output) + for _output in output : + if isinstance(_output,NectarCAMContainer) : + self.writer.write(table_name = str(_output.__class__.__name__), + containers = _output, + ) + elif isinstance(_output,TriggerMapContainer) : + for i,key in enumerate(_output.containers.keys()) : + self.writer.write(table_name = f"{_output.containers[key].__class__.__name__}_{i}/{key.name}", + containers = _output.containers[key], + ) + else : + raise TypeError("component output must be an instance of TriggerMapContainer or NectarCAMContainer") + + self.writer.close() + super().finish() + self.log.warning("Shutting down.") + + @property + def event_source(self): + """ + Getter method for the _event_source attribute. + """ + return copy.copy(self._event_source) + + @event_source.setter + def event_source(self, value): + """ + Setter method to set a new NectarCAMEventSource to the _reader attribute. + + Args: + value: a NectarCAMEventSource instance. + """ + if isinstance(value, NectarCAMEventSource): + self._event_source = value + else: + raise TypeError("The reader must be a NectarCAMEventSource") + + @property + def _npixels(self): + """ + Getter method for the _npixels attribute. + """ + return self.__npixels + + @property + def _pixels_id(self): + """ + Getter method for the _pixels_id attribute. + """ + return self.__pixels_id + + @property + def npixels(self): + """ + Getter method for the npixels attribute. + """ + return copy.deepcopy(self.__npixels) + + @property + def pixels_id(self): + """ + Getter method for the pixels_id attribute. + """ + return copy.deepcopy(self.__pixels_id) + + +def main(): + """run the tool""" + tool = EventsLoopNectarCAMCalibrationTool() + tool.run() + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/src/nectarchain/makers/extractor/__init__.py b/src/nectarchain/makers/extractor/__init__.py new file mode 100644 index 00000000..2c75a364 --- /dev/null +++ b/src/nectarchain/makers/extractor/__init__.py @@ -0,0 +1 @@ +# from .charge_extractor import * diff --git a/src/nectarchain/makers/extractor/charge_extractor.py b/src/nectarchain/makers/extractor/charge_extractor.py new file mode 100644 index 00000000..670234c5 --- /dev/null +++ b/src/nectarchain/makers/extractor/charge_extractor.py @@ -0,0 +1,128 @@ +import logging + +import numpy as np +from ctapipe.image import ImageExtractor +from numba import guvectorize +from scipy.interpolate import InterpolatedUnivariateSpline +from scipy.signal import find_peaks + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + +__all__ = ["gradient_extractor"] + + +class gradient_extractor(ImageExtractor): + fixed_window = np.uint8(16) + height_peak = np.uint8(10) + + def __call__(self, waveforms, telid, selected_gain_channel, substract_ped=False): + shape = waveforms.shape + waveforms = waveforms.reshape(shape[0] * shape[1], shape[2]) + + if substract_ped: + log.info("substracting pedestal") + # calculate pedestal + ped_mean = np.mean(waveforms[:, :16], axis=1).reshape( + shape[0] * shape[1], 1 + ) + + y = waveforms - (ped_mean) @ ( + np.ones(1, shape[2]) + ) # waveforms without pedestal + else: + log.info("do not substract pedestal") + y = waveforms + + waveforms.reshape(shape[0], shape[1], shape[2]) + + ( + peak_time, + charge_integral, + charge_sum, + charge_integral_fixed_window, + charge_sum_fixed_window, + ) = extract_charge(waveforms, self.height_peak, self.fixed_window) + + return peak_time, charge_integral + + +@guvectorize( + [ + "(np.uint16[:], np.uint8, np.uint8, np.uint16[:])", + ], + "(n,p,s),(),()->(n,p)", + nopython=True, + cache=True, +) +def extract_charge(y, height_peak, fixed_window): + x = np.linspace(0, len(y), len(y)) + xi = np.linspace(0, len(y), 251) + ius = InterpolatedUnivariateSpline(x, y) + yi = ius(xi) + peaks, _ = find_peaks(yi, height=height_peak) + # find the peak + if len(peaks) > 0: + # find the max peak + # max_peak = max(yi[peaks]) + max_peak_index = np.argmax(yi[peaks], axis=0) + # Check if there is not a peak but a plateaux + # 1. divide for the maximum to round to the first digit + # 2. create a new array with normalized and rounded values, yi_rounded + yi_rounded = np.around(yi[peaks] / max(yi[peaks]), 1) + maxima_peak_index = np.argwhere(yi_rounded == np.amax(yi_rounded)) + if len(maxima_peak_index) > 1: + # saturated event + max_peak_index = int(np.median(maxima_peak_index)) + # width_peak = 20 + if (xi[peaks[max_peak_index]] > 20) & ( + xi[peaks[max_peak_index]] < 40 + ): # Search the adaptive integration window + # calculate total gradients (not used, only for plot) + yi_grad_tot = np.gradient(yi, 1) + maxposition = peaks[max_peak_index] + # calcualte grandients starting from the max peak and going to the left to find the left margin of the window + yi_left = yi[:maxposition] + yi_grad_left = np.gradient(yi_left[::-1], 0.9) + change_grad_pos_left = ( + np.where(yi_grad_left[:-1] * yi_grad_left[1:] < 0)[0] + 1 + )[0] + # calcualte grandients starting from the max peak and going to the right to find the right margin of the window + yi_right = yi[maxposition:] + yi_grad_right = np.gradient(yi_right, 0.5) + change_grad_pos_right = ( + np.where(yi_grad_right[:-1] * yi_grad_right[1:] < 0)[0] + 1 + )[0] + charge_integral = ius.integral( + xi[peaks[max_peak_index] - (change_grad_pos_left)], + xi[peaks[max_peak_index] + change_grad_pos_right], + ) + charge_sum = yi[ + (peaks[max_peak_index] - (change_grad_pos_left)) : peaks[max_peak_index] + + change_grad_pos_right + ].sum() / ( + change_grad_pos_left + change_grad_pos_right + ) # simple sum integration + adaptive_window = change_grad_pos_right + change_grad_pos_left + window_right = (fixed_window - 6) / 2 + window_left = (fixed_window - 6) / 2 + 6 + charge_integral_fixed_window = ius.integral( + xi[peaks[max_peak_index] - (window_left)], + xi[peaks[max_peak_index] + window_right], + ) + charge_sum_fixed_window = yi[ + (peaks[max_peak_index] - (window_left)) : peaks[max_peak_index] + + window_right + ].sum() / ( + fixed_window + ) # simple sum integration + else: + log.info("No peak found, maybe it is a pedestal or noisy run!") + return ( + adaptive_window, + charge_integral, + charge_sum, + charge_integral_fixed_window, + charge_sum_fixed_window, + ) diff --git a/src/nectarchain/makers/extractor/tests/test_utils.py b/src/nectarchain/makers/extractor/tests/test_utils.py new file mode 100644 index 00000000..fc7dcf4d --- /dev/null +++ b/src/nectarchain/makers/extractor/tests/test_utils.py @@ -0,0 +1,22 @@ +import pytest + + +class TestCtapipeExtractor: + @pytest.mark.skip("numba conflict") + # Tests that the function returns the image and peak_time values from a valid DL1CameraContainer object. + def test_get_image_peak_time_valid_object(self): + from ctapipe.containers import DL1CameraContainer + + from nectarchain.makers.extractor.utils import CtapipeExtractor + + # Create a valid DL1CameraContainer object + container = DL1CameraContainer() + container.image = [1, 2, 3, 4, 5] + container.peak_time = [10, 4, 5, 6, 9] + + # Call the function under test + result_image, result_peak_time = CtapipeExtractor.get_image_peak_time(container) + + # Check the result + assert result_image == [1, 2, 3, 4, 5] + assert result_peak_time == [10, 4, 5, 6, 9] diff --git a/src/nectarchain/makers/extractor/utils.py b/src/nectarchain/makers/extractor/utils.py new file mode 100644 index 00000000..0fab1212 --- /dev/null +++ b/src/nectarchain/makers/extractor/utils.py @@ -0,0 +1,25 @@ +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.containers import DL1CameraContainer + + +class CtapipeExtractor: + """ + A class to extract the image and peak time from a DL1CameraContainer object. + """ + + def get_image_peak_time(cameraContainer): + """ + Extracts the image and peak time from a DL1CameraContainer object. + + Parameters: + cameraContainer (DL1CameraContainer): The DL1CameraContainer object to extract the image and peak time from. + + Returns: + tuple: A tuple containing the image and peak time values from the container. + """ + return cameraContainer.image, cameraContainer.peak_time diff --git a/src/nectarchain/makers/tests/test_chargesMakers.py b/src/nectarchain/makers/tests/test_chargesMakers.py new file mode 100644 index 00000000..04409263 --- /dev/null +++ b/src/nectarchain/makers/tests/test_chargesMakers.py @@ -0,0 +1,197 @@ +import logging + +import numpy as np +from ctapipe.containers import EventType + +from nectarchain.data.container import ChargesContainer, ChargesContainerIO +from nectarchain.makers import ChargesMaker, WaveformsMaker + +logging.basicConfig( + format="%(asctime)s %(name)s %(levelname)s %(message)s", level=logging.DEBUG +) +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + + +class TestChargesMaker: + run_number = 3938 + max_events = 100 + + def test_instance(self): + chargesMaker = ChargesMaker( + run_number=TestChargesMaker.run_number, + max_events=TestChargesMaker.max_events, + ) + assert isinstance(chargesMaker, ChargesMaker) + + def test_shape_valid(self): + chargesMaker = ChargesMaker( + run_number=TestChargesMaker.run_number, + max_events=TestChargesMaker.max_events, + ) + chargesContainer = chargesMaker.make()[0] + + assert chargesContainer.nevents <= TestChargesMaker.max_events + assert chargesContainer.run_number == TestChargesMaker.run_number + assert chargesContainer.ucts_timestamp.shape == (chargesContainer.nevents,) + assert chargesContainer.ucts_busy_counter.shape == (chargesContainer.nevents,) + assert chargesContainer.ucts_event_counter.shape == (chargesContainer.nevents,) + assert chargesContainer.event_type.shape == (chargesContainer.nevents,) + assert chargesContainer.event_id.shape == (chargesContainer.nevents,) + assert chargesContainer.trig_pattern_all.shape[0] == chargesContainer.nevents + assert chargesContainer.trig_pattern_all.shape[2] == 4 + + assert chargesContainer.trig_pattern.shape[0] == chargesContainer.nevents + assert chargesContainer.multiplicity.shape == (chargesContainer.nevents,) + + assert chargesContainer.charges_hg.mean() != 0 + assert chargesContainer.charges_lg.mean() != 0 + assert chargesContainer.peak_hg.mean() != 0 + assert chargesContainer.peak_lg.mean() != 0 + assert chargesContainer.charges_hg.shape == ( + chargesContainer.nevents, + chargesContainer.npixels, + ) + assert chargesContainer.charges_lg.shape == ( + chargesContainer.nevents, + chargesContainer.npixels, + ) + assert chargesContainer.peak_hg.shape == ( + chargesContainer.nevents, + chargesContainer.npixels, + ) + assert chargesContainer.peak_lg.shape == ( + chargesContainer.nevents, + chargesContainer.npixels, + ) + + assert chargesContainer.broken_pixels_hg.shape == ( + chargesContainer.nevents, + chargesContainer.npixels, + ) + assert chargesContainer.broken_pixels_lg.shape == ( + chargesContainer.nevents, + chargesContainer.npixels, + ) + + def test_make_restart_eventsource(self): + chargesMaker = ChargesMaker( + run_number=TestChargesMaker.run_number, + max_events=TestChargesMaker.max_events, + ) + chargesContainer_list = chargesMaker.make(restart_from_begining=True) + assert isinstance(chargesContainer_list[0], ChargesContainer) + + def test_make_LocalPeakWindowSum(self): + chargesMaker = ChargesMaker( + run_number=TestChargesMaker.run_number, + max_events=TestChargesMaker.max_events, + ) + chargesContainer_list = chargesMaker.make( + method="LocalPeakWindowSum", window_shift=-4, window_length=16 + ) + assert isinstance(chargesContainer_list[0], ChargesContainer) + + def test_all_multiple_trigger(self): + trigger1 = EventType.FLATFIELD + trigger2 = EventType.SKY_PEDESTAL + chargesMaker = ChargesMaker( + run_number=TestChargesMaker.run_number, + max_events=TestChargesMaker.max_events, + ) + chargesContainer_list = chargesMaker.make(trigger_type=[trigger1, trigger2]) + for chargesContainer in chargesContainer_list: + assert isinstance(chargesContainer, ChargesContainer) + + def test_all_trigger_None(self): + chargesMaker = ChargesMaker( + run_number=TestChargesMaker.run_number, + max_events=TestChargesMaker.max_events, + ) + chargesContainer_list = chargesMaker.make() + assert isinstance(chargesContainer_list[0], ChargesContainer) + + def test_create_from_waveforms(self): + waveformsMaker = WaveformsMaker( + run_number=TestChargesMaker.run_number, + max_events=TestChargesMaker.max_events, + ) + waveformsContainer_list = waveformsMaker.make() + chargesContainer = ChargesMaker.create_from_waveforms( + waveformsContainer_list[0], + method="LocalPeakWindowSum", + window_shift=-4, + window_length=16, + ) + assert isinstance(chargesContainer, ChargesContainer) + + def test_select_charges(self): + chargesMaker = ChargesMaker( + run_number=TestChargesMaker.run_number, + max_events=TestChargesMaker.max_events, + ) + chargesContainer_list = chargesMaker.make() + pixel_id = np.array([3, 67, 87]) + assert isinstance( + ChargesMaker.select_charges_hg(chargesContainer_list[0], pixel_id), + np.ndarray, + ) + assert isinstance( + ChargesMaker.select_charges_lg(chargesContainer_list[0], pixel_id), + np.ndarray, + ) + + def test_histo(self): + chargesMaker = ChargesMaker( + run_number=TestChargesMaker.run_number, + max_events=TestChargesMaker.max_events, + ) + chargesContainer_list = chargesMaker.make() + histo = ChargesMaker.histo_hg(chargesContainer_list[0]) + assert isinstance(histo, np.ndarray) + assert histo.mean() != 0 + assert histo.shape[0] == 2 + assert histo.shape[1] == chargesContainer_list[0].npixels + histo = ChargesMaker.histo_lg(chargesContainer_list[0]) + assert isinstance(histo, np.ndarray) + assert histo.mean() != 0 + assert histo.shape[0] == 2 + assert histo.shape[1] == chargesContainer_list[0].npixels + + def test_sort_ChargesContainer(self): + chargesMaker = ChargesMaker( + run_number=TestChargesMaker.run_number, + max_events=TestChargesMaker.max_events, + ) + chargesContainer_list = chargesMaker.make() + sortWfs = ChargesMaker.sort(chargesContainer_list[0], method="event_id") + assert np.array_equal( + sortWfs.event_id, np.sort(chargesContainer_list[0].event_id) + ) + + def test_write_load_container(self): + chargesMaker = ChargesMaker( + run_number=TestChargesMaker.run_number, + max_events=TestChargesMaker.max_events, + ) + chargesContainer_list = chargesMaker.make() + ChargesContainerIO.write( + "/tmp/test_charge_container/", chargesContainer_list[0], overwrite=True + ) + loaded_charge = ChargesContainerIO.load( + f"/tmp/test_charge_container", run_number=TestChargesMaker.run_number + ) + assert np.array_equal( + chargesContainer_list[0].charges_hg, loaded_charge.charges_hg + ) + + +if __name__ == "__main__": + import logging + + logging.basicConfig( + format="%(asctime)s %(name)s %(levelname)s %(message)s", level=logging.DEBUG + ) + log = logging.getLogger(__name__) + log.handlers = logging.getLogger("__main__").handlers + TestChargesMaker().test_write_load_container() diff --git a/src/nectarchain/makers/tests/test_core.py b/src/nectarchain/makers/tests/test_core.py new file mode 100644 index 00000000..98ae818f --- /dev/null +++ b/src/nectarchain/makers/tests/test_core.py @@ -0,0 +1,43 @@ +import logging + +from nectarchain.data.container import ChargesContainer +from nectarchain.makers import ArrayDataMaker, ChargesMaker, WaveformsMaker + +logging.basicConfig( + format="%(asctime)s %(name)s %(levelname)s %(message)s", level=logging.DEBUG +) +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + + +class TestArrayDataMaker: + run_number = 3938 + max_events = 100 + + def test_merge(self): + chargesMaker = ChargesMaker( + run_number=TestArrayDataMaker.run_number, + max_events=TestArrayDataMaker.max_events, + ) + charges_1 = chargesMaker.make() + chargesMaker_2 = ChargesMaker( + run_number=TestArrayDataMaker.run_number, + max_events=TestArrayDataMaker.max_events, + ) + charges_2 = chargesMaker_2.make() + + merged = ArrayDataMaker.merge(charges_1, charges_2) + assert isinstance(merged, ChargesContainer) + + def test_merge_different_container(self): + chargesMaker = ChargesMaker( + run_number=TestArrayDataMaker.run_number, + max_events=TestArrayDataMaker.max_events, + ) + charges_1 = chargesMaker.make() + wfsMaker_2 = WaveformsMaker( + run_number=TestArrayDataMaker.run_number, + max_events=TestArrayDataMaker.max_events, + ) + wfs_2 = wfsMaker_2.make() + merged = ArrayDataMaker.merge(charges_1, wfs_2) diff --git a/src/nectarchain/makers/tests/test_waveformsMakers.py b/src/nectarchain/makers/tests/test_waveformsMakers.py new file mode 100644 index 00000000..35aaec16 --- /dev/null +++ b/src/nectarchain/makers/tests/test_waveformsMakers.py @@ -0,0 +1,167 @@ +import logging + +import numpy as np +from ctapipe.containers import EventType + +from nectarchain.data.container import ( + ArrayDataContainer, + WaveformsContainer, + WaveformsContainerIO, +) +from nectarchain.makers.waveformsMakers import WaveformsMaker + +logging.basicConfig( + format="%(asctime)s %(name)s %(levelname)s %(message)s", level=logging.DEBUG +) +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + + +class TestWaveformsMaker: + run_number = 3938 + max_events = 100 + + def test_instance(self): + waveformsMaker = WaveformsMaker( + run_number=TestWaveformsMaker.run_number, + max_events=TestWaveformsMaker.max_events, + ) + assert isinstance(waveformsMaker, WaveformsMaker) + + def test_shape_valid(self): + waveformsMaker = WaveformsMaker( + run_number=TestWaveformsMaker.run_number, + max_events=TestWaveformsMaker.max_events, + ) + waveformsContainer = waveformsMaker.make()[0] + + assert waveformsContainer.nevents <= TestWaveformsMaker.max_events + assert waveformsContainer.run_number == TestWaveformsMaker.run_number + assert waveformsContainer.ucts_timestamp.shape == (waveformsContainer.nevents,) + assert waveformsContainer.ucts_busy_counter.shape == ( + waveformsContainer.nevents, + ) + assert waveformsContainer.ucts_event_counter.shape == ( + waveformsContainer.nevents, + ) + assert waveformsContainer.event_type.shape == (waveformsContainer.nevents,) + assert waveformsContainer.event_id.shape == (waveformsContainer.nevents,) + assert ( + waveformsContainer.trig_pattern_all.shape[0] == waveformsContainer.nevents + ) + assert waveformsContainer.trig_pattern_all.shape[2] == 4 + + assert waveformsContainer.trig_pattern.shape[0] == waveformsContainer.nevents + assert waveformsContainer.multiplicity.shape == (waveformsContainer.nevents,) + + assert waveformsContainer.wfs_hg.mean() != 0 + assert waveformsContainer.wfs_lg.mean() != 0 + assert waveformsContainer.wfs_hg.shape == ( + waveformsContainer.nevents, + waveformsContainer.npixels, + waveformsContainer.nsamples, + ) + assert waveformsContainer.wfs_lg.shape == ( + waveformsContainer.nevents, + waveformsContainer.npixels, + waveformsContainer.nsamples, + ) + assert waveformsContainer.broken_pixels_hg.shape == ( + waveformsContainer.nevents, + waveformsContainer.npixels, + ) + assert waveformsContainer.broken_pixels_lg.shape == ( + waveformsContainer.nevents, + waveformsContainer.npixels, + ) + + def test_all_multiple_trigger(self): + trigger1 = EventType.FLATFIELD + trigger2 = EventType.SKY_PEDESTAL + waveformsMaker = WaveformsMaker( + run_number=TestWaveformsMaker.run_number, + max_events=TestWaveformsMaker.max_events, + ) + waveformsContainer_list = waveformsMaker.make( + trigger_type=[trigger1, trigger2], restart_from_begining=True + ) + for waveformsContainer in waveformsContainer_list: + assert isinstance(waveformsContainer, WaveformsContainer) + assert waveformsContainer.wfs_hg.mean() != 0 + + def test_all_trigger_None(self): + waveformsMaker = WaveformsMaker( + run_number=TestWaveformsMaker.run_number, + max_events=TestWaveformsMaker.max_events, + ) + waveformsContainer_list = waveformsMaker.make() + assert isinstance(waveformsContainer_list[0], WaveformsContainer) + + def test_select_waveforms_hg(self): + waveformsMaker = WaveformsMaker( + run_number=TestWaveformsMaker.run_number, + max_events=TestWaveformsMaker.max_events, + ) + waveformsContainer_list = waveformsMaker.make() + pixel_id = np.array([3, 67, 87]) + assert isinstance( + WaveformsMaker.select_waveforms_hg(waveformsContainer_list[0], pixel_id), + np.ndarray, + ) + assert isinstance( + WaveformsMaker.select_waveforms_lg(waveformsContainer_list[0], pixel_id), + np.ndarray, + ) + + def test_sort_WaveformsContainer(self): + waveformsMaker = WaveformsMaker( + run_number=TestWaveformsMaker.run_number, + max_events=TestWaveformsMaker.max_events, + ) + waveformsContainer_list = waveformsMaker.make() + sortWfs = WaveformsMaker.sort(waveformsContainer_list[0], method="event_id") + assert np.array_equal( + sortWfs.event_id, np.sort(waveformsContainer_list[0].event_id) + ) + + def test_write_load_container(self): + waveformsMaker = WaveformsMaker( + run_number=TestWaveformsMaker.run_number, + max_events=TestWaveformsMaker.max_events, + ) + waveformsContainer_list = waveformsMaker.make() + WaveformsContainerIO.write( + "/tmp/test_wfs_container/", waveformsContainer_list[0], overwrite=True + ) + loaded_wfs = WaveformsContainerIO.load( + f"/tmp/test_wfs_container", TestWaveformsMaker.run_number + ) + assert np.array_equal(waveformsContainer_list[0].wfs_hg, loaded_wfs.wfs_hg) + + def test_create_from_events_list(self): + waveformsMaker = WaveformsMaker( + run_number=TestWaveformsMaker.run_number, + max_events=TestWaveformsMaker.max_events, + ) + events_list = [] + for i, event in enumerate(waveformsMaker._reader): + events_list.append(event) + waveformsContainer = WaveformsMaker.create_from_events_list( + events_list, + waveformsMaker.run_number, + waveformsMaker.npixels, + waveformsMaker.nsamples, + waveformsMaker.subarray, + waveformsMaker.pixels_id, + ) + assert isinstance(waveformsContainer, WaveformsContainer) + + +if __name__ == "__main__": + import logging + + logging.basicConfig( + format="%(asctime)s %(name)s %(levelname)s %(message)s", level=logging.DEBUG + ) + log = logging.getLogger(__name__) + log.handlers = logging.getLogger("__main__").handlers diff --git a/src/nectarchain/makers/waveformsMakers.py b/src/nectarchain/makers/waveformsMakers.py new file mode 100644 index 00000000..04105fb8 --- /dev/null +++ b/src/nectarchain/makers/waveformsMakers.py @@ -0,0 +1,32 @@ +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 +from ctapipe.containers import EventType +from ctapipe.instrument import SubarrayDescription +from ctapipe_io_nectarcam import constants +from ctapipe_io_nectarcam.containers import NectarCAMDataContainer +from ctapipe.core.traits import ComponentNameList +from tqdm import tqdm + +from ..data.container import WaveformsContainer +from .core import EventsLoopNectarCAMCalibrationTool +from .component import NectarCAMComponent + +__all__ = ["WaveformsNectarCAMCalibrationTool"] + + +class WaveformsNectarCAMCalibrationTool(EventsLoopNectarCAMCalibrationTool): + """class use to make the waveform extraction from event read from r0 data""" + componentsList = ComponentNameList( + NectarCAMComponent, + default_value = ["WaveformsComponent"], + help="List of Component names to be apply, the order will be respected" + ).tag(config=True) + \ No newline at end of file diff --git a/src/nectarchain/tests/test_version.py b/src/nectarchain/tests/test_version.py index b62c0c66..c0c1392d 100644 --- a/src/nectarchain/tests/test_version.py +++ b/src/nectarchain/tests/test_version.py @@ -1,4 +1,4 @@ def test_version(): from nectarchain import __version__ - assert __version__ != '0.0.0' \ No newline at end of file + assert __version__ != "0.0.0" diff --git a/src/nectarchain/user_scripts/ggrolleron/gain_PhotoStat_computation.py b/src/nectarchain/user_scripts/ggrolleron/gain_PhotoStat_computation.py index 3ad357ab..44bd47ec 100644 --- a/src/nectarchain/user_scripts/ggrolleron/gain_PhotoStat_computation.py +++ b/src/nectarchain/user_scripts/ggrolleron/gain_PhotoStat_computation.py @@ -1,113 +1,148 @@ import logging -import sys import os +import sys from pathlib import Path + import numpy as np -os.makedirs(os.environ.get('NECTARCHAIN_LOG'),exist_ok = True) +os.makedirs(os.environ.get("NECTARCHAIN_LOG"), exist_ok=True) -#to quiet numba -logging.basicConfig(format='%(asctime)s %(name)s %(levelname)s %(message)s') +# to quiet numba +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") logging.getLogger("numba").setLevel(logging.WARNING) import argparse -from nectarchain.calibration.NectarGain import PhotoStatGainFFandPed +from astropy.table import QTable + +from nectarchain.makers.calibration.gain.PhotoStatisticMakers import PhotoStatisticMaker 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}/' - ) - -#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 - ) - -#tag for SPE fit results propagation -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='' - ) - -parser.add_argument('--overwrite', - action='store_true', - default=False, - help='to force overwrite files on disk' - ) - -#for plotting correlation -parser.add_argument('--correlation', - action='store_true', - default=True, - help='to plot correlation between SPE gain computation and Photo-statistic gain resluts' - ) - -#extractor arguments -parser.add_argument('--chargeExtractorPath', - help='charge extractor path where charges are saved', - type=str - ) - -parser.add_argument('--FFchargeExtractorWindowLength', - help='charge extractor window length in ns', - 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) - - -def main(args) : - figpath = os.environ.get('NECTARCHAIN_FIGURES') - - - photoStat_FFandPed = PhotoStatGainFFandPed(args.FF_run_number, args.ped_run_number, SPEresults = args.SPE_fit_results,method = args.chargeExtractorPath, FFchargeExtractorWindowLength = args.FFchargeExtractorWindowLength) - photoStat_FFandPed.run() - 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) + 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}/", +) + +# 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, +) + +# tag for SPE fit results propagation +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="", +) + +parser.add_argument( + "--overwrite", + action="store_true", + default=False, + help="to force overwrite files on disk", +) + +# for plotting correlation +parser.add_argument( + "--correlation", + action="store_true", + default=True, + help="to plot correlation between SPE gain computation and Photo-statistic gain resluts", +) + +# extractor arguments +parser.add_argument( + "--chargeExtractorPath", + help="charge extractor path where charges are saved", + type=str, +) + +parser.add_argument( + "--FFchargeExtractorWindowLength", + help="charge extractor window length in ns", + 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, +) + + +def main(args): + figpath = os.environ.get("NECTARCHAIN_FIGURES") + + 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, + ) + 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, + ) 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 : - fig = photoStat_FFandPed.plot_correlation() - 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 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 __name__ == "__main__": args = parser.parse_args() logginglevel = logging.FATAL - if args.verbosity == "warning" : + if args.verbosity == "warning": logginglevel = logging.WARNING - elif args.verbosity == "info" : + elif args.verbosity == "info": logginglevel = logging.INFO - elif args.verbosity == "debug" : + elif args.verbosity == "debug": logginglevel = logging.DEBUG os.makedirs(f"{os.environ.get('NECTARCHAIN_LOG')}/{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") + 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", + ) log = logging.getLogger(__name__) log.setLevel(logginglevel) ##tips to add message to stdout handler = logging.StreamHandler(sys.stdout) handler.setLevel(logginglevel) - formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s') + formatter = logging.Formatter( + "%(asctime)s - %(name)s - %(levelname)s - %(message)s" + ) handler.setFormatter(formatter) log.addHandler(handler) - main(args) \ No newline at end of file + main(args) diff --git a/src/nectarchain/user_scripts/ggrolleron/gain_PhotoStat_computation.sh b/src/nectarchain/user_scripts/ggrolleron/gain_PhotoStat_computation.sh index aa4c8dc8..e6809ed3 100644 --- a/src/nectarchain/user_scripts/ggrolleron/gain_PhotoStat_computation.sh +++ b/src/nectarchain/user_scripts/ggrolleron/gain_PhotoStat_computation.sh @@ -2,4 +2,4 @@ #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 -p 3938 -f 3937 --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 -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 c349e0b2..80b90c47 100644 --- a/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_combined_computation.py +++ b/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_combined_computation.py @@ -1,143 +1,176 @@ import logging -import sys import os -from pathlib import Path +import sys import time +from pathlib import Path -os.makedirs(os.environ.get('NECTARCHAIN_LOG'),exist_ok = True) +os.makedirs(os.environ.get("NECTARCHAIN_LOG"), exist_ok=True) -#to quiet numba +# to quiet numba logging.getLogger("numba").setLevel(logging.WARNING) import argparse -#import seaborn as sns -from nectarchain.calibration.container import ChargeContainer -from nectarchain.calibration.NectarGain import NectarGainSPESingleSignalfromHHVFit +# import seaborn as sns +from nectarchain.data.container import ChargeContainer +from nectarchain.makers.calibration.gain.FlatFieldSPEMakers import ( + FlatFieldSingleNominalSPEMaker, +) 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}/' - ) - -#run numbers -parser.add_argument('-r', '--run_number', - help='spe run', - type=int) - -parser.add_argument('--overwrite', - action='store_true', - default=False, - help='to force overwrite files on disk' - ) - -#output figures path extension -parser.add_argument('--output_fig_tag', - type = str, - default='', - help='tag to set output figure path' - ) - -#pixels selected -parser.add_argument('-p','--pixels', - 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', - 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 arguments -parser.add_argument('--chargeExtractorPath', - help='charge extractor path where charges are saved', - type=str - ) - -#for VVH combined fit -parser.add_argument('--combined', - 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' - ) -parser.add_argument('--VVH_fitted_results', - help='previoulsy fitted VVH data path for nominal SPE fit by fixing some shared parameters', - type=str - ) -#tag for SPE fit results propagation -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='' - ) -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) - - -def main(args) : - figpath = f"{os.environ.get('NECTARCHAIN_FIGURES')}/" + 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}/", +) + +# run numbers +parser.add_argument("-r", "--run_number", help="spe run", type=int) + +parser.add_argument( + "--overwrite", + action="store_true", + default=False, + help="to force overwrite files on disk", +) + +# output figures and path extension +parser.add_argument( + "--display", action="store_true", default=False, help="whether to save plot or not" +) +parser.add_argument( + "--output_fig_tag", type=str, default="", help="tag to set output figure path" +) + +# pixels selected +parser.add_argument( + "-p", "--pixels", 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", 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 arguments +parser.add_argument( + "--chargeExtractorPath", + help="charge extractor path where charges are saved", + type=str, +) + +# for VVH combined fit +parser.add_argument( + "--combined", + 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", +) +parser.add_argument( + "--VVH_fitted_results", + help="previoulsy fitted VVH data path for nominal SPE fit by fixing some shared parameters", + type=str, +) +# tag for SPE fit results propagation +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="", +) +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, +) + + +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 "" - charge_run = ChargeContainer.from_file(f"{os.environ.get('NECTARCAMDATA')}/charges/{args.chargeExtractorPath}/",args.run_number) + charge_run = ChargeContainer.from_file( + f"{os.environ.get('NECTARCAMDATA')}/charges/{args.chargeExtractorPath}/", + args.run_number, + ) - if args.combined : + if args.combined: raise NotImplementedError("combined fit not implemented yet") - else : - gain_Std = NectarGainSPESingleSignalfromHHVFit(signal = charge_run, - nectarGainSPEresult=args.VVH_fitted_results, - same_luminosity=args.same_luminosity - ) + 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.run(pixel = args.pixels, multiproc = args.multiproc, 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}") + 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}", + ) + 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) - log.info(f"convergence rate : {len(gain_Std._output_table[gain_Std._output_table['is_valid']])/gain_Std.npixels}") + 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, + ) + log.info( + f"convergence rate : {len(gain_Std._results[gain_Std._results['is_valid']])/gain_Std.npixels}" + ) + if __name__ == "__main__": args = parser.parse_args() logginglevel = logging.FATAL - if args.verbosity == "warning" : + if args.verbosity == "warning": logginglevel = logging.WARNING - elif args.verbosity == "info" : + elif args.verbosity == "info": logginglevel = logging.INFO - elif args.verbosity == "debug" : + elif args.verbosity == "debug": logginglevel = logging.DEBUG - os.makedirs(f"{os.environ.get('NECTARCHAIN_LOG')}/{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") - + os.makedirs( + f"{os.environ.get('NECTARCHAIN_LOG','/tmp/nectarchain_log')}/{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", + ) + log = logging.getLogger(__name__) log.setLevel(logginglevel) ##tips to add message to stdout handler = logging.StreamHandler(sys.stdout) handler.setLevel(logginglevel) - formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s') + formatter = logging.Formatter( + "%(asctime)s - %(name)s - %(levelname)s - %(message)s" + ) handler.setFormatter(formatter) log.addHandler(handler) - main(args) \ No newline at end of file + main(args) 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 7a4dd733..609a4138 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,9 @@ #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 --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" +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 --SPE_fit_results_tag 2634 --chargeExtractorPath LocalPeakWindowSum_4-12 --overwrite --multiproc --nproc 6 --chunksize 50 --same_luminosity --combined +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 2634 (taken at 1400V) (luminosity can not be supposed the same) -python gain_SPEfit_combined_computation.py -r 3936 --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-2634-LocalPeakWindowSum_4-12/output_table.ecsv" +#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" diff --git a/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_computation.py b/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_computation.py index 2fc0177e..514a4758 100644 --- a/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_computation.py +++ b/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_computation.py @@ -1,132 +1,162 @@ import logging -import sys import os -from pathlib import Path +import sys import time +from pathlib import Path -os.makedirs(os.environ.get('NECTARCHAIN_LOG'),exist_ok = True) +os.makedirs(os.environ.get("NECTARCHAIN_LOG"), exist_ok=True) -#to quiet numba +# to quiet numba logging.getLogger("numba").setLevel(logging.WARNING) import argparse -#import seaborn as sns -from nectarchain.calibration.container import ChargeContainer -from nectarchain.calibration.NectarGain import NectarGainSPESingleSignalStd,NectarGainSPESingleSignal +# import seaborn as sns +from nectarchain.data.container import ChargeContainer +from nectarchain.makers.calibration.gain.FlatFieldSPEMakers import ( + FlatFieldSingleHHVSPEMaker, + FlatFieldSingleHHVStdSPEMaker, +) 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}/' - ) - -#run numbers -parser.add_argument('-r', '--run_number', - help='spe run', - type=int) - -parser.add_argument('--overwrite', - action='store_true', - default=False, - help='to force overwrite files on disk' - ) - -parser.add_argument('--voltage_tag', - type = str, - default='', - help='tag for voltage specifcication (1400V or nominal), used to setup the output path. See help for more details' - ) - -#output figures path extension -parser.add_argument('--output_fig_tag', - type = str, - default='', - help='tag to set output figure path' - ) - -#pixels selected -parser.add_argument('-p','--pixels', - nargs="+", - default=None, - help='pixels selected', - type=int) - -#for let free pp and n : -parser.add_argument('--free_pp_n', - action='store_true', - default=False, - help='to let free pp and n' - ) - -#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 - ) - -#verbosity argument -parser.add_argument('-v',"--verbosity", - help='set the verbosity level of logger', - default="info", - choices=["fatal","debug","info","warning"], - type=str) - - -def main(args) : + 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}/", +) + +# run numbers +parser.add_argument("-r", "--run_number", help="spe run", type=int) + +parser.add_argument( + "--overwrite", + action="store_true", + default=False, + help="to force overwrite files on disk", +) + +parser.add_argument( + "--voltage_tag", + type=str, + default="", + help="tag for voltage specifcication (1400V or nominal), used to setup the output path. See help for more details", +) + +# output figures and path extension +parser.add_argument( + "--display", action="store_true", default=False, help="whether to save plot or not" +) +parser.add_argument( + "--output_fig_tag", type=str, default="", help="tag to set output figure path" +) + +# pixels selected +parser.add_argument( + "-p", "--pixels", nargs="+", default=None, help="pixels selected", type=int +) + +# for let free pp and n : +parser.add_argument( + "--free_pp_n", action="store_true", default=False, help="to let free pp and n" +) + +# 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, +) + +# verbosity argument +parser.add_argument( + "-v", + "--verbosity", + help="set the verbosity level of logger", + default="info", + choices=["fatal", "debug", "info", "warning"], + type=str, +) + + +def main(args): figpath = f"{os.environ.get('NECTARCHAIN_FIGURES')}/" figpath_ext = "" if args.output_fig_tag == "" else f"-{args.output_fig_tag}" multipath = "MULTI-" if args.multiproc else "" SPEpath = "SPE" if args.free_pp_n else "SPEStd" - charge_run_1400V = ChargeContainer.from_file(f"{os.environ.get('NECTARCAMDATA')}/charges/{args.chargeExtractorPath}/",args.run_number) + charge_run_1400V = ChargeContainer.from_file( + f"{os.environ.get('NECTARCAMDATA')}/charges/{args.chargeExtractorPath}/", + args.run_number, + ) - if args.free_pp_n : - gain_Std = NectarGainSPESingleSignal(signal = charge_run_1400V) + if args.free_pp_n: + gain_Std = FlatFieldSingleHHVSPEMaker.create_from_chargeContainer( + signal=charge_run_1400V + ) + + else: + gain_Std = FlatFieldSingleHHVStdSPEMaker.create_from_chargeContainer( + signal=charge_run_1400V + ) - else : - gain_Std = NectarGainSPESingleSignalStd(signal = charge_run_1400V) t = time.time() - gain_Std.run(pixel = args.pixels, multiproc = args.multiproc, nproc = args.nproc, chunksize = args.chunksize, figpath = figpath+f"/{multipath}{args.voltage_tag}-{SPEpath}-{args.run_number}-{args.chargeExtractorPath}{figpath_ext}") + 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._output_table[gain_Std._output_table['is_valid']])/gain_Std.npixels if args.pixels is None else len(gain_Std._output_table[gain_Std._output_table['is_valid']])/len(args.pixels) + 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}") + if __name__ == "__main__": args = parser.parse_args() logginglevel = logging.FATAL - if args.verbosity == "warning" : + if args.verbosity == "warning": logginglevel = logging.WARNING - elif args.verbosity == "info" : + elif args.verbosity == "info": logginglevel = logging.INFO - elif args.verbosity == "debug" : + elif args.verbosity == "debug": logginglevel = logging.DEBUG os.makedirs(f"{os.environ.get('NECTARCHAIN_LOG')}/{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") + 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", + ) log = logging.getLogger(__name__) log.setLevel(logginglevel) ##tips to add message to stdout handler = logging.StreamHandler(sys.stdout) handler.setLevel(logginglevel) - formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s') + formatter = logging.Formatter( + "%(asctime)s - %(name)s - %(levelname)s - %(message)s" + ) handler.setFormatter(formatter) log.addHandler(handler) - main(args) \ No newline at end of file + main(args) diff --git a/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_computation.sh b/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_computation.sh index c290c10c..932d6721 100644 --- a/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_computation.sh +++ b/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_computation.sh @@ -1,10 +1,13 @@ #HOW TO : #to perform SPE fit of a HHV run -python gain_SPEfit_computation.py -r 2634 --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 +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 --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 +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 #to perform SPE fit of a run at nominal voltage -python gain_SPEfit_computation.py -r 3936 --voltage_tag "nominal" --chargeExtractorPath LocalPeakWindowSum_4-12 --overwrite --multiproc --nproc 6 --chunksize 50 +python gain_SPEfit_computation.py -r 3936 --display --voltage_tag "nominal" --chargeExtractorPath LocalPeakWindowSum_4-12 --overwrite --multiproc --nproc 6 --chunksize 50 + +#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 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 10112ce5..24249618 100644 --- a/src/nectarchain/user_scripts/ggrolleron/load_wfs_compute_charge.py +++ b/src/nectarchain/user_scripts/ggrolleron/load_wfs_compute_charge.py @@ -1,134 +1,160 @@ -from pathlib import Path -import sys -import os import argparse -import json import glob - +import json import logging +import os +import sys +from pathlib import Path + 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") +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", +) log = logging.getLogger(__name__) -from nectarchain.calibration.container import WaveformsContainer,WaveformsContainers,ChargeContainer,ChargeContainers +from nectarchain.data.container import ( + ChargeContainer, + ChargeContainers, + WaveformsContainer, + WaveformsContainers, +) parser = argparse.ArgumentParser( - prog = 'load_wfs_compute_charge', - description = 'This program load waveforms from fits.fz run files and compute charge') - -#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) - -#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', - 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', - type=int) - -#boolean arguments -parser.add_argument('--reload_wfs', - action='store_true', - default=False, - help='to force re-computation of waveforms from fits.fz files' - ) -parser.add_argument('--overwrite', - action='store_true', - default=False, - help='to force overwrite files on disk' - ) - -parser.add_argument('--split', - action='store_true', - default=False, - help='split waveforms extraction with 1 file per fits.fz raw data file' - ) - -#extractor arguments -parser.add_argument('--extractorMethod', - choices=["FullWaveformSum", - "FixedWindowSum", - "GlobalPeakWindowSum", - "LocalPeakWindowSum", - "SlidingWindowMaxSum", - "TwoPassWindowSum"], - default="LocalPeakWindowSum", - help='charge extractor method', - type=str - ) -parser.add_argument('--extractor_kwargs', - default={'window_width' : 16, 'window_shift' : 4}, - help='charge extractor kwargs', - type=json.loads - ) - -#verbosity argument -parser.add_argument('-v',"--verbosity", - help='set the verbosity level of logger', - default="info", - choices=["fatal","debug","info","warning"], - type=str) + prog="load_wfs_compute_charge", + description="This program load waveforms from fits.fz run files and compute charge", +) + +# 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 +) + +# 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", + 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", + type=int, +) + +# boolean arguments +parser.add_argument( + "--reload_wfs", + action="store_true", + default=False, + help="to force re-computation of waveforms from fits.fz files", +) +parser.add_argument( + "--overwrite", + action="store_true", + default=False, + help="to force overwrite files on disk", +) + +parser.add_argument( + "--split", + action="store_true", + default=False, + help="split waveforms extraction with 1 file per fits.fz raw data file", +) + +# extractor arguments +parser.add_argument( + "--extractorMethod", + choices=[ + "FullWaveformSum", + "FixedWindowSum", + "GlobalPeakWindowSum", + "LocalPeakWindowSum", + "SlidingWindowMaxSum", + "TwoPassWindowSum", + ], + default="LocalPeakWindowSum", + help="charge extractor method", + type=str, +) +parser.add_argument( + "--extractor_kwargs", + default={"window_width": 16, "window_shift": 4}, + help="charge extractor kwargs", + type=json.loads, +) + +# verbosity argument +parser.add_argument( + "-v", + "--verbosity", + help="set the verbosity level of logger", + default="info", + choices=["fatal", "debug", "info", "warning"], + type=str, +) args = parser.parse_args() -#control shape of arguments lists -for arg in ['spe','ff','ped'] : +# 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) + 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) + 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) : + +def load_wfs_no_split(i, runs_list, max_events, nevents, overwrite): """method to load waveforms without splitting Args: @@ -136,18 +162,21 @@ def load_wfs_no_split(i,runs_list,max_events,nevents,overwrite) : runs_list (list): the run number list max_events (list): max_events list nevents (list): nevents list - overwrite (bool): to overwrite + 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 = WaveformsContainer(runs_list[i], max_events=max_events[i], nevents=nevents[i]) wfs.load_wfs() - wfs.write(f"{os.environ['NECTARCAMDATA']}/waveforms/",overwrite = overwrite) + 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) : + +def load_wfs_charge_split( + i, runs_list, max_events, overwrite, charge_childpath, extractor_kwargs +): """_summary_ Args: @@ -155,37 +184,41 @@ def load_wfs_charge_split(i,runs_list,max_events,overwrite,charge_childpath,extr runs_list (list): the run number list max_events (list): max_events list nevents (list): nevents list - overwrite (bool): to overwrite + 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) + 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) : + 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) + 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)) + 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 + return wfs, charge -def load_wfs_charge_split_from_wfsFiles(wfsFiles,charge_childpath,extractor_kwargs) : + +def load_wfs_charge_split_from_wfsFiles(wfsFiles, charge_childpath, extractor_kwargs): """_summary_ Args: @@ -197,13 +230,15 @@ def load_wfs_charge_split_from_wfsFiles(wfsFiles,charge_childpath,extractor_kwar 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): + 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)) + 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 @@ -214,17 +249,20 @@ def load_wfs_charge_split_from_wfsFiles(wfsFiles,charge_childpath,extractor_kwar del wfs.event_id del wfs.trig_pattern_all del wfs - #gc.collect() - + # gc.collect() + log.info("merging charge") charge = charge.merge() - return None,charge + return None, charge -def load_wfs_compute_charge(runs_list : list, - reload_wfs : bool = False, - overwrite : bool= False, - charge_extraction_method : str = "FullWaveformSum", - **kwargs) -> None : + +def load_wfs_compute_charge( + runs_list: list, + reload_wfs: bool = False, + overwrite: bool = False, + charge_extraction_method: str = "FullWaveformSum", + **kwargs, +) -> None: """this method is used to load waveforms from zfits files and compute charge with an user specified method Args: @@ -237,119 +275,159 @@ def load_wfs_compute_charge(runs_list : list, 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))]) + # 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",{}) + charge_childpath = kwargs.get("charge_childpath", charge_extraction_method) + extractor_kwargs = kwargs.get("extractor_kwargs", {}) - split = kwargs.get("split",False) - + split = kwargs.get("split", False) - for i in range(len(runs_list)) : + 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 : + 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") + 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" + ) 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) + 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) - else : - wfs = load_wfs_no_split(i,runs_list,max_events,nevents,overwrite) - - - if not(split) : + else: + 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) + + if not (split): log.info(f"computation of charge with {charge_childpath}") - charge = ChargeContainer.from_waveforms(wfs,method = charge_childpath,**extractor_kwargs) + charge = ChargeContainer.from_waveforms( + wfs, method=charge_childpath, **extractor_kwargs + ) del wfs - charge.write(f"{os.environ['NECTARCAMDATA']}/charges/{path}/",overwrite = overwrite) + 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))]) +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") - - + charge_extraction_method = kwargs.get("extractorMethod", "FullWaveformSum") - load_wfs_compute_charge(runs_list = runs_list, - charge_extraction_method = charge_extraction_method, - nevents = nevents, - max_events = max_events, - **kwargs) + load_wfs_compute_charge( + runs_list=runs_list, + charge_extraction_method=charge_extraction_method, + nevents=nevents, + max_events=max_events, + **kwargs, + ) -if __name__ == '__main__': - - #run of interest - #spe_run_number = [2633,2634,3784] - #ff_run_number = [2608] - #ped_run_number = [2609] - #spe_nevents = [49227,49148,-1] +if __name__ == "__main__": + # run of interest + # spe_run_number = [2633,2634,3784] + # ff_run_number = [2608] + # ped_run_number = [2609] + # spe_nevents = [49227,49148,-1] args = parser.parse_args() logginglevel = logging.FATAL - if args.verbosity == "warning" : + if args.verbosity == "warning": logginglevel = logging.WARNING - elif args.verbosity == "info" : + elif args.verbosity == "info": logginglevel = logging.INFO - elif args.verbosity == "debug" : + elif args.verbosity == "debug": logginglevel = logging.DEBUG os.makedirs(f"{os.environ.get('NECTARCHAIN_LOG')}/{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") + 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", + ) log = logging.getLogger(__name__) log.setLevel(logginglevel) ##tips to add message to stdout handler = logging.StreamHandler(sys.stdout) handler.setLevel(logginglevel) - formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s') + formatter = logging.Formatter( + "%(asctime)s - %(name)s - %(levelname)s - %(message)s" + ) handler.setFormatter(formatter) log.addHandler(handler) @@ -357,23 +435,23 @@ def main(spe_run_number : list = [], log.info(f"arguments are : {arg}") key_to_pop = [] - for key in arg.keys() : - if arg[key] is None : + for key in arg.keys(): + if arg[key] is None: key_to_pop.append(key) - for key in key_to_pop : + 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']}" - - arg['path'] = path - + + 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']}" + + arg["path"] = path + main(**arg) 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 new file mode 100644 index 00000000..94b7fa53 --- /dev/null +++ b/src/nectarchain/user_scripts/ggrolleron/notebooks/waveforms_charges_extraction_new_framework.py @@ -0,0 +1,96 @@ +# %% +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 = 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 + +# %% +tool.initialize() + +# %% +tool.setup() + +# %% +tool.start() + +# %% +tool.finish() + +# %% +container = WaveformsContainers() +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) + container.containers[trigger] = next(tableReader) + except Exception as err: + print(err) + +# %% +container.containers + +# %% +container = ChargesContainers() +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) + container.containers[trigger] = next(tableReader) + except Exception as err: + print(err) + +# %% +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 + +# %% +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 5575a3d4..031c2841 100644 --- a/src/nectarchain/user_scripts/ggrolleron/test.py +++ b/src/nectarchain/user_scripts/ggrolleron/test.py @@ -1,82 +1,139 @@ -import numpy as np -#import pandas as pd +import logging +import os +import sys +import time +from pathlib import Path + +# import pandas as pd import matplotlib as mpl -import matplotlib.pyplot as plt import matplotlib.cm as cm -from pathlib import Path -import sys,os -import logging -logging.basicConfig(format='%(asctime)s %(name)s %(levelname)s %(message)s',level=logging.INFO) +import matplotlib.pyplot as plt +import numpy as np + +logging.basicConfig( + format="%(asctime)s %(name)s %(levelname)s %(message)s", level=logging.INFO +) log = logging.getLogger(__name__) ##tips to add message to stdout handler = logging.StreamHandler(sys.stdout) handler.setLevel(logging.INFO) -formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s') +formatter = logging.Formatter("%(asctime)s - %(name)s - %(levelname)s - %(message)s") handler.setFormatter(formatter) log.addHandler(handler) +import glob +from nectarchain.data.container import ( + ChargeContainer, + ChargeContainers, + WaveformsContainer, + WaveformsContainers, +) +from nectarchain.data.container.utils import DataManagement -from nectarchain.calibration.container import ChargeContainer,WaveformsContainer,WaveformsContainers,ChargeContainers -from nectarchain.calibration.container.utils import DataManagement -def test_multicontainers() : +def test_check_wfs(): run_number = [3731] - #waveforms = WaveformsContainers(run_number[0],max_events=20000) - #log.info("waveforms created") - #waveforms.load_wfs() - #log.info("waveforms loaded") + files = glob.glob( + f"{os.environ['NECTARCAMDATA']}/waveforms/waveforms_run{run_number[0]}_*.fits" + ) + if len(files) == 0: + raise FileNotFoundError(f"no splitted waveforms found") + else: + charge_container = ChargeContainers() + for j, file in enumerate(files): + log.debug(f"loading wfs from file {file}") + wfs = WaveformsContainer.load(file) + print(f"min : {wfs.wfs_hg.min()}") + # fig,ax = wfs.plot_waveform_hg(0) + # for i in range(wfs.nevents) : + # wfs.plot_waveform_hg(i,figure = fig,ax = ax) + + log.debug(f"computation of charge for file {file}") + charge = ChargeContainer.from_waveforms( + wfs, method="LocalPeakWindowSum", window_width=16, window_shift=4 + ) + hist = charge.histo_hg() + charge_container.append(charge) + + +def test_extractor(): + run_number = [2633] + + wfs = WaveformsContainer(run_number[0], max_events=200) # - #waveforms.write(f"{os.environ['NECTARCAMDATA']}/waveforms/",overwrite = True) - #log.info('waveforms written') + wfs.load_wfs() + # + # spe_run_1000V.write(f"{os.environ['NECTARCAMDATA']}/waveforms/",overwrite = True) - #waveforms = WaveformsContainers.load(f"{os.environ['NECTARCAMDATA']}/waveforms/waveforms_run{run_number[0]}") - #log.info("waveforms loaded") + # wfs = WaveformsContainer.load(f"{os.environ['NECTARCAMDATA']}/waveforms/waveforms_run{run_number[0]}.fits") - #charge = ChargeContainers.from_waveforms(waveforms,method = "LocalPeakWindowSum", window_width = 16, window_shift = 4) - #log.info("charge computed") - # - path = "LocalPeakWindowSum_4-12" - #charge.write(f"{os.environ['NECTARCAMDATA']}/charges/{path}/",overwrite = True) - #log.info("charge written") + # charge = ChargeContainer.compute_charge(spe_run_1000V,1,method = "gradient_extractor",) - charge = ChargeContainers.from_file(f"{os.environ['NECTARCAMDATA']}/charges/{path}/",run_number = run_number[0]) - log.info("charge loaded") + t = time.time() + charge = ChargeContainer.from_waveforms( + wfs, method="LocalPeakWindowSum", window_width=16, window_shift=4 + ) + log.info(f"LocalPeakWindowSum duration : {time.time() - t} seconds") - charge_merged = charge.merge() - log.info('charge merged') + charge = ChargeContainer.from_waveforms( + wfs, method="GlobalPeakWindowSum", window_width=16, window_shift=4 + ) + log.info(f"GlobalPeakWindowSum duration : {time.time() - t} seconds") + charge = ChargeContainer.from_waveforms( + wfs, method="FullWaveformSum", window_width=16, window_shift=4 + ) + log.info(f"FullWaveformSum duration : {time.time() - t} seconds") + charge = ChargeContainer.from_waveforms( + wfs, method="FixedWindowSum", window_width=16, window_shift=4, peak_index=30 + ) + log.info(f"FullWindowSum duration : {time.time() - t} seconds") + charge = ChargeContainer.from_waveforms( + wfs, method="SlidingWindowMaxSum", window_width=16, window_shift=4 + ) + log.info(f"SlidingWindowMaxSum duration : {time.time() - t} seconds") + # charge = ChargeContainer.from_waveforms(wfs, method = "NeighborPeakWindowSum", window_width = 16, window_shift = 4) + # log.info(f"NeighborPeakWindowSum duration : {time.time() - t} seconds") + # charge = ChargeContainer.from_waveforms(wfs, method = "BaselineSubtractedNeighborPeakWindowSum", baseline_start = 2, baseline_end = 12, window_width = 16, window_shift = 4) + # log.info(f"BaselineSubtractedNeighborPeakWindowSum duration : {time.time() - t} seconds") -""" -def test_simplecontainer() : - run_number = [2633] - ped_run_number = [2630] - FF_run_number = [2609] + charge = ChargeContainer.from_waveforms( + wfs, method="TwoPassWindowSum", window_width=16, window_shift=4 + ) + log.info(f"TwoPassWindowSum duration : {time.time() - t} seconds") - spe_run_1000V = WaveformsContainer(run_number[0],max_events = 1000) + # charge.write(f"{os.environ['NECTARCAMDATA']}/charges/std/",overwrite = False) - spe_run_1000V.load_wfs() - spe_run_1000V.write(f"{os.environ['NECTARCAMDATA']}/waveforms/",overwrite = True) +def test_simplecontainer(): + run_number = [2633] - spe_run_1000V.load(f"{os.environ['NECTARCAMDATA']}/waveforms/waveforms_run{run_number[0]}.fits") + spe_run_1000V = WaveformsContainer(run_number[0], max_events=1000) - charge = ChargeContainer.compute_charge(spe_run_1000V,1,method = "gradient_extractor") + spe_run_1000V.load_wfs() + # spe_run_1000V.write(f"{os.environ['NECTARCAMDATA']}/waveforms/",overwrite = True) + # spe_run_1000V.load(f"{os.environ['NECTARCAMDATA']}/waveforms/waveforms_run{run_number[0]}.fits") - charge = ChargeContainer.from_waveform(spe_run_1000V) + charge = ChargeContainer.compute_charge( + spe_run_1000V, + 1, + method="LocalPeakWindowSum", + extractor_kwargs={"window_width": 16, "window_shift": 4}, + ) + charge = ChargeContainer.from_waveforms(spe_run_1000V) - charge.write(f"{os.environ['NECTARCAMDATA']}/charges/std/",overwrite = True) + # charge.write(f"{os.environ['NECTARCAMDATA']}/charges/std/",overwrite = True) -""" -if __name__ == "__main__" : - test_multicontainers() +if __name__ == "__main__": + test_check_wfs() - print("work completed") \ No newline at end of file + print("work completed") diff --git a/src/nectarchain/user_scripts/jlenain/submitGainFitDIRAC/submitGainFitter.py b/src/nectarchain/user_scripts/jlenain/submitGainFitDIRAC/submitGainFitter.py index 3002a647..a65113f1 100755 --- a/src/nectarchain/user_scripts/jlenain/submitGainFitDIRAC/submitGainFitter.py +++ b/src/nectarchain/user_scripts/jlenain/submitGainFitDIRAC/submitGainFitter.py @@ -10,7 +10,7 @@ from astropy import time # nectarchain imports -from nectarchain.calibration.container.utils import DataManagement as dm +from nectarchain.data.container.utils import DataManagement as dm # DIRAC imports from DIRAC.Interfaces.API.Dirac import Dirac