diff --git a/src/nectarchain/data/container/pedestal_container.py b/src/nectarchain/data/container/pedestal_container.py index b0c11d72..65fa6991 100644 --- a/src/nectarchain/data/container/pedestal_container.py +++ b/src/nectarchain/data/container/pedestal_container.py @@ -50,7 +50,14 @@ class NectarCAMPedestalContainer(NectarCAMContainer): pedestal_std_hg (np.ndarray): An array of standard deviations of high gain pedestals. pedestal_std_lg (np.ndarray): An array of standard deviations of low gain - pedestals. + pedestal_charge_mean_hg (np.ndarray): An array of high gain + mean pedestal charges. + pedestal_charge_mean_lg (np.ndarray): An array of low gain + mean pedestal charges. + pedestal_charge_std_hg (np.ndarray): An array of standard deviations + of high gain pedestal charges. + pedestal_charge_std_lg (np.ndarray): An array of standard deviations + of low gain pedestal charges. """ nsamples = Field( @@ -101,6 +108,33 @@ class NectarCAMPedestalContainer(NectarCAMContainer): description="low gain pedestals standard deviations", ) + pedestal_charge_mean_hg = Field( + type=np.ndarray, + dtype=np.float64, + ndim=1, + description="high gain mean pedestal charges", + ) + + pedestal_charge_mean_lg = Field( + type=np.ndarray, + dtype=np.float64, + ndim=1, + description="low gain mean pedestal charges", + ) + + pedestal_charge_std_hg = Field( + type=np.ndarray, + dtype=np.float64, + ndim=1, + description="high gain pedestal charges standard deviations", + ) + pedestal_charge_std_lg = Field( + type=np.ndarray, + dtype=np.float64, + ndim=1, + description="low gain pedestal charges standard deviations", + ) + pixel_mask = Field( type=np.ndarray, dtype=np.int8, diff --git a/src/nectarchain/data/container/tests/test_pedestal.py b/src/nectarchain/data/container/tests/test_pedestal.py index 07ba89b6..fba63c0e 100644 --- a/src/nectarchain/data/container/tests/test_pedestal.py +++ b/src/nectarchain/data/container/tests/test_pedestal.py @@ -19,6 +19,10 @@ def generate_mock_pedestal_container(): pedestal_mean_lg = np.float64(np.random.uniform(240, 260, size=(npixels, nsamples))) pedestal_std_hg = np.float64(np.random.normal(size=(npixels, nsamples))) pedestal_std_lg = np.float64(np.random.normal(size=(npixels, nsamples))) + pedestal_charge_mean_hg = np.float64(np.random.uniform(1.0e4, 1.5e4, size=npixels)) + pedestal_charge_mean_lg = np.float64(np.random.uniform(1.0e4, 1.5e4, size=npixels)) + pedestal_charge_std_hg = 30 * np.float64(np.random.normal(size=npixels)) + pedestal_charge_std_lg = 30 * np.float64(np.random.normal(size=npixels)) pixel_mask = np.int8(np.random.randint(0, 2, size=(nchannels, npixels))) # create pedestal container @@ -32,6 +36,10 @@ def generate_mock_pedestal_container(): pedestal_mean_lg=pedestal_mean_lg, pedestal_std_hg=pedestal_std_hg, pedestal_std_lg=pedestal_std_lg, + pedestal_charge_mean_hg=pedestal_charge_mean_hg, + pedestal_charge_mean_lg=pedestal_charge_mean_lg, + pedestal_charge_std_hg=pedestal_charge_std_hg, + pedestal_charge_std_lg=pedestal_charge_std_lg, pixel_mask=pixel_mask, ) pedestal_container.validate() @@ -47,6 +55,10 @@ def generate_mock_pedestal_container(): "pedestal_mean_lg": pedestal_mean_lg, "pedestal_std_hg": pedestal_std_hg, "pedestal_std_lg": pedestal_std_lg, + "pedestal_charge_mean_hg": pedestal_charge_mean_hg, + "pedestal_charge_mean_lg": pedestal_charge_mean_lg, + "pedestal_charge_std_hg": pedestal_charge_std_hg, + "pedestal_charge_std_lg": pedestal_charge_std_lg, "pixel_mask": pixel_mask, } @@ -72,6 +84,18 @@ def test_create_pedestal_container_mem(self): ) assert np.allclose(pedestal_container.pedestal_std_hg, dict["pedestal_std_hg"]) assert np.allclose(pedestal_container.pedestal_std_lg, dict["pedestal_std_lg"]) + assert np.allclose( + pedestal_container.pedestal_charge_mean_hg, dict["pedestal_charge_mean_hg"] + ) + assert np.allclose( + pedestal_container.pedestal_charge_mean_lg, dict["pedestal_charge_mean_lg"] + ) + assert np.allclose( + pedestal_container.pedestal_charge_std_hg, dict["pedestal_charge_std_hg"] + ) + assert np.allclose( + pedestal_container.pedestal_charge_std_lg, dict["pedestal_charge_std_lg"] + ) assert np.allclose(pedestal_container.pixel_mask, dict["pixel_mask"]) # FIXME diff --git a/src/nectarchain/makers/calibration/pedestal_makers.py b/src/nectarchain/makers/calibration/pedestal_makers.py index 73dd8d93..7ae4929d 100644 --- a/src/nectarchain/makers/calibration/pedestal_makers.py +++ b/src/nectarchain/makers/calibration/pedestal_makers.py @@ -15,7 +15,6 @@ log = logging.getLogger(__name__) log.handlers = logging.getLogger("__main__").handlers - __all__ = ["PedestalNectarCAMCalibrationTool"] @@ -31,6 +30,63 @@ class PedestalNectarCAMCalibrationTool(NectarCAMCalibrationTool): def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) + @staticmethod + def mean_std_multisample(nevents, means, stds): + """ + Method that calculates means and std of the combination of multiple subsamples. + Works for both: + + - pedestal data (means/stds shaped (n_pixels, n_samples)) + - charge data (means/stds shaped (n_pixels,)) + + Parameters + ---------- + nevents : list of `~numpy.ndarray` + Number of events for each sample (per pixel) + means : list of `~numpy.ndarray` + Mean values + stds : list of `~numpy.ndarray` + Std values + + Returns + ------- + mean : `~numpy.ndarray` + Mean values of combined sample + std : `~numpy.ndarray` + Std values of combined sample + nevent : `~numpy.ndarray` + Number of events of combined sample (per pixel) + """ + + # convert lists to numpy arrays + # axis 0 corresponds to the subsamples + nevents = np.array(nevents) + means = np.array(means) + stds = np.array(stds) + + total_nevents = np.sum(nevents, axis=0) + + # Handle both 1D and 2D cases cleanly + if means.ndim == 3: + # (n_subsamples, n_pixels, n_samples) + nevents_expanded = nevents[:, :, np.newaxis] + total_nevents_expanded = total_nevents[:, np.newaxis] + elif means.ndim == 2: + # (n_subsamples, n_pixels) + nevents_expanded = nevents + total_nevents_expanded = total_nevents + else: + log.error("Unexpected shape for means array") + + mean = np.sum(nevents_expanded * means, axis=0) / total_nevents_expanded + num = np.sum( + (nevents_expanded - 1) * stds**2 + nevents_expanded * (means - mean) ** 2, + axis=0, + ) + std = np.sqrt(num / (total_nevents_expanded - 1)) + + return mean, std, total_nevents + def _init_output_path(self): """ Initialize output path @@ -76,96 +132,93 @@ def _combine_results(self): # Loop over sliced results to fill the combined results self.log.info("Combine sliced results") - for i, _pedestalContainer in enumerate(pedestalContainers): + nevents_list = [] + mean_lists = {"ped_hg": [], "ped_lg": [], "charge_hg": [], "charge_lg": []} + std_lists = {"ped_hg": [], "ped_lg": [], "charge_hg": [], "charge_lg": []} + first = True + for _pedestalContainer in pedestalContainers: pedestalContainer = list(_pedestalContainer.containers.values())[0] - if i == 0: - # initialize fields for the combined results based on first slice + + # usable pixel mask + usable_pixels = pedestalContainer.pixel_mask == 0 + usable_pixels = np.logical_and(usable_pixels[0], usable_pixels[1]) + + nevents_list.append(pedestalContainer.nevents * usable_pixels) + mean_lists["ped_hg"].append(pedestalContainer.pedestal_mean_hg) + mean_lists["ped_lg"].append(pedestalContainer.pedestal_mean_lg) + std_lists["ped_hg"].append(pedestalContainer.pedestal_std_hg) + std_lists["ped_lg"].append(pedestalContainer.pedestal_std_lg) + mean_lists["charge_hg"].append(pedestalContainer.pedestal_charge_mean_hg) + mean_lists["charge_lg"].append(pedestalContainer.pedestal_charge_mean_lg) + std_lists["charge_hg"].append(pedestalContainer.pedestal_charge_std_hg) + std_lists["charge_lg"].append(pedestalContainer.pedestal_charge_std_lg) + + if first: nsamples = pedestalContainer.nsamples - nevents = np.zeros(len(pedestalContainer.nevents)) pixels_id = pedestalContainer.pixels_id ucts_timestamp_min = pedestalContainer.ucts_timestamp_min ucts_timestamp_max = pedestalContainer.ucts_timestamp_max - pedestal_mean_hg = np.zeros( - np.shape(pedestalContainer.pedestal_mean_hg) - ) - pedestal_mean_lg = np.zeros( - np.shape(pedestalContainer.pedestal_mean_lg) - ) - pedestal_std_hg = np.zeros(np.shape(pedestalContainer.pedestal_std_hg)) - pedestal_std_lg = np.zeros(np.shape(pedestalContainer.pedestal_std_lg)) - else: - # otherwise consider the overall time interval - ucts_timestamp_min = np.minimum( - ucts_timestamp_min, pedestalContainer.ucts_timestamp_min - ) - ucts_timestamp_max = np.maximum( - ucts_timestamp_max, pedestalContainer.ucts_timestamp_max - ) - # for all slices - # derive from pixel mask a mask that sets usable pixels - # accept only pixels for which no flags were raised - usable_pixels = pedestalContainer.pixel_mask == 0 - # use a pixel only if it has no flag on either channel - usable_pixels = np.logical_and(usable_pixels[0], usable_pixels[1]) - # cumulated number of events - nevents += pedestalContainer.nevents * usable_pixels - # add mean, std sum elements - pedestal_mean_hg += ( - pedestalContainer.pedestal_mean_hg - * pedestalContainer.nevents[:, np.newaxis] - * usable_pixels[:, np.newaxis] - ) - pedestal_mean_lg += ( - pedestalContainer.pedestal_mean_lg - * pedestalContainer.nevents[:, np.newaxis] - * usable_pixels[:, np.newaxis] + first = False + + # update timestamps + ucts_timestamp_min = np.minimum( + ucts_timestamp_min, pedestalContainer.ucts_timestamp_min ) - pedestal_std_hg += ( - pedestalContainer.pedestal_std_hg**2 - * pedestalContainer.nevents[:, np.newaxis] - * usable_pixels[:, np.newaxis] + ucts_timestamp_max = np.maximum( + ucts_timestamp_max, pedestalContainer.ucts_timestamp_max ) - pedestal_std_lg += ( - pedestalContainer.pedestal_std_lg**2 - * pedestalContainer.nevents[:, np.newaxis] - * usable_pixels[:, np.newaxis] + + # Compute combined stats for both gains + results = {} + for q in ["ped_hg", "ped_lg", "charge_hg", "charge_lg"]: + mean, std, nevents = self.mean_std_multisample( + nevents_list, mean_lists[q], std_lists[q] ) - # calculate final values of mean and std - pedestal_mean_hg /= nevents[:, np.newaxis] - pedestal_mean_lg /= nevents[:, np.newaxis] - pedestal_std_hg /= nevents[:, np.newaxis] - pedestal_std_hg = np.sqrt(pedestal_std_hg) - pedestal_std_lg /= nevents[:, np.newaxis] - pedestal_std_lg = np.sqrt(pedestal_std_lg) + results[q] = {"mean": mean, "std": std} + + mean_hg, std_hg = results["ped_hg"]["mean"], results["ped_hg"]["std"] + mean_lg, std_lg = results["ped_lg"]["mean"], results["ped_lg"]["std"] + charge_mean_hg, charge_std_hg = ( + results["charge_hg"]["mean"], + results["charge_hg"]["std"], + ) + charge_mean_lg, charge_std_lg = ( + results["charge_lg"]["mean"], + results["charge_lg"]["std"], + ) + # flag bad pixels in overall results based on same criteria as for individual - # slides - # reconstitute dictionary with cumulated results consistently with - # PedestalComponent ped_stats = {} - array_shape = np.append([N_GAINS], np.shape(pedestal_mean_hg)) + array_shape = np.append([N_GAINS], np.shape(mean_hg)) for statistic in ["mean", "std"]: ped_stat = np.zeros(array_shape) if statistic == "mean": - ped_stat[HIGH_GAIN] = pedestal_mean_hg - ped_stat[LOW_GAIN] = pedestal_mean_lg + ped_stat[HIGH_GAIN] = mean_hg + ped_stat[LOW_GAIN] = mean_lg elif statistic == "std": - ped_stat[HIGH_GAIN] = pedestal_std_hg - ped_stat[LOW_GAIN] = pedestal_std_lg + ped_stat[HIGH_GAIN] = std_hg + ped_stat[LOW_GAIN] = std_lg # Store the result in the dictionary ped_stats[statistic] = ped_stat # use flagging method from PedestalComponent pixel_mask = self.components[0].flag_bad_pixels(ped_stats, nevents) + # reconstitute dictionary with cumulated results consistently with + # PedestalComponent output = NectarCAMPedestalContainer( nsamples=nsamples, nevents=nevents, pixels_id=pixels_id, ucts_timestamp_min=ucts_timestamp_min, ucts_timestamp_max=ucts_timestamp_max, - pedestal_mean_hg=pedestal_mean_hg, - pedestal_mean_lg=pedestal_mean_lg, - pedestal_std_hg=pedestal_std_hg, - pedestal_std_lg=pedestal_std_lg, + pedestal_mean_hg=mean_hg, + pedestal_mean_lg=mean_lg, + pedestal_std_hg=std_hg, + pedestal_std_lg=std_lg, + pedestal_charge_mean_hg=charge_mean_hg, + pedestal_charge_mean_lg=charge_mean_lg, + pedestal_charge_std_hg=charge_std_hg, + pedestal_charge_std_lg=charge_std_lg, pixel_mask=pixel_mask, ) diff --git a/src/nectarchain/makers/calibration/tests/test_pedestal_tool.py b/src/nectarchain/makers/calibration/tests/test_pedestal_tool.py index ffcbbdf5..e5ceb012 100644 --- a/src/nectarchain/makers/calibration/tests/test_pedestal_tool.py +++ b/src/nectarchain/makers/calibration/tests/test_pedestal_tool.py @@ -81,7 +81,11 @@ def test_base(self): assert np.allclose(output.pedestal_mean_hg, 245.0, atol=20.0) assert np.allclose(output.pedestal_mean_lg, 245.0, atol=20.0) assert np.allclose(output.pedestal_std_hg, 10, atol=10) - assert np.allclose(output.pedestal_std_lg, 2.5, atol=2.3) + assert np.allclose(output.pedestal_std_lg, 3.0, atol=2.5) + assert np.allclose(output.pedestal_charge_mean_hg, 14700.0, rtol=0.1) + assert np.allclose(output.pedestal_charge_mean_lg, 14700.0, rtol=0.1) + assert np.allclose(output.pedestal_charge_std_hg, 50.0, atol=50) + assert np.allclose(output.pedestal_charge_std_lg, 50.0, atol=50) # Check output on disk pedestalContainers = NectarCAMPedestalContainers.from_hdf5(outfile) @@ -158,8 +162,18 @@ def test_base(self): assert np.allclose(pedestalContainer.pedestal_mean_hg, 245.0, atol=20.0) assert np.allclose(pedestalContainer.pedestal_mean_lg, 245.0, atol=20.0) assert np.allclose(pedestalContainer.pedestal_std_hg, 10, atol=10) + assert np.allclose(pedestalContainer.pedestal_std_lg, 3, atol=2.5) assert np.allclose( - pedestalContainer.pedestal_std_lg, 2.5, atol=2.0 if i == 0 else 2.3 + pedestalContainer.pedestal_charge_mean_hg, 14700.0, rtol=0.1 + ) + assert np.allclose( + pedestalContainer.pedestal_charge_mean_lg, 14700.0, rtol=0.1 + ) + assert np.allclose( + pedestalContainer.pedestal_charge_std_hg, 50.0, atol=50 + ) + assert np.allclose( + pedestalContainer.pedestal_charge_std_lg, 50.0, atol=50 ) def test_timesel(self): @@ -215,9 +229,11 @@ def test_timesel(self): assert np.allclose(output.pedestal_mean_hg, 245.0, atol=20.0) assert np.allclose(output.pedestal_mean_lg, 245.0, atol=20.0) assert np.allclose(output.pedestal_std_hg, 10, atol=10) - assert np.allclose( - output.pedestal_std_lg, 2.5, atol=2.0 if i == 0 else 2.3 - ) + assert np.allclose(output.pedestal_std_lg, 3, atol=2.5) + assert np.allclose(output.pedestal_charge_mean_hg, 14700.0, rtol=0.1) + assert np.allclose(output.pedestal_charge_mean_lg, 14700.0, rtol=0.1) + assert np.allclose(output.pedestal_charge_std_hg, 50.0, atol=50) + assert np.allclose(output.pedestal_charge_std_lg, 50.0, atol=50) def test_WaveformsStdFilter(self): """ @@ -266,12 +282,13 @@ def test_WaveformsStdFilter(self): assert np.allclose(output.pedestal_mean_hg, 245.0, atol=20.0) assert np.allclose(output.pedestal_mean_lg, 245.0, atol=20.0) # verify that fluctuations are reduced - assert np.allclose( - output.pedestal_std_hg, 3.0, atol=2.0 if i == 0 else 4.0 - ) - assert np.allclose( - output.pedestal_std_lg, 2.5, atol=2.0 if i == 0 else 2.6 - ) + assert np.allclose(output.pedestal_std_hg, 4.0, atol=4.0) + assert np.allclose(output.pedestal_std_lg, 2.5, atol=3.0) + assert np.allclose(output.pedestal_charge_mean_hg, 14700.0, rtol=0.1) + assert np.allclose(output.pedestal_charge_mean_lg, 14700.0, rtol=0.1) + # verify that fluctuations are reduced + assert np.allclose(output.pedestal_charge_std_hg, 40.0, atol=50) + assert np.allclose(output.pedestal_charge_std_lg, 40.0, atol=50) def test_ChargeDistributionFilter(self): """ @@ -321,9 +338,11 @@ def test_ChargeDistributionFilter(self): assert np.allclose(output.pedestal_mean_hg, 245.0, atol=20.0) assert np.allclose(output.pedestal_mean_lg, 245.0, atol=20.0) assert np.allclose(output.pedestal_std_hg, 10.0, atol=10.0) - assert np.allclose( - output.pedestal_std_lg, 2.5 if i == 0 else 2.2, atol=3.0 - ) + assert np.allclose(output.pedestal_std_lg, 3, atol=2.6) + assert np.allclose(output.pedestal_charge_mean_hg, 14700.0, rtol=0.1) + assert np.allclose(output.pedestal_charge_mean_lg, 14700.0, rtol=0.1) + assert np.allclose(output.pedestal_charge_std_hg, 50.0, atol=50) + assert np.allclose(output.pedestal_charge_std_lg, 50.0, atol=50) def test_pixel_mask(self): """ diff --git a/src/nectarchain/makers/component/pedestal_component.py b/src/nectarchain/makers/component/pedestal_component.py index 819eeb36..de69c62c 100644 --- a/src/nectarchain/makers/component/pedestal_component.py +++ b/src/nectarchain/makers/component/pedestal_component.py @@ -202,38 +202,50 @@ def __init__(self, subarray, config=None, parent=None, *args, **kwargs): ) @staticmethod - def calculate_stats(waveformsContainers, wfs_mask, statistics): + def calculate_stats(container, mask, statistics): """Calculate statistics for the pedestals from a waveforms container. Parameters ---------- - waveformsContainers : `~nectarchain.data.container.WaveformsContainer` - Waveforms container - wfs_mask : `numpy.ndarray` - Mask to apply to exclude outliers with shape (n_pixels,n_samples) + containers : `~nectarchain.data.container` + Waveforms or Charges container + mask : `numpy.ndarray` + Mask to apply to exclude outliers with shape (n_events,n_pixels,n_samples) statistics : `list` Names of the statistics (numpy.ma attributes) to compute Returns ------- ped_stats : `dict` - A dictionary containing 3D (n_chan,n_pixels,n_samples) arrays for each - statistic. + A dictionary containing 3D (n_chan,n_pixels,n_samples) + or 2D (n_chan,n_pixels) arrays for each statistic. """ + # Detect container type by attribute name + if hasattr(container, "wfs_hg"): + data_hg = container.wfs_hg + data_lg = container.wfs_lg + is_waveform = True + else: + data_hg = container.charges_hg + data_lg = container.charges_lg + is_waveform = False + + # Adapt mask shape + # Waveforms: (n_events, n_pixels, n_samples) + # Charges: (n_events, n_pixels) + if not is_waveform: + mask = mask[:, :, 0] + ped_stats = {} for stat in statistics: # Calculate the statistic along axis = 0, that is over events - ped_stat_hg = getattr(ma, stat)( - ma.masked_array(waveformsContainers.wfs_hg, wfs_mask), axis=0 - ) - ped_stat_lg = getattr(ma, stat)( - ma.masked_array(waveformsContainers.wfs_lg, wfs_mask), axis=0 - ) + ped_stat_hg = getattr(ma, stat)(ma.masked_array(data_hg, mask), axis=0) + ped_stat_lg = getattr(ma, stat)(ma.masked_array(data_lg, mask), axis=0) # Create a 3D array for the statistic - array_shape = np.append([N_GAINS], np.shape(waveformsContainers.wfs_hg[0])) + array_shape = np.append([N_GAINS], np.shape(data_hg[0])) ped_stat = np.zeros(array_shape) ped_stat[HIGH_GAIN] = ped_stat_hg ped_stat[LOW_GAIN] = ped_stat_lg @@ -306,7 +318,7 @@ def flag_bad_pixels(self, ped_stats, nevents): # Flag on standard deviation per pixel # Standard deviation of pedestal in channel/pixel above threshold log.info( - f"Flag pixels with pedestal standard deviation in a chennel/pixel above " + f"Flag pixels with pedestal standard deviation in a channel/pixel above " f"the maximum acceptable value {self.pixel_mask_std_pixel_max}" ) flag_pixel_std = np.int8( @@ -496,6 +508,8 @@ def finish(self, *args, **kwargs): self._waveformsContainers = waveformsContainers.containers[ EventType.SKY_PEDESTAL ] + # log.info('JPL: waveformsContainers=',waveformsContainers.containers[ + # EventType.SKY_PEDESTAL].nsamples) # Check if waveforms container is empty if self._waveformsContainers is None: @@ -510,12 +524,8 @@ def finish(self, *args, **kwargs): # container with no results return None else: - # If we want to filter based on charges distribution - # make sure that the charge distribution container is filled - if ( - self.filter_method == "ChargeDistributionFilter" - and self._chargesContainers is None - ): + # Make sure that the charge distribution container is filled + if self._chargesContainers is None: log.debug("Compute charges from waveforms") chargesComponent_kwargs = {} chargesComponent_configurable_traits = ( @@ -575,6 +585,9 @@ def finish(self, *args, **kwargs): self._ped_stats = self.calculate_stats( self._waveformsContainers, self._wfs_mask, statistics ) + self._charge_stats = self.calculate_stats( + self._chargesContainers, self._wfs_mask, statistics + ) # calculate the number of events per pixel used to compute the quantitites # start wit total number of events @@ -599,6 +612,10 @@ def finish(self, *args, **kwargs): pedestal_mean_lg=self._ped_stats["mean"][LOW_GAIN], pedestal_std_hg=self._ped_stats["std"][HIGH_GAIN], pedestal_std_lg=self._ped_stats["std"][LOW_GAIN], + pedestal_charge_mean_hg=self._charge_stats["mean"][HIGH_GAIN], + pedestal_charge_mean_lg=self._charge_stats["mean"][LOW_GAIN], + pedestal_charge_std_hg=self._charge_stats["std"][HIGH_GAIN], + pedestal_charge_std_lg=self._charge_stats["std"][LOW_GAIN], pixel_mask=pixel_mask, ) diff --git a/src/nectarchain/user_scripts/ltibaldo/extract_pedestals.py b/src/nectarchain/user_scripts/ltibaldo/extract_pedestals.py new file mode 100644 index 00000000..238761c8 --- /dev/null +++ b/src/nectarchain/user_scripts/ltibaldo/extract_pedestals.py @@ -0,0 +1,70 @@ +import logging +import os +import pathlib +from multiprocessing import Pool + +import numpy as np + +logging.basicConfig( + format="%(asctime)s %(name)s %(levelname)s %(message)s", level=logging.INFO +) +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + +from nectarchain.makers.calibration import PedestalNectarCAMCalibrationTool + +events_per_slice = 500 +nthreads = 30 + +run_list = np.concatenate( + ( + # np.arange(6882, 6891), + # np.arange(6672, 6681), + # np.arange(7144, 7153), + # np.arange(6954, 6963), + # np.arange(7020, 7029), + # np.arange(6543, 6552), + # np.arange(7077, 7086), + # np.arange(7153, 7181), + # np.arange(7182, 7190), + # np.arange(6891, 6927), + # np.arange(6681, 6717), + # np.arange(6552, 6588), + # np.arange(6963, 6999), + # np.arange(7086, 7122), + # np.arange(7029, 7065), + np.arange(6552, 6588), + np.arange(7086, 7110), + ) +) + + +def process_run(run_number): + outfile = os.environ["NECTARCAMDATA"] + "/runs/pedestal_cfilt3s_{}.h5".format( + run_number + ) + tool = PedestalNectarCAMCalibrationTool( + progress_bar=True, + run_number=run_number, + max_events=1999, + events_per_slice=events_per_slice, + log_level=0, + output_path=outfile, + overwrite=True, + filter_method="ChargeDistributionFilter", + charge_sigma_low_thr=3.0, + charge_sigma_high_thr=3.0, + ) + + tool.initialize() + tool.setup() + + tool.start() + tool.finish() + + +args = [int(x) for x in run_list] +pool = Pool(processes=nthreads) +pool.map(process_run, args) +pool.close() +pool.join() diff --git a/src/nectarchain/user_scripts/ltibaldo/extract_waveforms.py b/src/nectarchain/user_scripts/ltibaldo/extract_waveforms.py index 9767b39b..2e2cf33d 100644 --- a/src/nectarchain/user_scripts/ltibaldo/extract_waveforms.py +++ b/src/nectarchain/user_scripts/ltibaldo/extract_waveforms.py @@ -8,17 +8,17 @@ from nectarchain.makers import WaveformsNectarCAMCalibrationTool -run_number = 3938 +run_number = 6954 tool = WaveformsNectarCAMCalibrationTool( progress_bar=True, run_number=run_number, max_events=499, log_level=20, - events_per_slice = 100, + events_per_slice=100, ) tool.initialize() tool.setup() tool.start() -output = tool.finish() \ No newline at end of file +output = tool.finish() diff --git a/src/nectarchain/user_scripts/ltibaldo/pedestal_NSB.py b/src/nectarchain/user_scripts/ltibaldo/pedestal_NSB.py new file mode 100644 index 00000000..e7d5d08c --- /dev/null +++ b/src/nectarchain/user_scripts/ltibaldo/pedestal_NSB.py @@ -0,0 +1,399 @@ +import os + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import seaborn as sns +import tables + +# Imports from ctapipe +from ctapipe.coordinates import EngineeringCameraFrame +from ctapipe.instrument import CameraGeometry +from ctapipe.visualization import CameraDisplay +from ctapipe_io_nectarcam.constants import N_PIXELS +from scipy import stats + +sns.set(style="whitegrid") + +rms_threshold = 47 + +pedfiles = [ # OK + { + "T": -5, + "data": [ + {"NSB": 0, "runs": np.arange(6882, 6891)}, + {"NSB": 10.6, "runs": np.arange(6891, 6900)}, + {"NSB": 20.4, "runs": np.arange(6900, 6909)}, + {"NSB": 39.8, "runs": np.arange(6909, 6918)}, + {"NSB": 78.8, "runs": np.arange(6918, 6927)}, + ], + }, + { + "T": 0, + "data": [ # OK + {"NSB": 0, "runs": np.arange(6672, 6681)}, + {"NSB": 10.6, "runs": np.arange(6681, 6690)}, + {"NSB": 20.4, "runs": np.arange(6690, 6699)}, + {"NSB": 39.8, "runs": np.arange(6699, 6708)}, + {"NSB": 78.8, "runs": np.arange(6708, 6717)}, + ], + }, + { + "T": 5, + "data": [ # OK + {"NSB": 0, "runs": np.arange(6543, 6552)}, + {"NSB": 10.6, "runs": np.arange(6552, 6561)}, + {"NSB": 20.4, "runs": np.arange(6661, 6570)}, + {"NSB": 39.8, "runs": np.arange(6570, 6579)}, + {"NSB": 78.8, "runs": np.arange(6579, 6588)}, + ], + }, + { + "T": 10, + "data": [ # OK + {"NSB": 0, "runs": np.arange(7144, 7153)}, + {"NSB": 10.6, "runs": np.arange(7153, 7162)}, + {"NSB": 20.4, "runs": np.arange(7162, 7171)}, + {"NSB": 39.8, "runs": np.arange(7171, 7180)}, + {"NSB": 78.8, "runs": np.append([7180], np.arange(7182, 7190))}, + ], + }, + { + "T": 14, + "data": [ # OK + {"NSB": 0, "runs": np.arange(6954, 6963)}, + {"NSB": 10.6, "runs": np.arange(6963, 6972)}, + {"NSB": 20.4, "runs": np.arange(6972, 6981)}, + {"NSB": 39.8, "runs": np.arange(6981, 6990)}, + {"NSB": 78.8, "runs": np.arange(6990, 6999)}, + ], + }, + { + "T": 20, + "data": [ # OK + {"NSB": 0, "runs": np.arange(7077, 7086)}, + {"NSB": 10.6, "runs": np.arange(7086, 7095)}, + {"NSB": 20.4, "runs": np.arange(7095, 7104)}, + {"NSB": 39.8, "runs": np.arange(7104, 7112)}, + {"NSB": 78.8, "runs": np.arange(7113, 7122)}, + ], + }, + { + "T": 25, + "data": [ # OK + {"NSB": 0, "runs": np.arange(7020, 7029)}, + {"NSB": 10.6, "runs": np.arange(7029, 7038)}, + {"NSB": 20.4, "runs": np.arange(7038, 7047)}, + {"NSB": 39.8, "runs": np.arange(7047, 7056)}, + {"NSB": 78.8, "runs": np.arange(7056, 7064)}, + ], + }, +] + +outfigroot = os.environ["NECTARCHAIN_FIGURES"] +pixel_display = [100, 144, 240, 723, 816, 1034, 1516] +fill_value = np.nan + +pixel_process = np.arange(N_PIXELS) +# pixel_process = pixel_display + +slopes_width_hg = np.full([len(pedfiles), N_PIXELS], fill_value=fill_value) +slopes_width_lg = np.full([len(pedfiles), N_PIXELS], fill_value=fill_value) + +for itemp, uberset in enumerate(pedfiles): + temp = uberset["T"] + print(f"Analyzing data for temperature {temp}") + + # create output directory if needed + outdir = os.path.join(outfigroot, f"NSB_T{temp}deg") + if not os.path.exists(outdir): + os.makedirs(outdir) + + slopes_ped_hg = np.full(N_PIXELS, fill_value=fill_value) + slopes_rms_hg = np.full(N_PIXELS, fill_value=fill_value) + slopes_ped_lg = np.full(N_PIXELS, fill_value=fill_value) + slopes_rms_lg = np.full(N_PIXELS, fill_value=fill_value) + + removed_pixels = 0 + flagged_pixels = 0 + + for pixel_id in pixel_process: + print(f"Working on pixel {pixel_id}") + # fill panda dataframe + nsbs = [] + channels = [] + avgpeds = [] + pedsrms = [] + + for dataset in uberset["data"]: + for run_number in dataset["runs"]: + run_number = int(run_number) + filename = ( + os.environ["NECTARCAMDATA"] + + f"/runs/pedestal_cfilt3s_{run_number}.h5" + ) + h5file = tables.open_file(filename) + table = h5file.root["data_combined"]["NectarCAMPedestalContainer_0"][0] + for channel in ["hg", "lg"]: + pedestal = table[f"pedestal_mean_{channel}"][ + table["pixels_id"] == pixel_id + ][0] + rms = table[f"pedestal_charge_std_{channel}"][ + table["pixels_id"] == pixel_id + ][0] + avgped = np.average(pedestal) + combrms = rms + nsbs.append(dataset["NSB"]) + channels.append(channel) + avgpeds.append(avgped) + pedsrms.append(combrms) + h5file.close() + + if np.any(np.isnan(avgpeds)) or np.any(np.isnan(combrms)): + flagged_pixels += 1 + print("Bad pixel") + elif np.any( + np.array(pedsrms)[np.array(nsbs) == 0.0] > rms_threshold + ): # filter on NSB OFF + removed_pixels += 1 + print("Removed pixel") + else: + df = pd.DataFrame( + { + "I_NSB (mA)": nsbs, + "channel": channels, + "pedestal (ADC)": avgpeds, + "pedestal width^2 (ADC^2)": np.power(pedsrms, 2), + } + ) + + for k, q in enumerate(["pedestal (ADC)", "pedestal width^2 (ADC^2)"]): + if pixel_id in pixel_display: + lm = sns.lmplot( + x="I_NSB (mA)", + y=q, + hue="channel", + data=df, + scatter_kws={"alpha": 0.0}, + legend=False, + ) + + ax = lm.axes[0][0] + fig = ax.get_figure() + fig.subplots_adjust(top=0.9) + # clean legend to avoid duplication + for _artist in ax.lines + ax.collections + ax.patches + ax.images: + _artist.set_label(s=None) + + sns.violinplot( + x="I_NSB (mA)", + y=q, + hue="channel", + split=True, + native_scale=True, + data=df, + ax=ax, + ) + + # extract slope and intercept for display and further usage + colors = sns.color_palette() + for s, channel in enumerate(["hg", "lg"]): + filtered_df = df[df["channel"] == channel] + slope, intercept, r_value, pv, se = stats.linregress( + filtered_df["I_NSB (mA)"], + filtered_df[q], + ) + if k == 0 and s == 0: + slopes_ped_hg[pixel_id] = slope + elif k == 0 and s == 1: + slopes_ped_lg[pixel_id] = slope + elif k == 1 and s == 0: + slopes_rms_hg[pixel_id] = slope + elif k == 1 and s == 1: + slopes_rms_lg[pixel_id] = slope + if pixel_id in pixel_display: + ax.annotate( + f"y = {slope:.4f} I + {intercept:.4f}", + (0.05, 0.85 - s * 0.05), + color=colors[s], + xycoords="axes fraction", + ) + + if pixel_id in pixel_display: + ax.set_title(f"Pixel {pixel_id}, T {temp} deg") + fig.savefig(f"{outdir}/pixel_{pixel_id}_{k}.png") + del fig + plt.close() + + print( + f"{flagged_pixels} pixels were flagged as bad by the pedestal estimation tool" + ) + print( + f"{removed_pixels} pixels were removed because they had a pedestal RMS " + f"exceeding {rms_threshold}" + ) + + slopes_width_hg[itemp] = slopes_rms_hg + slopes_width_lg[itemp] = slopes_rms_lg + + camera = CameraGeometry.from_name("NectarCam-003") + camera = camera.transform_to(EngineeringCameraFrame()) + + for k in range(2): # quantity to treat + fig = plt.figure(figsize=(8, 8)) + fig.subplots_adjust(wspace=0.35, hspace=0.3) + for s in range(2): # HG and LG channels + if k == 0 and s == 0: + slopes = slopes_ped_hg + elif k == 0 and s == 1: + slopes = slopes_ped_lg + elif k == 1 and s == 0: + slopes = slopes_rms_hg + elif k == 1 and s == 1: + slopes = slopes_rms_lg + # histograms + ax = plt.subplot(2, 2, s + 1) + plt.hist(slopes, bins="auto", histtype="step") + if k == 1: + ax.set_yscale("log") + if k == 0 and s == 0: + ax.set_title(f"Pedestal slope (ADC/mA), hg, T {temp} deg") + elif k == 0 and s == 1: + ax.set_title(f"Pedestal slope (ADC/mA), lg, T {temp} deg") + elif k == 1 and s == 0: + ax.set_title(f"Pedestal width slope (ADC^2/mA), hg, T {temp} deg") + elif k == 1 and s == 1: + ax.set_title(f"Pedestal width slope (ADC^2/mA), lg, T {temp} deg") + plt.axvline(np.nanmean(slopes), linestyle=":") + plt.axvline(np.nanmean(slopes) - np.nanstd(slopes), linestyle=":") + plt.axvline(np.nanmean(slopes) + np.nanstd(slopes), linestyle=":") + ax.annotate( + r"mean: {:.3f} $\pm$ {:.3f}".format( + np.nanmean(slopes), np.nanstd(slopes) + ), + (0.05, 0.85 - s * 0.05), + xycoords="axes fraction", + ) + # camera display + ax = plt.subplot(2, 2, s + 3) + disp = CameraDisplay(camera) + disp.image = slopes + # if k == 0: + # disp.set_limits_minmax(-0.3, 0.3) + # elif k == 1: + # disp.set_limits_minmax(-0.07, 0.07) + disp.cmap = "Spectral_r" + disp.add_colorbar() + disp.update() + ax.set_title("") + fig.savefig(f"{outdir}/camera_{k}.png") + +# Ped width^2-I slope as a function of temp +# create output directory if needed +outdir = os.path.join(outfigroot, "NSB_temperature") +if not os.path.exists(outdir): + os.makedirs(outdir) + +slopes_hg = np.full(N_PIXELS, fill_value=fill_value) +slopes_lg = np.full(N_PIXELS, fill_value=fill_value) +x = [data["T"] for data in pedfiles] + +print("Working on temperature variations") +bad_pixels_hg = 0 +bad_pixels_lg = 0 +for pixel_id in pixel_process: + print(f"Working on pixel {pixel_id}") + + for s in range(2): # channels + isgood = True + if s == 0: + y = slopes_width_hg[:, pixel_id] + label = "hg" + if np.any(np.isnan(y)): + isgood = False + bad_pixels_hg += 1 + elif s == 1: + y = slopes_width_lg[:, pixel_id] + label = "lg" + if np.any(np.isnan(y)): + isgood = False + bad_pixels_lg += 1 + + if isgood: + slope, intercept, r_value, pv, se = stats.linregress(x, y) + + if s == 0: + slopes_hg[pixel_id] = slope + elif s == 1: + slopes_lg[pixel_id] = slope + + if pixel_id in pixel_display: + df = pd.DataFrame( + { + "T (deg)": x, + "pedestal width^2-I slope (ADC^2/mA)": y, + } + ) + lm = sns.lmplot( + x="T (deg)", + y="pedestal width^2-I slope (ADC^2/mA)", + data=df, + scatter_kws={"alpha": 0.5}, + legend=False, + ) + + ax = lm.axes[0][0] + fig = ax.get_figure() + fig.subplots_adjust(top=0.9) + + ax.annotate( + f"y = {slope:.4f} T + {intercept:.4f}", + (0.05, 0.85 - s * 0.05), + xycoords="axes fraction", + ) + ax.set_title(f"Pixel {pixel_id}, {label}") + fig.savefig(f"{outdir}/pixel_{pixel_id}_{label}.png") + del fig + plt.close() + +print(f"Bad pixels: HG {bad_pixels_hg}, LG {bad_pixels_lg}") + +camera = CameraGeometry.from_name("NectarCam-003") +camera = camera.transform_to(EngineeringCameraFrame()) + +fig = plt.figure(figsize=(8, 8)) +fig.subplots_adjust(wspace=0.35, hspace=0.3) +for s in range(2): # HG and LG channels + if s == 0: + slopes = slopes_hg + elif s == 1: + slopes = slopes_lg + # histograms + ax = plt.subplot(2, 2, s + 1) + plt.hist(slopes, bins="auto", histtype="step") + if s == 0: + ax.set_title("Pedestal width slope (ADC^2/mA/deg), hg") + elif k == 1 and s == 1: + ax.set_title("Pedestal width slope (ADC^2/mA/deg), lg") + ax.set_yscale("log") + plt.axvline(np.nanmean(slopes), linestyle=":") + plt.axvline(np.nanmean(slopes) - np.nanstd(slopes), linestyle=":") + plt.axvline(np.nanmean(slopes) + np.nanstd(slopes), linestyle=":") + ax.annotate( + r"mean: {:.3f} $\pm$ {:.3f}".format(np.nanmean(slopes), np.nanstd(slopes)), + (0.05, 0.85 - s * 0.05), + xycoords="axes fraction", + ) + # camera display + ax = plt.subplot(2, 2, s + 3) + disp = CameraDisplay(camera) + disp.image = slopes + if s == 0: + disp.set_limits_minmax(-30.0, 10.0) + if s == 1: + disp.set_limits_minmax(-0.15, 0.15) + disp.cmap = "Spectral_r" + disp.add_colorbar() + disp.update() + ax.set_title("") +fig.savefig(f"{outdir}/camera.png") diff --git a/src/nectarchain/user_scripts/ltibaldo/pedestal_temperature.py b/src/nectarchain/user_scripts/ltibaldo/pedestal_temperature.py new file mode 100644 index 00000000..8c5a290c --- /dev/null +++ b/src/nectarchain/user_scripts/ltibaldo/pedestal_temperature.py @@ -0,0 +1,202 @@ +import os + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import seaborn as sns +import tables + +# Imports from ctapipe +from ctapipe.coordinates import EngineeringCameraFrame +from ctapipe.instrument import CameraGeometry +from ctapipe.visualization import CameraDisplay +from ctapipe_io_nectarcam.constants import N_PIXELS +from scipy import stats + +sns.set(style="whitegrid") + +rms_threshold = 47 + +data = [ + {"T": -5, "runs": np.arange(6882, 6891)}, + {"T": 0, "runs": np.arange(6672, 6681)}, + {"T": 5, "runs": np.arange(6543, 6552)}, + {"T": 10, "runs": np.arange(7144, 7153)}, + {"T": 14, "runs": np.arange(6954, 6963)}, + {"T": 20, "runs": np.arange(7077, 7086)}, + {"T": 25, "runs": np.arange(7020, 7029)}, +] + +pixel_display = [100, 144, 663, 1058, 491, 631, 656, 701, 756, 757, 1612, 1629] + +fill_value = np.nan +slopes_ped_hg = np.full(N_PIXELS, fill_value=fill_value) +slopes_rms_hg = np.full(N_PIXELS, fill_value=fill_value) +slopes_ped_lg = np.full(N_PIXELS, fill_value=fill_value) +slopes_rms_lg = np.full(N_PIXELS, fill_value=fill_value) + +removed_pixels = 0 +flagged_pixels = 0 + +for pixel_id in np.arange(N_PIXELS): + print(f"Working on pixel {pixel_id}") + # fill panda dataframe + temperatures = [] + channels = [] + avgpeds = [] + pedsrms = [] + + for dataset in data: + for run_number in dataset["runs"]: + run_number = int(run_number) + filename = ( + os.environ["NECTARCAMDATA"] + f"/runs/pedestal_cfilt3s_{run_number}.h5" + ) + h5file = tables.open_file(filename) + table = h5file.root["data_combined"]["NectarCAMPedestalContainer_0"][0] + for channel in ["hg", "lg"]: + pedestal = table[f"pedestal_mean_{channel}"][ + table["pixels_id"] == pixel_id + ][0] + rms = table[f"pedestal_charge_std_{channel}"][ + table["pixels_id"] == pixel_id + ][0] + avgped = np.average(pedestal) + combrms = rms + temperatures.append(dataset["T"]) + channels.append(channel) + avgpeds.append(avgped) + pedsrms.append(combrms) + h5file.close() + + if np.any(np.isnan(avgpeds)) or np.any(np.isnan(combrms)): + flagged_pixels += 1 + print("Bad pixel") + elif np.any(np.array(pedsrms) > rms_threshold): # filter on NSB OFF + removed_pixels += 1 + print("Removed pixel") + else: + df = pd.DataFrame( + { + "T (deg)": temperatures, + "channel": channels, + "pedestal (ADC)": avgpeds, + "pedestal width (ADC)": pedsrms, + } + ) + + for k, q in enumerate(["pedestal (ADC)", "pedestal width (ADC)"]): + if pixel_id in pixel_display: + lm = sns.lmplot( + x="T (deg)", + y=q, + hue="channel", + data=df, + scatter_kws={"alpha": 0.0}, + legend=False, + ) + + ax = lm.axes[0][0] + fig = ax.get_figure() + fig.subplots_adjust(top=0.9) + # clean legend to avoid duplication + for _artist in ax.lines + ax.collections + ax.patches + ax.images: + _artist.set_label(s=None) + + sns.violinplot( + x="T (deg)", + y=q, + hue="channel", + split=True, + native_scale=True, + data=df, + ax=ax, + ) + + # extract slope and intercept for display and further usage + colors = sns.color_palette() + for s, channel in enumerate(["hg", "lg"]): + filtered_df = df[df["channel"] == channel] + slope, intercept, r_value, pv, se = stats.linregress( + filtered_df["T (deg)"], + filtered_df[q], + ) + if k == 0 and s == 0: + slopes_ped_hg[pixel_id] = slope + elif k == 0 and s == 1: + slopes_ped_lg[pixel_id] = slope + elif k == 1 and s == 0: + slopes_rms_hg[pixel_id] = slope + elif k == 1 and s == 1: + slopes_rms_lg[pixel_id] = slope + if pixel_id in pixel_display: + ax.annotate( + f"y = {slope:.4f} T + {intercept:.4f}", + (0.05, 0.85 - s * 0.05), + color=colors[s], + xycoords="axes fraction", + ) + + if pixel_id in pixel_display: + ax.set_title(f"Pixel {pixel_id}, NSB OFF") + fig.savefig( + f"{os.environ['NECTARCHAIN_FIGURES']}/pixel{pixel_id}_{k}.png" + ) + del fig + plt.close() + +print(f"{flagged_pixels} pixels were flagged as bad by the pedestal estimation tool") +print( + f"{removed_pixels} pixels were removed because they had a pedestal RMS " + f"exceeding {rms_threshold}" +) + +camera = CameraGeometry.from_name("NectarCam-003") +camera = camera.transform_to(EngineeringCameraFrame()) + +for k in range(2): # quantity to treat + fig = plt.figure(figsize=(8, 8)) + fig.subplots_adjust(wspace=0.35, hspace=0.3) + for s in range(2): # HG and LG channels + if k == 0 and s == 0: + slopes = slopes_ped_hg + elif k == 0 and s == 1: + slopes = slopes_ped_lg + elif k == 1 and s == 0: + slopes = slopes_rms_hg + elif k == 1 and s == 1: + slopes = slopes_rms_lg + # histograms + ax = plt.subplot(2, 2, s + 1) + plt.hist(slopes, bins="auto", histtype="step") + if k == 1: + ax.set_yscale("log") + if k == 0 and s == 0: + ax.set_title("Pedestal slope (ADC/deg), hg, NSB OFF") + elif k == 0 and s == 1: + ax.set_title("Pedestal slope (ADC/deg), lg, NSB OFF") + elif k == 1 and s == 0: + ax.set_title("Pedestal width slope (ADC/deg), hg, NSB OFF") + elif k == 1 and s == 1: + ax.set_title("Pedestal width slope (ADC/deg), lg, NSB OFF") + plt.axvline(np.nanmean(slopes), linestyle=":") + plt.axvline(np.nanmean(slopes) - np.nanstd(slopes), linestyle=":") + plt.axvline(np.nanmean(slopes) + np.nanstd(slopes), linestyle=":") + ax.annotate( + r"mean: {:.3f} $\pm$ {:.3f}".format(np.nanmean(slopes), np.nanstd(slopes)), + (0.05, 0.85 - s * 0.05), + xycoords="axes fraction", + ) + # camera display + ax = plt.subplot(2, 2, s + 3) + disp = CameraDisplay(camera) + disp.image = slopes + if k == 0: + disp.set_limits_minmax(-0.3, 0.3) + elif k == 1: + disp.set_limits_minmax(-0.07, 0.07) + disp.cmap = "Spectral_r" + disp.add_colorbar() + disp.update() + ax.set_title("") + fig.savefig(f"{os.environ['NECTARCHAIN_FIGURES']}/camera_{k}.png") diff --git a/src/nectarchain/user_scripts/ltibaldo/plot_pedestal_output.py b/src/nectarchain/user_scripts/ltibaldo/plot_pedestal_output.py index b0404100..6b5bbf2e 100644 --- a/src/nectarchain/user_scripts/ltibaldo/plot_pedestal_output.py +++ b/src/nectarchain/user_scripts/ltibaldo/plot_pedestal_output.py @@ -22,7 +22,7 @@ # fill results for plotting i = 0 for result in h5file.root.__members__: - table = h5file.root[result]["NectarCAMPedestalContainer"][0] + table = h5file.root[result]["NectarCAMPedestalContainer_0"][0] wf = table["pedestal_mean_hg"][table["pixels_id"] == pixel_id][0] std = table["pedestal_std_hg"][table["pixels_id"] == pixel_id][0] if result == "data_combined":