diff --git a/lstchain/calib/camera/sampling_interval_coefficient_calculate.py b/lstchain/calib/camera/sampling_interval_coefficient_calculate.py new file mode 100644 index 0000000000..0ac0e35932 --- /dev/null +++ b/lstchain/calib/camera/sampling_interval_coefficient_calculate.py @@ -0,0 +1,283 @@ +import numpy as np +from astropy.io import fits +from numba import jit, prange +import logging + +from ctapipe.core import Component +from ctapipe_io_lst.constants import ( + N_PIXELS, N_SAMPLES, N_CAPACITORS_CHANNEL +) +import matplotlib.pyplot as plt +from ctapipe.instrument import CameraGeometry +from ctapipe.visualization import CameraDisplay + +logger = logging.getLogger(__name__) + +THRESHOLD_PULSE_ADC = 100 +THRESHOLD_PULS_PIX_NUM = 1800 +CHANNEL = ["HG", "LG"] +CRITERIA = [0.07, 0.05] + + +__all__ = ['SamplingIntervalCoefficientPeakCount', 'SamplingIntervalCoefficientCalculate'] + + +class SamplingIntervalCoefficientPeakCount(Component): + """ + The SamplingIntervalCoefficientPeakCount class to create a peak count table for DRS4 sampling interval calibration. + """ + def __init__(self, gain): + + self.peak_count = np.zeros([N_PIXELS, N_CAPACITORS_CHANNEL], dtype=np.uint16) + self.fc_count = np.zeros([N_PIXELS, N_CAPACITORS_CHANNEL], dtype=np.uint16) + self.gain = gain + + def increment_peak_count(self, event, tel_id, r0_r1_calibrator): + + waveform = event.r1.tel[tel_id].waveform + + first_capacitors = r0_r1_calibrator.first_cap + r1_sample_start = r0_r1_calibrator.r1_sample_start.tel[tel_id] + + # There are some events without test pulses for almost all modules. + # In those events, a few modules have pulse signals of which pulse shapes are not expected... + # For the moment, such events are just removed in the anlysis + pulse_event_flag = np.sum(np.max(waveform[self.gain], axis=1) > THRESHOLD_PULSE_ADC) > THRESHOLD_PULS_PIX_NUM + + if pulse_event_flag : + # find pulse peak position + pulse_flag = np.max(waveform[self.gain], axis=1) > THRESHOLD_PULSE_ADC + pulse_peak = np.argmax(waveform[self.gain], axis=1) + fc = first_capacitors[tel_id][self.gain] + pulse_peak_abs_pos = (fc + r1_sample_start + pulse_peak) % N_CAPACITORS_CHANNEL + + # increment peak count array in case that test pulses exist in the readout window + self.peak_count[np.arange(N_PIXELS), pulse_peak_abs_pos] += pulse_flag + self.fc_count[np.arange(N_PIXELS), fc % N_CAPACITORS_CHANNEL] += pulse_flag + + +class SamplingIntervalCoefficientCalculate(Component): + """ + The SamplingIntervalCoefficientCalculate class to create a sampling interval coefficient table for LST readout system using chip DRS4. + """ + def __init__(self, gain): + + self.peak_count_stack ={} + self.fc_count_stack = {} + self.sampling_interval_coefficient = {} + self.gain = gain + + def stack_peak_count_fits(self, file_list): + # stack peak count fits files, and sort stacked peak counts by gain, run_id + + for single_file in file_list: + + run_id, gain_in_file, peak_count, fc_count = load_peak_count_fits_file(single_file) + + # select only a given gain + if gain_in_file == self.gain: + if not run_id in self.peak_count_stack.keys(): + self.peak_count_stack[run_id] = np.zeros([N_PIXELS, N_CAPACITORS_CHANNEL], dtype=np.uint16) + self.fc_count_stack[run_id] = np.zeros([N_PIXELS, N_CAPACITORS_CHANNEL], dtype=np.uint16) + + self.peak_count_stack[run_id] += peak_count + self.fc_count_stack[run_id] += fc_count + + def convert_to_samp_interval_coefficient(self): + # convert peak counts to sampling interval coefficient + + for run_id in self.peak_count_stack.keys(): + self.sampling_interval_coefficient[run_id] = np.zeros([N_PIXELS, N_CAPACITORS_CHANNEL + N_SAMPLES]) + + stack_events_num = np.sum(self.peak_count_stack[run_id], axis=1) + stack_events_num = stack_events_num.reshape(1, N_PIXELS).T + + self.sampling_interval_coefficient[run_id][:,:N_CAPACITORS_CHANNEL] = self.peak_count_stack[run_id] / stack_events_num * N_CAPACITORS_CHANNEL + self.sampling_interval_coefficient[run_id][:,N_CAPACITORS_CHANNEL:] = self.sampling_interval_coefficient[run_id][:, :N_SAMPLES] + + def set_charge_array(self): + # set charge array to compute charge resolutions using calculated coefficients with different observation runs + + # before calibration + self.charge_array_before_corr = [] + + # after calibration + self.charge_array_after_corr ={} + + for run_id in self.peak_count_stack.keys(): + self.charge_array_after_corr[run_id] = [] + + def calc_charge(self, event, tel_id, r0_r1_calibrator, extractor): + + waveform = event.r1.tel[tel_id].waveform + + # There are some events without test pulses for almost all modules. + # In those events, a few modules have pulse signals of which pulse shapes are not expected... + # For the moment, such events are just removed in the anlysis + pulse_event_flag = np.sum(np.max(waveform[self.gain], axis=1) > THRESHOLD_PULSE_ADC) > THRESHOLD_PULS_PIX_NUM + + if pulse_event_flag: + + first_capacitors = r0_r1_calibrator.first_cap + r1_sample_start = r0_r1_calibrator.r1_sample_start.tel[tel_id] + r1_sample_end = r0_r1_calibrator.r1_sample_end.tel[tel_id] + + # find pulse peak position + pulse_peak = np.argmax(waveform[self.gain], axis=1) + fc = first_capacitors[tel_id][self.gain] + + # calculate charge (before calibration) + charge, peak_time = extractor(waveform[self.gain], tel_id, np.full(N_PIXELS, self.gain, dtype=int)) + self.charge_array_before_corr.append(charge.tolist()) + + # calculate charge (after calibration with sampling interval coefficients obtained by each observation run) + for run_id in self.peak_count_stack.keys(): + sampling_interval_coefficient_event = np.zeros([N_PIXELS, r1_sample_end - r1_sample_start]) + + for pixel in range(N_PIXELS): + sampling_interval_coefficient_event[pixel] = ( + self.sampling_interval_coefficient[run_id][pixel][fc[pixel]%N_CAPACITORS_CHANNEL + r1_sample_start: + fc[pixel]%N_CAPACITORS_CHANNEL + r1_sample_end] + ) + + charge, peak_time = extractor(waveform[self.gain] * sampling_interval_coefficient_event, tel_id, np.full(N_PIXELS, self.gain, dtype=int)) + + self.charge_array_after_corr[run_id].append(charge.tolist()) + + def calc_charge_reso(self): + # calculate charge resolution to evaluate the sampling interval correction + + # before correction + self.charge_reso_array_before_corr = ( + np.std(np.array(self.charge_array_before_corr).T, axis=1)/np.mean(np.array(self.charge_array_before_corr).T, axis=1) + ) + + # after correction + self.charge_reso_array_after_corr ={} + + for run_id in self.peak_count_stack.keys(): + self.charge_reso_array_after_corr[run_id] = ( + np.std(np.array(self.charge_array_after_corr[run_id]).T, axis=1)/np.mean(np.array(self.charge_array_after_corr[run_id]).T, axis=1) + ) + + def verify(self): + # verify the obtained sampling interval coefficients + # Coefficients which give the best charge resolution are adopted as the final ones + # If the charge resolution is worse than required criteria (<7% for high gain, <5% for low gain), coefficients are filled with ones. + + self.charge_reso_final = np.zeros(N_PIXELS) + self.used_run = np.zeros(N_PIXELS, dtype=np.uint16) + self.sampling_interval_coefficient_final = np.zeros([N_PIXELS, N_CAPACITORS_CHANNEL]) + + n_runs = len(self.charge_reso_array_after_corr.keys()) + run_id_array = np.zeros(n_runs, dtype=np.uint16) + charge_reso_array_after_corr_all = np.zeros([n_runs, N_PIXELS]) + + for i, ikey in enumerate(self.charge_reso_array_after_corr.keys()): + run_id_array[i] = ikey + charge_reso_array_after_corr_all[i] = self.charge_reso_array_after_corr[ikey] + + for pixel in range(N_PIXELS): + min_charge_reso_arg = np.argmin(charge_reso_array_after_corr_all.T[pixel]) + + self.charge_reso_final[pixel] = charge_reso_array_after_corr_all[min_charge_reso_arg, pixel] + self.used_run[pixel] = run_id_array[min_charge_reso_arg] + + criteria_flag = self.charge_reso_final[pixel] < CRITERIA[self.gain] + + if criteria_flag: + self.sampling_interval_coefficient_final[pixel] = self.sampling_interval_coefficient[self.used_run[pixel]][pixel, :N_CAPACITORS_CHANNEL] + else: + logger.info(f'[charge resolution] {CHANNEL[self.gain]} pixel {pixel}: {self.charge_reso_final[pixel]} > {CRITERIA[self.gain]}') + self.sampling_interval_coefficient_final[pixel] = np.ones(N_CAPACITORS_CHANNEL) + + def plot_results(self, verify_data_path, output_file): + + plt.rcParams['font.size']=15 + camera = CameraGeometry.from_name('LSTCam-002') + + num_run = len(self.charge_reso_array_after_corr.keys()) + + plt.figure(figsize=(20, 6*(num_run+1))) + plt.tight_layout() + plt.subplots_adjust(top=0.94) + plt.suptitle("{} (data for the verification: {})".format(CHANNEL[self.gain], verify_data_path)) + + for i, run_id in enumerate(self.charge_reso_array_after_corr.keys()): + + # charge resolution scatter plot (before correction vs after correction) + plt.subplot(num_run + 1, 3, i*3+1) + plt.scatter(self.charge_reso_array_before_corr, self.charge_reso_array_after_corr[run_id]) + plt.plot([0, 0.15], [0, 0.15], color='black', ls='--') + plt.hlines(CRITERIA[self.gain], 0, 0.15, color='red', ls='--') + plt.xlim(0,.15) + plt.ylim(0,.15) + plt.title('charge resolution (run {})'.format(run_id)) + plt.xlabel('before calibration') + plt.ylabel('after calibration') + plt.grid(color='gray', alpha=0.5) + + # charge resolution histogram + plt.subplot(num_run + 1, 3, i*3+2) + plt.hist(self.charge_reso_array_after_corr[run_id], range=[0, 0.1], bins=50) + plt.vlines(CRITERIA[self.gain], 0.05, 500, color='red', ls='--') + plt.title('charge resolution (after calib.) (run {})'.format(run_id)) + plt.xlabel('charge resolution') + plt.yscale('log') + plt.grid(color='gray', alpha=0.5) + + # camera map + plt.subplot(num_run + 1, 3, i*3+3) + disp = CameraDisplay(camera) + disp.image = self.charge_reso_array_after_corr[run_id] + + plt.grid(color='gray', alpha=0.5) + disp.set_limits_minmax(0, CRITERIA[self.gain]) + disp.highlight_pixels(self.charge_reso_array_after_corr[run_id]>CRITERIA[self.gain], color='red') + disp.add_colorbar() + + # Final result + plt.subplot(num_run + 1, 3, num_run*3+1) + plt.scatter(self.charge_reso_array_before_corr, self.charge_reso_final) + plt.plot([0, 0.15], [0, 0.15], color='black', ls='--') + plt.hlines(CRITERIA[self.gain], 0, 0.15, color='red', ls='--') + plt.xlim(0,.15) + plt.ylim(0,.15) + plt.title('charge resolution (Final)'.format(run_id)) + plt.xlabel('before calibration') + plt.ylabel('after calibration') + plt.grid(color='gray', alpha=0.5) + + plt.subplot(num_run + 1, 3, num_run*3+2) + plt.hist(self.charge_reso_final, range=[0, 0.1], bins=50) + plt.vlines(CRITERIA[self.gain], 0.05, 500, color='red', ls='--') + plt.title('charge resolution (after calib.) (Final)'.format(run_id)) + plt.xlabel('charge resolution') + plt.yscale('log') + plt.grid(color='gray', alpha=0.5) + + plt.subplot(num_run + 1, 3, num_run*3+3) + disp = CameraDisplay(camera) + disp.image = self.charge_reso_final + + plt.grid(color='gray', alpha=0.5) + disp.set_limits_minmax(0, CRITERIA[self.gain]) + disp.highlight_pixels(self.charge_reso_final>CRITERIA[self.gain], color='red') + disp.add_colorbar() + + for i, bad_pixel in enumerate(np.where(self.charge_reso_final>CRITERIA[self.gain])[0]): + plt.text(-1.2, 1.5-0.2*i, 'pix{}: {:.3f}'.format(bad_pixel, self.charge_reso_final[bad_pixel]), color='red') + + + plt.savefig(output_file) + +def load_peak_count_fits_file(input_file): + + hdulist = fits.open(input_file) + + run_id = hdulist[0].header['run_id'] + gain = hdulist[0].header['gain'] + peak_count = hdulist[1].data + fc_count = hdulist[2].data + + return run_id, gain, peak_count, fc_count diff --git a/lstchain/tools/lstchain_create_sampling_interval_coefficient_file.py b/lstchain/tools/lstchain_create_sampling_interval_coefficient_file.py new file mode 100644 index 0000000000..be736cc838 --- /dev/null +++ b/lstchain/tools/lstchain_create_sampling_interval_coefficient_file.py @@ -0,0 +1,133 @@ +import numpy as np +from astropy.io import fits +import h5py +from tqdm import tqdm +import os + +from ctapipe.core import Provenance, traits +from ctapipe.core import Tool +from ctapipe.image.extractor import ImageExtractor + +from ctapipe_io_lst import LSTEventSource +from ctapipe_io_lst.calibration import LSTR0Corrections +from ctapipe_io_lst.constants import N_GAINS, N_PIXELS, N_CAPACITORS_CHANNEL + +from lstchain.calib.camera.sampling_interval_coefficient_calculate import SamplingIntervalCoefficientCalculate +from lstchain.paths import run_info_from_filename + +class SamplingIntervalCoefficientHDFWriter(Tool): + """ + Tool that generates a sampling interval coefficient hdf file. + For getting help run: + lstchain_create_sampling_interval_coefficient_file --help + """ + + name = "SamplingIntervalCoefficientHDFWriter" + description = "Generate a sampling interval coefficient table" + + input_peak_count_fits = traits.Path( + help="Path to the generated peak count fits files by lstchain_create_peak_count_file_for_sampling_interval_calib", + directory_ok=True, + ).tag(config=True) + + output = traits.Path( + help="Path to the output fits file", + directory_ok=False, + ).tag(config=True) + + glob = traits.Unicode( + help="Filename pattern to glob files in the directory", + default_value="*.fits" + ).tag(config=True) + + gain = traits.Int( + help="Select gain channel (0: High Gain, 1: Low Gain)", + ).tag(config=True) + + progress_bar = traits.Bool( + help="show progress bar during processing", + default_value=True, + ).tag(config=True) + + aliases = { + ("i", "input"): "EventSource.input_url", + ("o", "output"): "SamplingIntervalCoefficientHDFWriter.output", + ("g", "gain"): "SamplingIntervalCoefficientHDFWriter.gain", + "input_peak_count_fits": "SamplingIntervalCoefficientHDFWriter.input_peak_count_fits", + "pedestal": "LSTR0Corrections.drs4_pedestal_path", + "max-events": "EventSource.max_events", + "r1_sample_start": "LSTR0Corrections.r1_sample_start", + "r1_sample_end": "LSTR0Corrections.r1_sample_end", + "progress_bar": "SamplingIntervalCoefficientHDFWriter.progress_bar", + } + + classes = [LSTEventSource, LSTR0Corrections, SamplingIntervalCoefficientCalculate] + + def setup(self): + + self.sampling_interval_coefficient_calculate = SamplingIntervalCoefficientCalculate(self.gain) + self.eventsource = LSTEventSource(parent=self) + self.image_extractor = ImageExtractor.from_name('LocalPeakWindowSum', subarray = self.eventsource.subarray) + + self.path_list = [str(self.input_peak_count_fits)] + + if self.input_peak_count_fits.is_dir(): + self.path_list = sorted(self.input_peak_count_fits.glob(self.glob)) + + def start(self): + + self.sampling_interval_coefficient_calculate.stack_peak_count_fits(self.path_list) + self.sampling_interval_coefficient_calculate.convert_to_samp_interval_coefficient() + self.sampling_interval_coefficient_calculate.set_charge_array() + + for i, event in enumerate(tqdm( + self.eventsource, + desc = self.eventsource.__class__.__name__, + total = self.eventsource.max_events, + unit = "event", + disable = not self.progress_bar, + )): + + if event.index.event_id % 500 == 0: + self.log.debug(f'event id = {event.index.event_id}') + + # skip the first 1000 events which would not be calibrated correctly + if i < 1000: + continue + + self.sampling_interval_coefficient_calculate.calc_charge(event, tel_id = self.eventsource.tel_id, + r0_r1_calibrator = self.eventsource.r0_r1_calibrator, extractor = self.image_extractor) + + self.sampling_interval_coefficient_calculate.calc_charge_reso() + self.sampling_interval_coefficient_calculate.verify() + + def finish(self): + + CHANNEL = ["HG", "LG"] + + output_abs_path = os.path.abspath(self.output) + output_log_pdf_file = os.path.join(os.path.dirname(output_abs_path), "sampling_interval_coefficient_results_{}.pdf".format(CHANNEL[self.gain])) + + self.sampling_interval_coefficient_calculate.plot_results(self.eventsource.input_url, output_log_pdf_file) + + with h5py.File(self.output, 'a') as hf: + + if not 'sampling_interval_coefficient' in hf: + hf.create_dataset('sampling_interval_coefficient', data = np.zeros([N_GAINS, N_PIXELS, N_CAPACITORS_CHANNEL])) + hf['sampling_interval_coefficient'][self.gain] = self.sampling_interval_coefficient_calculate.sampling_interval_coefficient_final + + if not 'used_run' in hf: + hf.create_dataset('used_run', data = np.zeros([N_GAINS, N_PIXELS]), dtype=np.uint16) + hf['used_run'][self.gain] = self.sampling_interval_coefficient_calculate.used_run + + if not 'charge_resolution' in hf: + hf.create_dataset('charge_resolution', data = np.zeros([N_GAINS, N_PIXELS])) + hf['charge_resolution'][self.gain] = self.sampling_interval_coefficient_calculate.charge_reso_final + +def main(): + exe = SamplingIntervalCoefficientHDFWriter() + exe.run() + + +if __name__ == "__main__": + main()