Skip to content

Commit

Permalink
ChargesContainers I/O
Browse files Browse the repository at this point in the history 
cleaning
unit test for WaveformsMaker and ChargesMaker
  • Loading branch information
guillaume.grolleron committed Sep 15, 2023
1 parent fd0e78a commit f5d0180
Show file tree
Hide file tree
Showing 9 changed files with 400 additions and 240 deletions.
140 changes: 71 additions & 69 deletions src/nectarchain/data/container/chargesContainer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,10 @@

import numpy as np
from ctapipe.containers import Field
from abc import ABC
import os
from pathlib import Path
from astropy.io import fits

from .core import ArrayDataContainer

Expand All @@ -22,14 +26,13 @@ class ChargesContainer(ArrayDataContainer):
peak_lg = Field(
type = np.ndarray,
description = 'The low gain peak time')
_method = Field(
method = Field(
type = str,
description = 'The charge extraction method')
description = 'The charge extraction method used')


'''
def write(self,path : Path,**kwargs) :
class ChargesContainerIO(ABC) :
def write(path : Path, containers : ChargesContainer,**kwargs) :
"""method to write in an output FITS file the ChargeContainer.
Args:
Expand All @@ -41,64 +44,53 @@ def write(self,path : Path,**kwargs) :
log.info(f"saving in {path}")
os.makedirs(path,exist_ok = True)

#table = Table(self.charge_hg)
#table.meta["pixels_id"] = self._pixels_id
#table.write(Path(path)/f"charge_hg_run{self.run_number}.ecsv",format='ascii.ecsv',overwrite=kwargs.get('overwrite',False))
#
#table = Table(self.charge_lg)
#table.meta["pixels_id"] = self._pixels_id
#table.write(Path(path)/f"charge_lg_run{self.run_number}.ecsv",format='ascii.ecsv',overwrite=kwargs.get('overwrite',False))
#
#table = Table(self.peak_hg)
#table.meta["pixels_id"] = self._pixels_id
#table.write(Path(path)/f"peak_hg_run{self.run_number}.ecsv",format='ascii.ecsv',overwrite=kwargs.get('overwrite',False))
#
#table = Table(self.peak_lg)
#table.meta["pixels_id"] = self._pixels_id
#table.write(Path(path)/f"peak_lg_run{self.run_number}.ecsv",format='ascii.ecsv',overwrite=kwargs.get('overwrite',False))
hdr = fits.Header()
hdr['RUN'] = self._run_number
hdr['NEVENTS'] = self.nevents
hdr['NPIXELS'] = self.npixels
hdr['COMMENT'] = f"The charge containeur for run {self._run_number} with {self._method} method : primary is the pixels id, then you can find HG charge, LG charge, HG peak and LG peak, 2 last HDU are composed of event properties and trigger patern"
primary_hdu = fits.PrimaryHDU(self.pixels_id,header=hdr)
charge_hg_hdu = fits.ImageHDU(self.charge_hg,name = "HG charge")
charge_lg_hdu = fits.ImageHDU(self.charge_lg,name = "LG charge")
peak_hg_hdu = fits.ImageHDU(self.peak_hg, name = 'HG peak time')
peak_lg_hdu = fits.ImageHDU(self.peak_lg, name = 'LG peak time')
col1 = fits.Column(array = self.event_id, name = "event_id", format = '1I')
col2 = fits.Column(array = self.event_type, name = "event_type", format = '1I')
col3 = fits.Column(array = self.ucts_timestamp, name = "ucts_timestamp", format = '1K')
col4 = fits.Column(array = self.ucts_busy_counter, name = "ucts_busy_counter", format = '1I')
col5 = fits.Column(array = self.ucts_event_counter, name = "ucts_event_counter", format = '1I')
col6 = fits.Column(array = self.multiplicity, name = "multiplicity", format = '1I')
hdr['RUN'] = containers.run_number
hdr['NEVENTS'] = containers.nevents
hdr['NPIXELS'] = containers.npixels
hdr['METHOD'] = containers.method
hdr['CAMERA'] = containers.camera


hdr['COMMENT'] = f"The charge containeur for run {containers.run_number} with {containers.method} method : primary is the pixels id, then you can find HG charge, LG charge, HG peak and LG peak, 2 last HDU are composed of event properties and trigger patern"

primary_hdu = fits.PrimaryHDU(containers.pixels_id,header=hdr)
charge_hg_hdu = fits.ImageHDU(containers.charges_hg,name = "HG charge")
charge_lg_hdu = fits.ImageHDU(containers.charges_lg,name = "LG charge")
peak_hg_hdu = fits.ImageHDU(containers.peak_hg, name = 'HG peak time')
peak_lg_hdu = fits.ImageHDU(containers.peak_lg, name = 'LG peak time')

col1 = fits.Column(array = containers.broken_pixels_hg, name = "HG broken pixels", format = f'{containers.broken_pixels_hg.shape[1]}L')
col2 = fits.Column(array = containers.broken_pixels_lg, name = "LG broken pixels", format = f'{containers.broken_pixels_lg.shape[1]}L')
coldefs = fits.ColDefs([col1, col2])
broken_pixels = fits.BinTableHDU.from_columns(coldefs,name = 'trigger patern')

col1 = fits.Column(array = containers.event_id, name = "event_id", format = '1I')
col2 = fits.Column(array = containers.event_type, name = "event_type", format = '1I')
col3 = fits.Column(array = containers.ucts_timestamp, name = "ucts_timestamp", format = '1K')
col4 = fits.Column(array = containers.ucts_busy_counter, name = "ucts_busy_counter", format = '1I')
col5 = fits.Column(array = containers.ucts_event_counter, name = "ucts_event_counter", format = '1I')
col6 = fits.Column(array = containers.multiplicity, name = "multiplicity", format = '1I')

coldefs = fits.ColDefs([col1, col2, col3, col4, col5, col6])
event_properties = fits.BinTableHDU.from_columns(coldefs, name = 'event properties')

col1 = fits.Column(array = self.trig_pattern_all, name = "trig_pattern_all", format = f'{4 * self.CAMERA.n_pixels}L',dim = f'({self.CAMERA.n_pixels},4)')
col2 = fits.Column(array = self.trig_pattern, name = "trig_pattern", format = f'{self.CAMERA.n_pixels}L')
col1 = fits.Column(array = containers.trig_pattern_all, name = "trig_pattern_all", format = f'{4 * containers.trig_pattern_all.shape[1]}L',dim = f'({ containers.trig_pattern_all.shape[1]},4)')
col2 = fits.Column(array = containers.trig_pattern, name = "trig_pattern", format = f'{containers.trig_pattern_all.shape[1]}L')
coldefs = fits.ColDefs([col1, col2])
trigger_patern = fits.BinTableHDU.from_columns(coldefs, name = 'trigger patern')

hdul = fits.HDUList([primary_hdu, charge_hg_hdu, charge_lg_hdu,peak_hg_hdu,peak_lg_hdu,event_properties,trigger_patern])
hdul = fits.HDUList([primary_hdu, charge_hg_hdu, charge_lg_hdu,peak_hg_hdu,peak_lg_hdu,broken_pixels,event_properties,trigger_patern])
try :
hdul.writeto(Path(path)/f"charge_run{self.run_number}{suffix}.fits",overwrite=kwargs.get('overwrite',False))
log.info(f"charge saved in {Path(path)}/charge_run{self.run_number}{suffix}.fits")
hdul.writeto(Path(path)/f"charge_run{containers.run_number}{suffix}.fits",overwrite=kwargs.get('overwrite',False))
log.info(f"charge saved in {Path(path)}/charge_run{containers.run_number}{suffix}.fits")
except OSError as e :
log.warning(e)
except Exception as e :
log.error(e,exc_info = True)
raise e

@staticmethod
def from_file(path : Path,run_number : int,**kwargs) :
def load(path : Path,run_number : int,**kwargs) :
"""load ChargeContainer from FITS file previously written with ChargeContainer.write() method
Args:
Expand All @@ -116,24 +108,34 @@ def from_file(path : Path,run_number : int,**kwargs) :
filename = Path(path)/f"charge_run{run_number}.fits"

with fits.open(filename) as hdul :
pixels_id = hdul[0].data
nevents = hdul[0].header['NEVENTS']
npixels = hdul[0].header['NPIXELS']
charge_hg = hdul[1].data
charge_lg = hdul[2].data
peak_hg = hdul[3].data
peak_lg = hdul[4].data
cls = ChargeContainer(charge_hg,charge_lg,peak_hg,peak_lg,run_number,pixels_id,nevents,npixels)
cls.event_id = hdul[5].data["event_id"]
cls.event_type = hdul[5].data["event_type"]
cls.ucts_timestamp = hdul[5].data["ucts_timestamp"]
cls.ucts_busy_counter = hdul[5].data["ucts_busy_counter"]
cls.ucts_event_counter = hdul[5].data["ucts_event_counter"]
table_trigger = hdul[6].data
cls.trig_pattern_all = table_trigger["trig_pattern_all"]
return cls
'''
containers = ChargesContainer()
containers.run_number = hdul[0].header['RUN']
containers.nevents = hdul[0].header['NEVENTS']
containers.npixels = hdul[0].header['NPIXELS']
containers.method = hdul[0].header['METHOD']
containers.camera = hdul[0].header['CAMERA']

containers.pixels_id = hdul[0].data
containers.charges_hg = hdul[1].data
containers.charges_lg = hdul[2].data
containers.peak_hg = hdul[3].data
containers.peak_lg = hdul[4].data

broken_pixels = hdul[5].data
containers.broken_pixels_hg = broken_pixels["HG broken pixels"]
containers.broken_pixels_lg = broken_pixels["LG broken pixels"]



table_prop = hdul[6].data
containers.event_id = table_prop["event_id"]
containers.event_type = table_prop["event_type"]
containers.ucts_timestamp = table_prop["ucts_timestamp"]
containers.ucts_busy_counter = table_prop["ucts_busy_counter"]
containers.ucts_event_counter = table_prop["ucts_event_counter"]
containers.multiplicity = table_prop["multiplicity"]

table_trigger = hdul[7].data
containers.trig_pattern_all = table_trigger["trig_pattern_all"]
containers.trig_pattern = table_trigger["trig_pattern"]
return containers
1 change: 0 additions & 1 deletion src/nectarchain/data/container/core.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,6 @@
log = logging.getLogger(__name__)
log.handlers = logging.getLogger('__main__').handlers

import os
from ctapipe.containers import Container,Field
import numpy as np

Expand Down
62 changes: 30 additions & 32 deletions src/nectarchain/makers/chargesMakers.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,19 +6,19 @@
from argparse import ArgumentError
import numpy as np
import numpy.ma as ma
from pathlib import Path
import glob
import time
import os

from ctapipe.instrument import CameraGeometry

from ctapipe_io_nectarcam import constants
from ctapipe.containers import EventType
from ctapipe.image import ImageExtractor

from ctapipe.image.extractor import (FullWaveformSum,
FixedWindowSum,
GlobalPeakWindowSum,
LocalPeakWindowSum,
SlidingWindowMaxSum,
NeighborPeakWindowSum,
BaselineSubtractedNeighborPeakWindowSum,
TwoPassWindowSum)

from astropy.io import fits

from numba import guvectorize, float64, int64, bool_

Expand All @@ -28,9 +28,7 @@
from .extractor.utils import CtapipeExtractor




__all__ = ['ChargeContainer']
__all__ = ['ChargesMaker']

list_ctapipe_charge_extractor = ["FullWaveformSum",
"FixedWindowSum",
Expand Down Expand Up @@ -98,7 +96,7 @@ def make_histo(charge, all_range, mask_broken_pix, _mask, hist_ma_data):
pass


class ChargeMaker(ArrayDataMaker) :
class ChargesMaker(ArrayDataMaker) :
#constructors
def __init__(self,run_number : int,max_events : int = None,run_file = None,*args,**kwargs):
"""construtor
Expand All @@ -118,7 +116,7 @@ def __init__(self,run_number : int,max_events : int = None,run_file = None,*args
def _init_trigger_type(self,trigger_type,**kwargs) :
super()._init_trigger_type(trigger_type,**kwargs)
name = __class__._get_name_trigger(trigger_type)
log.info(f"initialization of the chargeMaker following trigger type : {name}")
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}"] = []
Expand All @@ -132,7 +130,7 @@ def make(self,
method: str = "FullWaveformSum",
*args,**kwargs):
kwargs["method"]=method
super().make(n_events=n_events,
return super().make(n_events=n_events,
trigger_type=trigger_type,
restart_from_begining=restart_from_begining,
*args,**kwargs)
Expand All @@ -144,18 +142,14 @@ def _make_event(self,
*args,
**kwargs
) :
wfs_hg_tmp=np.zeros((self.npixels,self.nsamples),dtype = np.uint16)
wfs_lg_tmp=np.zeros((self.npixels,self.nsamples),dtype = np.uint16)

super()._make_event(event = event,
trigger = trigger,
wfs_hg = wfs_hg_tmp,
wfs_lg = wfs_lg_tmp,
*args,
**kwargs)

wfs_hg_tmp,wfs_lg_tmp = super()._make_event(event = event,
trigger = trigger,
return_wfs = True,
*args,**kwargs)
name = __class__._get_name_trigger(trigger)

broken_pixels_hg,broken_pixels_lg = __class__._compute_broken_pixels(wfs_hg_tmp,wfs_lg_tmp)
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())

Expand Down Expand Up @@ -187,14 +181,15 @@ def _get_imageExtractor(method,subarray,**kwargs) :
imageExtractor = eval(method)(subarray,**extractor_kwargs)
return imageExtractor

def _make_output_container(self,trigger_type) :
def _make_output_container(self,trigger_type,method : str) :
output = []
for trigger in trigger_type :
chargesContainer = ChargesContainer(
run_number = self.run_number,
npixels = self.npixels,
run_number = self._run_number,
npixels = self._npixels,
camera = self.CAMERA_NAME,
pixels_id = self.pixels_id,
pixels_id = self._pixels_id,
method = method,
nevents = self.nevents(trigger),
charges_hg = self.charges_hg(trigger),
charges_lg = self.charges_lg(trigger),
Expand All @@ -221,7 +216,9 @@ def sort(chargesContainer :ChargesContainer, method = 'event_id') :
npixels = chargesContainer.npixels,
camera = chargesContainer.camera,
pixels_id = chargesContainer.pixels_id,
nevents = chargesContainer.nevents
nevents = chargesContainer.nevents,
method = chargesContainer.method

)
if method == 'event_id' :
index = np.argsort(chargesContainer.event_id)
Expand Down Expand Up @@ -261,7 +258,7 @@ def create_from_waveforms(waveformsContainer: WaveformsContainer, method: str =
for field in waveformsContainer.keys() :
if not(field in ["subarray","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)
Expand All @@ -272,6 +269,7 @@ def create_from_waveforms(waveformsContainer: WaveformsContainer, method: str =
chargesContainer.charges_lg = charges_lg
chargesContainer.peak_hg = peak_hg
chargesContainer.peak_lg = peak_lg
chargesContainer.method = method

return chargesContainer

Expand All @@ -298,10 +296,10 @@ def compute_charge(waveformContainer : WaveformsContainer,channel : int,method :
imageExtractor = __class__._get_imageExtractor(method = method,subarray = waveformContainer.subarray,**kwargs)

if channel == constants.HIGH_GAIN:
out = np.array([CtapipeExtractor.get_image_peak_time(imageExtractor(waveformContainer.wfs_hg[i],waveformContainer.TEL_ID,channel,waveformContainer.broken_pixels_hg)) for i in range(len(waveformContainer.wfs_hg))]).transpose(1,0,2)
out = np.array([CtapipeExtractor.get_image_peak_time(imageExtractor(waveformContainer.wfs_hg[i],__class__.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],waveformContainer.TEL_ID,channel,waveformContainer.broken_pixels_lg)) for i in range(len(waveformContainer.wfs_lg))]).transpose(1,0,2)
out = np.array([CtapipeExtractor.get_image_peak_time(imageExtractor(waveformContainer.wfs_lg[i],__class__.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}")
Expand Down
Loading

0 comments on commit f5d0180

Please sign in to comment.