diff --git a/.gitignore b/.gitignore index d1e9305f..84752161 100644 --- a/.gitignore +++ b/.gitignore @@ -86,4 +86,11 @@ target notebooks/*.html notebooks/*.png +nectarchain/calibration/test +**.log +**.pdf +**.csv + +#VScode +.vscode/ diff --git a/nectarchain/calibration/README.md b/nectarchain/calibration/README.md new file mode 100644 index 00000000..b48f235e --- /dev/null +++ b/nectarchain/calibration/README.md @@ -0,0 +1,7 @@ +Environment variable needed: +export NECTARCAMDATA="path to local NectarCam data, it can contain fits.fz run files, WaveformsContainer or ChargeContainer FITS files" + +Environment variables which can be defined: +export NECTARCHAIN_TEST="path to test for nectarchain" +export NECTARCHAIN_LOG="path to log for nectarchain" +export NECTARCHAIN_FIGURES="path to log figures for nectarchain" \ No newline at end of file diff --git a/nectarchain/calibration/container/__init__.py b/nectarchain/calibration/container/__init__.py new file mode 100644 index 00000000..8eae44c4 --- /dev/null +++ b/nectarchain/calibration/container/__init__.py @@ -0,0 +1,2 @@ +from .waveforms import * +from .charge import * \ No newline at end of file diff --git a/nectarchain/calibration/container/charge.py b/nectarchain/calibration/container/charge.py new file mode 100644 index 00000000..2707667c --- /dev/null +++ b/nectarchain/calibration/container/charge.py @@ -0,0 +1,348 @@ +from argparse import ArgumentError +import numpy as np +import numpy.ma as ma +from matplotlib import pyplot as plt +import copy +from pathlib import Path +import sys +import os +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.visualization import CameraDisplay +from ctapipe.coordinates import CameraFrame,EngineeringCameraFrame +from ctapipe.instrument import CameraGeometry +from ctapipe.image.extractor import (FullWaveformSum, + FixedWindowSum, + GlobalPeakWindowSum, + LocalPeakWindowSum, + SlidingWindowMaxSum, + NeighborPeakWindowSum, + BaselineSubtractedNeighborPeakWindowSum, + TwoPassWindowSum) + +from ctapipe_io_nectarcam import NectarCAMEventSource,constants +from ctapipe.io import EventSource, EventSeeker + +from astropy.table import Table +from astropy.io import fits + +from .waveforms import WaveformsContainer,WaveformsContainers + + + +#from .charge_extractor import * + +__all__ = ['ChargeContainer','ChargeContainers'] + +list_ctapipe_charge_extractor = ["FullWaveformSum", + "FixedWindowSum", + "GlobalPeakWindowSum", + "LocalPeakWindowSum", + "SlidingWindowMaxSum", + "NeighborPeakWindowSum", + "BaselineSubtractedNeighborPeakWindowSum", + "TwoPassWindowSum"] + +list_nectarchain_charge_extractor = ['gradient_extractor'] + + + + +class ChargeContainer() : + """class used to compute charge from waveforms container""" + TEL_ID = 0 + CAMERA = CameraGeometry.from_name("NectarCam-003") + + def __init__(self,charge_hg,charge_lg,peak_hg,peak_lg,run_number,pixels_id,nevents,npixels, method = "FullWaveformSum") : + self.charge_hg = charge_hg + self.charge_lg = charge_lg + self.peak_hg = peak_hg + self.peak_lg = peak_lg + self.__run_number = run_number + self.__pixels_id = pixels_id + self.__method = method + self.__nevents = nevents + self.__npixels = npixels + + + self.ucts_timestamp = np.zeros((self.nevents),dtype = np.uint64) + self.ucts_busy_counter = np.zeros((self.nevents),dtype = np.uint16) + self.ucts_event_counter = np.zeros((self.nevents),dtype = np.uint16) + self.event_type = np.zeros((self.nevents),dtype = np.uint8) + self.event_id = np.zeros((self.nevents),dtype = np.uint16) + self.trig_pattern_all = np.zeros((self.nevents,self.CAMERA.n_pixels,4),dtype = bool) + + @classmethod + def from_waveforms(cls,waveformContainer : WaveformsContainer,method : str = "FullWaveformSum",**kwargs) : + log.info(f"computing hg charge with {method} method") + charge_hg,peak_hg = ChargeContainer.compute_charge(waveformContainer,constants.HIGH_GAIN,method,**kwargs) + charge_hg = np.array(charge_hg,dtype = np.uint16) + log.info(f"computing lg charge with {method} method") + charge_lg,peak_lg = ChargeContainer.compute_charge(waveformContainer,constants.LOW_GAIN,method,**kwargs) + charge_lg = np.array(charge_lg,dtype = np.uint16) + + chargeContainer = cls(charge_hg,charge_lg,peak_hg,peak_lg,waveformContainer.run_number,waveformContainer.pixels_id,waveformContainer.nevents,waveformContainer.npixels ,method) + + chargeContainer.ucts_timestamp = waveformContainer.ucts_timestamp + chargeContainer.ucts_busy_counter = waveformContainer.ucts_busy_counter + chargeContainer.ucts_event_counter = waveformContainer.ucts_event_counter + chargeContainer.event_type = waveformContainer.event_type + chargeContainer.event_id = waveformContainer.event_id + chargeContainer.trig_pattern_all = waveformContainer.trig_pattern_all + + return chargeContainer + + + def write(self,path : Path,**kwargs) : + suffix = kwargs.get("suffix","") + if suffix != "" : suffix = f"_{suffix}" + + 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) + charge_lg_hdu = fits.ImageHDU(self.charge_lg) + peak_hg_hdu = fits.ImageHDU(self.peak_hg) + peak_lg_hdu = fits.ImageHDU(self.peak_lg) + + 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') + + coldefs = fits.ColDefs([col1, col2, col3, col4, col5, col6]) + event_properties = fits.BinTableHDU.from_columns(coldefs) + + 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') + coldefs = fits.ColDefs([col1, col2]) + trigger_patern = fits.BinTableHDU.from_columns(coldefs) + + hdul = fits.HDUList([primary_hdu, charge_hg_hdu, charge_lg_hdu,peak_hg_hdu,peak_lg_hdu,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") + 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) : + log.info(f"loading in {path} run number {run_number}") + + #table = Table.read(Path(path)/f"charge_hg_run{run_number}.ecsv") + #pixels_id = table.meta['pixels_id'] + #charge_hg = np.array([table[colname] for colname in table.colnames]).T + # + #table = Table.read(Path(path)/f"charge_lg_run{run_number}.ecsv") + #charge_lg = np.array([table[colname] for colname in table.colnames]).T + # + #table = Table.read(Path(path)/f"peak_hg_run{run_number}.ecsv") + #peak_hg = np.array([table[colname] for colname in table.colnames]).T + # + #table = Table.read(Path(path)/f"peak_lg_run{run_number}.ecsv") + #peak_lg = np.array([table[colname] for colname in table.colnames]).T + # + hdul = fits.open(Path(path)/f"charge_run{run_number}.fits") + 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 + + + @staticmethod + def compute_charge(waveformContainer : WaveformsContainer,channel : int,method : str = "FullWaveformSum" ,**kwargs) : + 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}") + ImageExtractor = eval(method)(waveformContainer.subarray,**kwargs) + if channel == constants.HIGH_GAIN: + return ImageExtractor(waveformContainer.wfs_hg,waveformContainer.TEL_ID,channel) + elif channel == constants.LOW_GAIN: + return ImageExtractor(waveformContainer.wfs_lg,waveformContainer.TEL_ID,channel) + else : + raise ArgumentError("channel must be 0 or 1") + + def histo_hg(self,n_bins : int = 1000,autoscale : bool = True) -> np.ndarray: + mask_broken_pix = np.array((self.charge_hg == self.charge_hg.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(self.charge_hg.T[~mask_broken_pix].T)) - 0.5,np.uint16(np.max(self.charge_hg.T[~mask_broken_pix].T)) + 1.5,1) + hist_ma = ma.masked_array(np.zeros((self.charge_hg.shape[1],all_range.shape[0]),dtype = np.uint16), mask=np.zeros((self.charge_hg.shape[1],all_range.shape[0]),dtype = bool)) + charge_ma = ma.masked_array(np.zeros((self.charge_hg.shape[1],all_range.shape[0])), mask=np.zeros((self.charge_hg.shape[1],all_range.shape[0]),dtype = bool)) + + new_data_mask = np.array([np.logical_or(charge_ma.mask.T[i],mask_broken_pix) for i in range(charge_ma.shape[1])]) + charge_ma.mask = new_data_mask + hist_ma.mask = new_data_mask + + + for i in range(self.charge_hg.shape[1]) : + log.debug(f'computing charge histogram for pixel {i}') + if mask_broken_pix[i] : + log.debug('This pixel is broken, I mask this one') + hist_ma.mask[i] = np.ones(hist_ma.mask[i].shape) + charge_ma.mask[i] = np.ones(charge_ma.mask[i].shape) + else : + log.debug("this pixel is not broken, let's continue computation") + hist,charge = np.histogram(self.charge_hg.T[i],bins=np.arange(np.uint16(np.min(self.charge_hg.T[i])) - 1, np.uint16(np.max(self.charge_hg.T[i])) + 2,1)) + charge_edges = np.array([np.mean(charge[i:i+2],axis = 0) for i in range(charge.shape[0]-1)]) + mask = (all_range >= charge_edges[0]) * (all_range <= charge_edges[-1]) + + #MASK THE DATA + hist_ma.mask[i] = np.logical_or(~mask, hist_ma.mask[i]) + charge_ma.mask[i] = np.logical_or(~mask, charge_ma.mask[i]) + + #FILL THE DATA + hist_ma.data[i][mask] = hist + charge_ma.data[i] = all_range + + return ma.masked_array((hist_ma,charge_ma)) + + + else : + hist = np.array([np.histogram(self.charge_hg.T[i],bins=n_bins)[0] for i in range(self.charge_hg.shape[1])]) + charge = np.array([np.histogram(self.charge_hg.T[i],bins=n_bins)[1] for i in range(self.charge_hg.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)) + + def histo_lg(self,n_bins: int = 1000,autoscale : bool = True) -> np.ndarray: + mask_broken_pix = np.array((self.charge_lg == self.charge_lg.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(self.charge_lg.T[~mask_broken_pix].T)) - 0.5,np.uint16(np.max(self.charge_lg.T[~mask_broken_pix].T)) + 1.5,1) + hist_ma = ma.masked_array(np.zeros((self.charge_lg.shape[1],all_range.shape[0]),dtype = np.uint16), mask=np.zeros((self.charge_lg.shape[1],all_range.shape[0]),dtype = bool)) + charge_ma = ma.masked_array(np.zeros((self.charge_lg.shape[1],all_range.shape[0])), mask=np.zeros((self.charge_lg.shape[1],all_range.shape[0]),dtype = bool)) + + new_data_mask = np.array([np.logical_or(charge_ma.mask.T[i],mask_broken_pix) for i in range(charge_ma.shape[1])]) + charge_ma.mask = new_data_mask + hist_ma.mask = new_data_mask + + + for i in range(self.charge_lg.shape[1]) : + log.debug(f'computing charge histogram for pixel {i}') + if mask_broken_pix[i] : + log.debug('This pixel is broken, I mask this one') + hist_ma.mask[i] = np.ones(hist_ma.mask[i].shape) + charge_ma.mask[i] = np.ones(charge_ma.mask[i].shape) + else : + log.debug("this pixel is not broken, let's continue computation") + hist,charge = np.histogram(self.charge_lg.T[i],bins=np.arange(np.uint16(np.min(self.charge_lg.T[i])) - 1, np.uint16(np.max(self.charge_lg.T[i])) + 2,1)) + charge_edges = np.array([np.mean(charge[i:i+2],axis = 0) for i in range(charge.shape[0]-1)]) + mask = (all_range >= charge_edges[0]) * (all_range <= charge_edges[-1]) + + #MASK THE DATA + hist_ma.mask[i] = np.logical_or(~mask, hist_ma.mask[i]) + charge_ma.mask[i] = np.logical_or(~mask, charge_ma.mask[i]) + + #FILL THE DATA + hist_ma.data[i][mask] = hist + charge_ma.data[i] = all_range + + else : + hist = np.array([np.histogram(self.charge_lg.T[i],bins=n_bins)[0] for i in range(self.charge_lg.shape[1])]) + charge = np.array([np.histogram(self.charge_lg.T[i],bins=n_bins)[1] for i in range(self.charge_lg.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)) + + @property + def run_number(self) : return self.__run_number + + @property + def pixels_id(self) : return self.__pixels_id + + @property + def npixels(self) : return self.__npixels + + @property + def nevents(self) : return self.__nevents + + + #physical properties + @property + def multiplicity(self) : return np.uint16(np.count_nonzero(self.trig_pattern,axis = 1)) + + @property + def trig_pattern(self) : return self.trig_pattern_all.any(axis = 1) + + + + +class ChargeContainers() : + def __init__(self, *args, **kwargs) : + self.chargeContainers = [] + self.__nChargeContainer = 0 + + @classmethod + def from_waveforms(cls,waveformContainers : WaveformsContainers,**kwargs) : + chargeContainers = cls() + for i in range(waveformContainers.nWaveformsContainer) : + chargeContainers.append(ChargeContainer.from_waveforms(waveformContainers.waveformsContainer[i])) + return chargeContainers + + def write(self,path : str, **kwargs) : + for i in range(self.__nChargeContainer) : + self.chargeContainers[i].write(path,suffix = f"{i:04d}" ,**kwargs) + + @property + def nChargeContainer(self) : return self.__nChargeContainer + + @property + def nevents(self) : + return np.sum([self.chargeContainers[i].nevents for i in range(self.__nChargeContainer)]) + + def append(self,chargeContainer : ChargeContainer) : + self.chargeContainers.append(chargeContainer) + self.__nChargeContainer += 1 \ No newline at end of file diff --git a/nectarchain/calibration/container/charge_extractor.py b/nectarchain/calibration/container/charge_extractor.py new file mode 100644 index 00000000..fc851fb1 --- /dev/null +++ b/nectarchain/calibration/container/charge_extractor.py @@ -0,0 +1,95 @@ +from ctapipe.image import ImageExtractor +import numpy as np +import logging +from scipy.interpolate import InterpolatedUnivariateSpline +from scipy.signal import find_peaks + +from numba import guvectorize + +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 \ No newline at end of file diff --git a/nectarchain/calibration/container/utils.py b/nectarchain/calibration/container/utils.py new file mode 100644 index 00000000..feeffc08 --- /dev/null +++ b/nectarchain/calibration/container/utils.py @@ -0,0 +1,177 @@ +from email.generator import Generator +import os +import glob +import re +import mechanize +import requests +import browser_cookie3 + +from DIRAC.Interfaces.API.Dirac import Dirac +from pathlib import Path +from typing import List,Tuple + + +import logging +logging.basicConfig(format='%(asctime)s %(name)s %(levelname)s %(message)s') +log = logging.getLogger(__name__) +log.handlers = logging.getLogger('__main__').handlers + +__all__ = ['DataManagment','ChainGenerator'] + +class DataManagment() : + @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 = DataManagment.get_GRID_location(run_number) + DataManagment.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 + """ + + 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_ + """ + + 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}" + + #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=&ModuleCount=&subtext=%23{run_number}',cookies = cookies) + + if "ELOG Login" in req.text : + log.debug('log to Elog with cookies impossible, try again with 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=&ModuleCount=&subtext=%23{run_number}',cookies = cookies) + + + lines = req.text.split('\r\n') + + url_data = None + for line in lines : + if '

' in line : + url_data = line.split("

")[0].split('FC:')[1] + break + + if output_lfns : + try : + #Dirac + dirac = Dirac() + loc = f"/vo.cta.in2p3.fr/nectarcam/{url_data.split('/')[-2]}/{url_data.split('/')[-1]}" + res = dirac.listCatalogDirectory(loc, printOutput=True) + + lfns = [] + 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:])) + diff --git a/nectarchain/calibration/container/waveforms.py b/nectarchain/calibration/container/waveforms.py new file mode 100644 index 00000000..48067026 --- /dev/null +++ b/nectarchain/calibration/container/waveforms.py @@ -0,0 +1,360 @@ +from argparse import ArgumentError +import numpy as np +from matplotlib import pyplot as plt +import copy +import os +from pathlib import Path + +from enum import Enum + +from astropy.io import fits +from astropy.table import QTable,Column,Table +import astropy.units as u + +from ctapipe.visualization import CameraDisplay +from ctapipe.coordinates import CameraFrame,EngineeringCameraFrame +from ctapipe.instrument import CameraGeometry,SubarrayDescription,TelescopeDescription + +from ctapipe_io_nectarcam import NectarCAMEventSource +from ctapipe.containers import EventType +from ctapipe.io import EventSource, EventSeeker + +from .utils import DataManagment,ChainGenerator + +import sys +import logging +logging.basicConfig(format='%(asctime)s %(name)s %(levelname)s %(message)s') +log = logging.getLogger(__name__) +log.handlers = logging.getLogger('__main__').handlers + +__all__ = ["WaveformsContainer","WaveformsContainers"] + + +class WaveformsContainer() : + """class used to load run and load waveforms from r0 data + """ + TEL_ID = 0 + CAMERA = CameraGeometry.from_name("NectarCam-003") + def __new__(cls,*args,**kwargs) : + obj = object.__new__(cls) + return obj + + def __init__(self,run_number : int,max_events : int = None,nevents : int = -1,run_file = None): + """construtor + + Args: + run_number (int): id of the run to be loaded + maxevents (int, optional): max of events to be loaded. Defaults to 0, to load everythings. + nevents (int, optional) : number of events in run if known (parameter used to save computing time) + merge_file (optional) : if True will load all fits.fz files of the run, else merge_file can be integer to select the run fits.fz file according to this number + + """ + + self.__run_number = run_number + #gerer ici le fait de traiter plusieurs fichiers ou simplement 1 par 1 + self.__reader = WaveformsContainer.load_run(run_number,max_events,run_file = run_file) + + #from reader members + self.__npixels = self.__reader.camera_config.num_pixels + self.__nsamples = self.__reader.camera_config.num_samples + self.__geometry = self.__reader.subarray.tel[WaveformsContainer.TEL_ID].camera + self.__subarray = self.__reader.subarray + self.__pixels_id = self.__reader.camera_config.expected_pixels_id + + #set camera properties + log.info(f"N pixels : {self.npixels}") + + #run properties + if nevents != -1 : + self.__nevents = nevents if max_events is None else min(max_events,nevents) #self.check_events() + else : + self.__nevents = self.check_events() + #reload file (bc check_events has drained reader generator) + self.__reader = WaveformsContainer.load_run(run_number,max_events,run_file = run_file) + log.info(f"N_events : {self.nevents}") + + #define zeros members which will be filled therafter + self.wfs_hg = np.zeros((self.nevents,self.npixels,self.nsamples),dtype = np.uint16) + self.wfs_lg = np.zeros((self.nevents,self.npixels,self.nsamples),dtype = np.uint16) + self.ucts_timestamp = np.zeros((self.nevents),dtype = np.uint64) + self.ucts_busy_counter = np.zeros((self.nevents),dtype = np.uint32) + self.ucts_event_counter = np.zeros((self.nevents),dtype = np.uint32) + self.event_type = np.zeros((self.nevents),dtype = np.uint8) + self.event_id = np.zeros((self.nevents),dtype = np.uint32) + self.trig_pattern_all = np.zeros((self.nevents,self.CAMERA.n_pixels,4),dtype = bool) + #self.trig_pattern = np.zeros((self.nevents,self.npixels),dtype = bool) + #self.multiplicity = np.zeros((self.nevents,self.npixels),dtype = np.uint16) + + + @staticmethod + def load_run(run_number : int,max_events : int = None, run_file = None) : + """Static method to load from $NECTARCAMDATA directory data for specified run with max_events + + Args: + run_number (int): run_id + maxevents (int, optional): max of events to be loaded. Defaults to 0, to load everythings. + merge_file (optional) : if True will load all fits.fz files of the run, else merge_file can be integer to select the run fits.fz file according to this number + Returns: + List[ctapipe_io_nectarcam.NectarCAMEventSource]: List of EventSource for each run files + """ + generic_filename,filenames = DataManagment.findrun(run_number) + if run_file is None : + 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 + + def check_events(self): + """method to check triggertype of events and counting number of events in all readers + it prints the trigger type when trigger type change over events (to not write one line by events) + + Returns: + int : number of events + """ + log.info("checking trigger type") + has_changed = True + previous_trigger = None + n_events = 0 + for i,event in enumerate(self.__reader): + if previous_trigger is not None : + has_changed = (previous_trigger.value != event.trigger.event_type.value) + previous_trigger = event.trigger.event_type + if has_changed : + log.info(f"event {i} is {previous_trigger} : {previous_trigger.value} trigger type") + n_events+=1 + return n_events + + + + def load_wfs(self,compute_trigger_patern = False): + wfs_hg_tmp=np.zeros((self.npixels,self.nsamples),dtype = np.uint16) + wfs_lg_tmp=np.zeros((self.npixels,self.nsamples),dtype = np.uint16) + + n_traited_events = 0 + for i,event in enumerate(self.__reader): + if i%100 == 0: + log.info(f"reading event number {i}") + + self.event_id[i] = np.uint16(event.index.event_id) + self.ucts_timestamp[i]=event.nectarcam.tel[WaveformsContainer.TEL_ID].evt.ucts_timestamp + self.event_type[i]=event.trigger.event_type.value + self.ucts_busy_counter[i] = event.nectarcam.tel[WaveformsContainer.TEL_ID].evt.ucts_busy_counter + self.ucts_event_counter[i] = event.nectarcam.tel[WaveformsContainer.TEL_ID].evt.ucts_event_counter + + + self.trig_pattern_all[i] = event.nectarcam.tel[WaveformsContainer.TEL_ID].evt.trigger_pattern.T + + for pix in range(self.npixels): + wfs_lg_tmp[pix]=event.r0.tel[0].waveform[1,self.pixels_id[pix]] + wfs_hg_tmp[pix]=event.r0.tel[0].waveform[0,self.pixels_id[pix]] + + self.wfs_hg[i] = wfs_hg_tmp + self.wfs_lg[i] = wfs_lg_tmp + + + + #if compute_trigger_patern and np.max(self.trig_pattern) == 0: + # self.compute_trigger_patern() + + + def write(self,path : str, **kwargs) : + """method to write in an output FITS file the WaveformsContainer. Two files are created, one FITS representing the data + and one HDF5 file representing the subarray configuration + + Args: + path (str): the directory where you want to save data + """ + suffix = kwargs.get("suffix","") + if suffix != "" : suffix = f"_{suffix}" + + log.info(f"saving in {path}") + os.makedirs(path,exist_ok = True) + + hdr = fits.Header() + hdr['RUN'] = self.__run_number + hdr['NEVENTS'] = self.nevents + hdr['NPIXELS'] = self.npixels + hdr['NSAMPLES'] = self.nsamples + hdr['SUBARRAY'] = self.subarray.name + + self.subarray.to_hdf(f"{Path(path)}/subarray_run{self.run_number}.hdf5",overwrite=kwargs.get('overwrite',False)) + + + + hdr['COMMENT'] = f"The waveforms containeur for run {self.__run_number} : primary is the pixels id, 2nd HDU : high gain waveforms, 3rd HDU : low gain waveforms, 4th HDU : event properties and 5th HDU trigger paterns." + + primary_hdu = fits.PrimaryHDU(self.pixels_id,header=hdr) + + wfs_hg_hdu = fits.ImageHDU(self.wfs_hg) + wfs_lg_hdu = fits.ImageHDU(self.wfs_lg) + + + col1 = fits.Column(array = self.event_id, name = "event_id", format = '1J') + 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 = '1J') + col5 = fits.Column(array = self.ucts_event_counter, name = "ucts_event_counter", format = '1J') + col6 = fits.Column(array = self.multiplicity, name = "multiplicity", format = '1I') + + coldefs = fits.ColDefs([col1, col2, col3, col4, col5, col6]) + event_properties = fits.BinTableHDU.from_columns(coldefs) + + 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') + coldefs = fits.ColDefs([col1, col2]) + trigger_patern = fits.BinTableHDU.from_columns(coldefs) + + hdul = fits.HDUList([primary_hdu, wfs_hg_hdu, wfs_lg_hdu,event_properties,trigger_patern]) + try : + hdul.writeto(Path(path)/f"waveforms_run{self.run_number}{suffix}.fits",overwrite=kwargs.get('overwrite',False)) + log.info(f"run saved in {Path(path)}/waveforms_run{self.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 load(path : str) : + """load WaveformsContainer from FITS file previously written with WaveformsContainer.write() method + Note : 2 files are loaded, the FITS one representing the waveforms data and a HDF5 file representing the subarray configuration. + This second file has to be next to the FITS file. + + Args: + path (str): path of the FITS file + + Returns: + WaveformsContainer: WaveformsContainer instance + """ + log.info(f"loading from {path}") + hdul = fits.open(Path(path)) + + cls = WaveformsContainer.__new__(WaveformsContainer) + + cls.__run_number = hdul[0].header['RUN'] + cls.__nevents = hdul[0].header['NEVENTS'] + cls.__npixels = hdul[0].header['NPIXELS'] + cls.__nsamples = hdul[0].header['NSAMPLES'] + + cls.__subarray = SubarrayDescription.from_hdf(Path(path.replace('waveforms_','subarray_').replace('fits','hdf5'))) + + + cls.__pixels_id = hdul[0].data + cls.wfs_hg = hdul[1].data + cls.wfs_lg = hdul[2].data + + table_prop = hdul[3].data + cls.event_id = table_prop["event_id"] + cls.event_type = table_prop["event_type"] + cls.ucts_timestamp = table_prop["ucts_timestamp"] + cls.ucts_busy_counter = table_prop["ucts_busy_counter"] + cls.ucts_event_counter = table_prop["ucts_event_counter"] + + table_trigger = hdul[4].data + cls.trig_pattern_all = table_trigger["trig_pattern_all"] + + return cls + + + + + + + def compute_trigger_patern(self) : + #mean.shape nevents * npixels + mean,std =np.mean(self.wfs_hg,axis = 2),np.std(self.wfs_hg,axis = 2) + self.trig_pattern = self.wfs_hg.max(axis = 2) > (mean + 3 * std) + + + + ##methods used to display + def display(self,evt,cmap = 'gnuplot2') : + image = self.wfs_hg.sum(axis=2) + disp = CameraDisplay(geometry=WaveformsContainer.CAMERA, image=image[evt], cmap=cmap) + disp.add_colorbar() + return disp + + def plot_waveform(self,evt) : + fig,ax = plt.subplots(1,1) + ax.plot(self.wfs_hg[evt].T) + return fig,ax + + + + @property + def reader(self) : return self.__reader + + @property + def npixels(self) : return self.__npixels + + @property + def nsamples(self) : return self.__nsamples + + @property + def geometry(self) : return self.__geometry + + @property + def subarray(self) : return self.__subarray + + @property + def pixels_id(self) : return self.__pixels_id + + @property + def nevents(self) : return self.__nevents + + @property + def run_number(self) : return self.__run_number + + + #physical properties + @property + def multiplicity(self) : return np.uint16(np.count_nonzero(self.trig_pattern,axis = 1)) + + @property + def trig_pattern(self) : return self.trig_pattern_all.any(axis = 1) + + + + + + + +class WaveformsContainers() : + """This class is to be used for computing waveforms of a run treating run files one by one + """ + + def __init__(self,run_number : int,max_events : int = None) : + log.info('Initialization of WaveformsContainers') + _,filenames = DataManagment.findrun(run_number) + self.waveformsContainer = [] + self.__nWaveformsContainer = 0 + for i,file in enumerate(filenames) : + self.waveformsContainer.append(WaveformsContainer(run_number,max_events=max_events,run_file=file)) + self.__nWaveformsContainer += 1 + if not(max_events is None) : max_events -= self.waveformsContainer[i].nevents + log.info(f'WaveformsContainer number {i} is created') + if max_events <= 0 : break + + def load_wfs(self,compute_trigger_patern = False) : + for i in range(self.__nWaveformsContainer) : + self.waveformsContainer[i].load_wfs(compute_trigger_patern = compute_trigger_patern) + + def write(self,path : str, **kwargs) : + for i in range(self.__nWaveformsContainer) : + self.waveformsContainer[i].write(path,suffix = f"{i:04d}" ,**kwargs) + + + @property + def nWaveformsContainer(self) : return self.__nWaveformsContainer + + @property + def nevents(self) : + return np.sum([self.waveformsContainer[i].nevents for i in range(self.__nWaveformsContainer)]) + + def append(self,waveformsContainer : WaveformsContainer) : + self.waveforms.append(waveformsContainer) + self.__nWaveformsContainer += 1 diff --git a/nectarchain/user_scripts/ggrolleron/load_wfs_compute_charge.py b/nectarchain/user_scripts/ggrolleron/load_wfs_compute_charge.py new file mode 100644 index 00000000..54cfd2c3 --- /dev/null +++ b/nectarchain/user_scripts/ggrolleron/load_wfs_compute_charge.py @@ -0,0 +1,241 @@ +import numpy as np +#import pandas as pd +import matplotlib as mpl +import matplotlib.pyplot as plt +import matplotlib.cm as cm +from pathlib import Path +import sys +import os +import argparse +import json + +import logging +logging.basicConfig(format='%(asctime)s %(name)s %(levelname)s %(message)s',level=logging.INFO,filename = f"{os.environ.get('NECTARCHAIN_LOG')}/{Path(__file__).stem}_{os.getpid()}.log") +log = logging.getLogger(__name__) +##tips to add message to stdout +handler = logging.StreamHandler(sys.stdout) +handler.setLevel(logging.DEBUG) +formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s') +handler.setFormatter(formatter) +log.addHandler(handler) + +from nectarchain.calibration.container import WaveformsContainer, WaveformsContainers +from nectarchain.calibration.container import ChargeContainer, ChargeContainers + +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=list) + +#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=list) + +#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' + ) + +#extractor arguments +parser.add_argument('--extractorMethod', + choices=["FullWaveformSum","LocalPeakWindowSum"], + 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 + ) + +args = parser.parse_args() + +#control shape of arguments lists +for arg in ['spe','ff','ped'] : + run_number = eval(f"args.{arg}_run_number") + max_events = eval(f"args.{arg}_max_events") + nevents = eval(f"args.{arg}_nevents") + + if not(max_events is None) and len(max_events) != len(run_number) : + e = Exception(f'{arg}_run_number and {arg}_max_events must have same length') + log.error(e,exc_info=True) + raise e + if not(nevents is None) and len(nevents) != len(run_number) : + e = Exception(f'{arg}_run_number and {arg}_nevents must have same length') + log.error(e,exc_info=True) + raise e + + +def load_wfs_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: + runs_list (list): list of runs for which you want to perfrom waveforms and charge extraction + reload_wfs (bool, optional): argument used to reload waveforms from pre-loaded waveforms (in fits format) or from zfits file. Defaults to False. + overwrite (bool, optional): to overwrite file on disk. Defaults to False. + charge_extraction_method (str, optional): ctapipe charge extractor. Defaults to "FullWaveformSum". + + Raises: + e : an error occurred during zfits loading from ctapipe EventSource + """ + + #print(runs_list) + #print(charge_extraction_method) + #print(overwrite) + #print(reload_wfs) + #print(kwargs) + + + max_events = kwargs.get("max_events",[None for i in range(len(runs_list))]) + nevents = kwargs.get("nevents",[-1 for i in range(len(runs_list))]) + + charge_childpath = kwargs.get("charge_childpath",charge_extraction_method) + extractor_kwargs = kwargs.get("extractor_kwargs",{}) + + + 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 : + 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 sps/hess/lpnhe/ggroller/projects/NECTARCAM/runs/waveforms/waveforms_run{runs_list[i]}.fits") + log.warning(f"reloading from zfits files") + 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) + except Exception as e : + log.error(e,exc_info = True) + raise e + else : + 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) + + log.info(f"computation of charge with {charge_childpath}") + charge = ChargeContainer.from_waveforms(wfs,method = charge_childpath,**extractor_kwargs) + charge.write(f"{os.environ['NECTARCAMDATA']}/charges/{path}/",overwrite = overwrite) + del wfs,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('spe_nevents',[-1 for i in range(len(ff_run_number))]) + ped_nevents = kwargs.pop('spe_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('spe_max_events',[None for i in range(len(ff_run_number))]) + ped_max_events = kwargs.pop('spe_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") + + + + 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] + + + arg = vars(args) + log.info(f"arguments are : {arg}") + + key_to_pop = [] + for key in arg.keys() : + if arg[key] is None : + key_to_pop.append(key) + + for key in key_to_pop : + arg.pop(key) + + log.info(f"arguments passed to main are : {arg}") + + path= args.extractorMethod+'_4-12' + arg['path'] = path + + main(**arg) + + diff --git a/nectarchain/user_scripts/ggrolleron/load_wfs_compute_charge.sh b/nectarchain/user_scripts/ggrolleron/load_wfs_compute_charge.sh new file mode 100644 index 00000000..89165512 --- /dev/null +++ b/nectarchain/user_scripts/ggrolleron/load_wfs_compute_charge.sh @@ -0,0 +1,7 @@ +#HOW TO : +#spe_run_number = [2633,2634,3784] +#ff_run_number = [2608] +#ped_run_number = [2609] +#spe_nevents = [49227,49148,-1] + +python load_wfs_compute_charge.py -s 2633 2634 3784 -f 2608 -p 2609 --spe_nevents 49227 49148 -1 --extractorMethod LocalPeakWindowSum --extractor_kwargs '{"window_width":16,"window_shift":4}' \ No newline at end of file diff --git a/nectarchain/user_scripts/ggrolleron/test.py b/nectarchain/user_scripts/ggrolleron/test.py new file mode 100644 index 00000000..bdf8125d --- /dev/null +++ b/nectarchain/user_scripts/ggrolleron/test.py @@ -0,0 +1,48 @@ +import numpy as np +#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) +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') +handler.setFormatter(formatter) +log.addHandler(handler) + + + +from nectarchain.calibration.container import ChargeContainer,WaveformsContainer +from nectarchain.calibration.container.utils import DataManagment + + + + +run_number = [2633] +ped_run_number = [2630] +FF_run_number = [2609] + +spe_run_1000V = WaveformsContainer(run_number[0],max_events = 1000) + +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.compute_charge(spe_run_1000V,1,method = "gradient_extractor") + + + +charge = ChargeContainer.from_waveform(spe_run_1000V) + + +charge.write(f"{os.environ['NECTARCAMDATA']}/charges/std/",overwrite = True) + + +print("work completed") \ No newline at end of file