diff --git a/lstchain/data/lstchain_lhfit_config.json b/lstchain/data/lstchain_lhfit_config.json index 3204da956b..39e742d683 100644 --- a/lstchain/data/lstchain_lhfit_config.json +++ b/lstchain/data/lstchain_lhfit_config.json @@ -291,6 +291,9 @@ ["type", "*", 0.0], ["type", "LST_LST_LSTCam", 0.0] ], + "spatial_selection": "dvr", + "dvr_pic_threshold": 8, + "dvr_pix_for_full_image": 500, "sigma_space": 3, "sigma_time": 4, "time_before_shower": [ diff --git a/lstchain/reco/r0_to_dl1.py b/lstchain/reco/r0_to_dl1.py index 0eef78ba35..baa8a93006 100644 --- a/lstchain/reco/r0_to_dl1.py +++ b/lstchain/reco/r0_to_dl1.py @@ -406,6 +406,7 @@ def r0_to_dl1( metadata=metadata, ) nsb_tuning = False + nsb_tuning_args = None if 'waveform_nsb_tuning' in config.keys(): nsb_tuning = config['waveform_nsb_tuning']['nsb_tuning'] if nsb_tuning: @@ -435,6 +436,7 @@ def r0_to_dl1( + str(np.asarray(nsb_original)[allowed_tel]) + 'GHz to {0:d}%'.format(int(nsb_tuning_ratio * 100 + 100.5)) + ' for telescopes ids ' + str(config['source_config']['LSTEventSource']['allowed_tels'])) + nsb_tuning_args = [nsb_tuning_ratio, nsb_original, pulse_template, charge_spe_cumulative_pdf] else: logger.warning('NSB tuning on waveform active in config but file is real data, option will be ignored') nsb_tuning = False @@ -443,10 +445,13 @@ def r0_to_dl1( if 'lh_fit_config' in config.keys(): lhfit_fitter_config = {'TimeWaveformFitter': config['lh_fit_config']} lhfit_fitter = TimeWaveformFitter(subarray=subarray, config=Config(lhfit_fitter_config)) - if lhfit_fitter_config['TimeWaveformFitter']['use_interleaved'] and not is_simu: + if lhfit_fitter_config['TimeWaveformFitter']['use_interleaved']: tmp_source = EventSource(input_url=input_filename, - config=Config(config["source_config"])) - lhfit_fitter.get_ped_from_interleaved(tmp_source) + config=Config(config["source_config"])) + if is_simu: + lhfit_fitter.get_ped_from_true_signal_less(tmp_source, nsb_tuning_args) + else: + lhfit_fitter.get_ped_from_interleaved(tmp_source) del tmp_source # initialize the writer of the interleaved events diff --git a/lstchain/reco/reconstructor.py b/lstchain/reco/reconstructor.py index 837a88914b..492694af50 100644 --- a/lstchain/reco/reconstructor.py +++ b/lstchain/reco/reconstructor.py @@ -4,18 +4,18 @@ import numpy as np import astropy.units as u -from lstchain.data.normalised_pulse_template import NormalizedPulseTemplate - from ctapipe.containers import EventType from ctapipe.core import TelescopeComponent -from ctapipe.core.traits import Bool, Float, FloatTelescopeParameter, Int +from ctapipe.core.traits import Bool, Float, FloatTelescopeParameter, Int, Unicode + +from lstchain.data.normalised_pulse_template import NormalizedPulseTemplate +from lstchain.image.modifier import tune_nsb_on_waveform from lstchain.io.lstcontainers import DL1LikelihoodParametersContainer from lstchain.reco.reconstructorCC import log_pdf as log_pdf logger = logging.getLogger(__name__) - class TimeWaveformFitter(TelescopeComponent): """ Class used to perform event reconstruction by fitting of a model on waveforms. @@ -24,9 +24,21 @@ class TimeWaveformFitter(TelescopeComponent): allow_none=False).tag(config=True) crosstalk = FloatTelescopeParameter(default_value=0, help='Average pixel crosstalk.', allow_none=False).tag(config=True) - sigma_space = Float(4, help='Size of the region on which the fit is performed relative to the image extension.', + spatial_selection = Unicode(default_value='dvr', + help='Method to select pixels to perform the likelihood fit. Can be:\n' + 'hillas : use the hillas length and width times sigma_space as an ellipsis.\n' + 'dvr : use data volume reduction logic with dvr_pic_threshold and' + ' dvr_pix_for_full_image', + allow_none=False).tag(config=True) + dvr_pic_threshold = Int(default_value=8, help='Pixel charge threshold for dvr like pixel selection.', + allow_none=False).tag(config=True) + dvr_pix_for_full_image = Int(default_value=500, help='Number of selected pixels above which all are kept.', + allow_none=False).tag(config=True) + sigma_space = Float(default_value=4, + help='Size of the region on which the fit is performed relative to the image extension.', allow_none=False).tag(config=True) - sigma_time = Float(3, help='Time window on which the fit is performed relative to the image temporal extension.', + sigma_time = Float(default_value=3, + help='Time window on which the fit is performed relative to the image temporal extension.', allow_none=False).tag(config=True) time_before_shower = FloatTelescopeParameter(default_value=10, help='Additional time at the start of the fit temporal window.', @@ -34,12 +46,13 @@ class TimeWaveformFitter(TelescopeComponent): time_after_shower = FloatTelescopeParameter(default_value=20, help='Additional time at the end of the fit temporal window.', allow_none=False).tag(config=True) - no_asymmetry = Bool(False, help='If true, the asymmetry of the spatial model is fixed to 0.', + no_asymmetry = Bool(default_value=False, help='If true, the asymmetry of the spatial model is fixed to 0.', allow_none=False).tag(config=True) - use_interleaved = Bool(None, help='If true, the std deviation of pedestals and dimmed pixels are estimated on ' - 'interleaved events', allow_none=True).tag(config=True) - n_peaks = Int(0, help='Maximum brightness (p.e.) for which the full likelihood computation is used. ' - 'If the Poisson term for Np.e.>n_peak is more than 1e-6 a Gaussian approximation is used.', + use_interleaved = Bool(default_value=None, help='If true, the std deviation of pedestals and dimmed pixels are ' + 'estimated on interleaved events', allow_none=True).tag(config=True) + n_peaks = Int(default_value=0, + help='Maximum brightness (p.e.) for which the full likelihood computation is used. ' + 'If the Poisson term for Np.e.>n_peak is more than 1e-6 a Gaussian approximation is used.', allow_none=False).tag(config=True) bound_charge_factor = FloatTelescopeParameter(default_value=4, help='Maximum relative change to the fitted charge parameter.', @@ -66,11 +79,12 @@ class TimeWaveformFitter(TelescopeComponent): default_seed_v_cm = FloatTelescopeParameter(default_value=40, help='Default starting value of v_cm when the seed extraction failed.', allow_none=False).tag(config=True) - verbose = Int(0, help='4 - used for tests: create debug plots\n' - '3 - create debug plots, wait for input after each event, increase minuit verbose level\n' - '2 - create debug plots, increase minuit verbose level\n' - '1 - increase minuit verbose level\n' - '0 - silent', allow_none=False).tag(config=True) + verbose = Int(default_value=0, + help='4 - used for tests: create debug plots\n' + '3 - create debug plots, wait for input after each event, increase minuit verbose level\n' + '2 - create debug plots, increase minuit verbose level\n' + '1 - increase minuit verbose level\n' + '0 - silent', allow_none=False).tag(config=True) def __init__(self, subarray, config=None, parent=None, **kwargs): super().__init__(subarray=subarray, config=config, parent=parent, **kwargs) @@ -81,7 +95,7 @@ def __init__(self, subarray, config=None, parent=None, **kwargs): self.template_dict[tel_id] = NormalizedPulseTemplate.load_from_eventsource( subarray.tel[tel_id].camera.readout) self.template_time_of_max_dict[tel_id] = self.template_dict[tel_id].compute_time_of_max() - poisson_peaks = np.arange(self.n_peaks+1, dtype=int) + poisson_peaks = np.arange(self.n_peaks + 1, dtype=int) poisson_peaks[0] = 1 self.factorial = np.cumprod(poisson_peaks, dtype='u8') # Find the transition charge between full likelihood computation and Gaussian approximation @@ -90,8 +104,8 @@ def __init__(self, subarray, config=None, parent=None, **kwargs): transition_charges = {} for config_crosstalk in self.crosstalk: # if n_peaks is set to 0, only the Gaussian approximation is used - transition_charges[config_crosstalk[2]] = 0.0 if self.n_peaks == 0\ - else self.find_transition_charge(config_crosstalk[2], 1e-2/self.n_peaks) + transition_charges[config_crosstalk[2]] = 0.0 if self.n_peaks == 0 \ + else self.find_transition_charge(config_crosstalk[2], 1e-2 / self.n_peaks) self.transition_charges = {} for tel_id in subarray.tel: self.transition_charges[tel_id] = transition_charges[self.crosstalk.tel[tel_id]] @@ -107,9 +121,7 @@ def __init__(self, subarray, config=None, parent=None, **kwargs): def get_ped_from_interleaved(self, source): """ - Parameters - ---------- - source + Use interleaved events to extract the pixel-wise pedestal variance from data. """ self.error = {} waveforms = {} @@ -121,13 +133,46 @@ def get_ped_from_interleaved(self, source): for tel_id in event.r1.tel.keys(): waveforms[tel_id].append(event.r1.tel[tel_id].waveform.squeeze()) for tel_id, tel_waveforms in waveforms.items(): - x=np.concatenate(np.asarray(tel_waveforms), axis=1) - std=[] + x = np.concatenate(np.asarray(tel_waveforms), axis=1) + std = [] for elt in x: std.append(np.nanstd(elt)) self.error[tel_id] = std self.allowed_pixels = (std > 0.5 * np.median(std)) + def get_ped_from_true_signal_less(self, source, nsb_tuning_args): + """ + Use pixels with no true signal in MC to extract the pedestal variance, assumed uniform in the camera. + """ + self.error = {} + waveforms = {} + n_pix = {} + dt = {} + for tel_id in self.subarray.tel: + waveforms[tel_id] = [] + readout = self.subarray.tel[tel_id].camera.readout + sampling_rate = readout.sampling_rate.to_value(u.GHz) + dt[tel_id] = (1.0 / sampling_rate) + waveforms[tel_id] = [] + for i, event in enumerate(source): + if i > 50 and len(waveforms) > 1000: + break + if event.trigger.event_type == EventType.SUBARRAY: + for tel_id in event.r1.tel.keys(): + if i == 0: + n_pix[tel_id] = event.r1.tel[tel_id].waveform.shape[0] + mask = event.simulation.tel[tel_id].true_image == 0 + wave = event.r1.tel[tel_id].waveform + if nsb_tuning_args is not None: + selected_gains = event.r1.tel[tel_id].selected_gain_channel + mask_high = (selected_gains == 0) + tune_nsb_on_waveform(wave, nsb_tuning_args[0], nsb_tuning_args[1][tel_id] * u.GHz, + dt[tel_id] * u.ns, nsb_tuning_args[2], mask_high, nsb_tuning_args[3]) + waveforms[tel_id].append(wave[mask].flatten()) + for tel_id, tel_waveforms in waveforms.items(): + self.error[tel_id] = np.full(n_pix[tel_id], np.nanstd(np.concatenate(tel_waveforms))) + self.allowed_pixels = np.ones(n_pix[tel_id], dtype=bool) + def call_setup(self, event, telescope_id, dl1_container): """ Extract all event dependent quantities used for the fit. @@ -155,7 +200,7 @@ def call_setup(self, event, telescope_id, dl1_container): pix_x = geometry.pix_x.to_value(unit) pix_y = geometry.pix_y.to_value(unit) r_max = geometry.guess_radius().to_value(unit) - pix_radius = np.sqrt(geometry.pix_area[0].to_value(unit ** 2)/np.pi) # find linear size of a pixel + pix_radius = np.sqrt(geometry.pix_area[0].to_value(unit ** 2) / np.pi) # find linear size of a pixel readout = self.subarray.tel[telescope_id].camera.readout sampling_rate = readout.sampling_rate.to_value(u.GHz) dt = (1.0 / sampling_rate) @@ -196,7 +241,7 @@ def call_setup(self, event, telescope_id, dl1_container): # With current likelihood computation, order and type of the parameters are important start_parameters = {'charge': dl1_container.intensity, 't_cm': dl1_container.intercept - - self.template_time_of_max_dict[telescope_id], + - self.template_time_of_max_dict[telescope_id], 'x_cm': start_x_cm.to_value(unit), 'y_cm': start_y_cm.to_value(unit), 'length': start_length, @@ -235,7 +280,8 @@ def call_setup(self, event, telescope_id, dl1_container): 'rl': (rl_min, rl_max) } - mask_pixel, mask_time = self.clean_data(pix_x, pix_y, pix_radius, times, start_parameters, telescope_id) + mask_pixel, mask_time = self.clean_data(image, geometry, pix_x, pix_y, pix_radius, times, start_parameters, + telescope_id) mask_pixel = mask_pixel & self.allowed_pixels spatial_ones = np.ones(np.sum(mask_pixel)) @@ -287,8 +333,8 @@ def __call__(self, event, telescope_id, dl1_container): focal_length = self.subarray.tel[telescope_id].optics.equivalent_focal_length angle_dist_eq = [(u.rad, u.m, lambda x: np.tan(x) * focal_length.to_value(u.m), lambda x: np.arctan(x / focal_length.to_value(u.m))), - (u.rad**2, u.m**2, lambda x: (np.tan(np.sqrt(x)) * focal_length.to_value(u.m))**2, - lambda x: (np.arctan(np.sqrt(x) / focal_length.to_value(u.m)))**2)] + (u.rad ** 2, u.m ** 2, lambda x: (np.tan(np.sqrt(x)) * focal_length.to_value(u.m)) ** 2, + lambda x: (np.arctan(np.sqrt(x) / focal_length.to_value(u.m))) ** 2)] with u.set_enabled_equivalencies(angle_dist_eq): self.start_parameters = None self.names_parameters = None @@ -299,11 +345,13 @@ def __call__(self, event, telescope_id, dl1_container): return self.predict(unit_cam, fit_params) - def clean_data(self, pix_x, pix_y, pix_radius, times, start_parameters, telescope_id): + def clean_data(self, image, geometry, pix_x, pix_y, pix_radius, times, start_parameters, telescope_id): """ Method used to select pixels and time samples used in the fitting procedure. - The spatial selection takes pixels in an ellipsis obtained from the seed Hillas parameters extended by one pixel - size and multiplied by a factor sigma_space. + The spatial selection can be done in one of two methods: + - takes pixels in an ellipsis obtained from the seed Hillas parameters extended by one pixel size and multiplied + by a factor sigma_space. + - use all pixels kept by the data volume reduction algorythm The temporal selection takes a time window centered on the seed time of center of mass and of duration equal to the time of propagation of the signal along the length of the ellipsis times a factor sigma_time. An additional fixed duration is also added before and after this time window through the time_before_shower and @@ -311,6 +359,10 @@ def clean_data(self, pix_x, pix_y, pix_radius, times, start_parameters, telescop Parameters ---------- + image : array_like + Charge in each pixel + geometry: `ctapipe.instrument.CameraGeometry` + Camera geometry pix_x, pix_y: array-like Pixels positions pix_radius: float @@ -326,20 +378,30 @@ def clean_data(self, pix_x, pix_y, pix_radius, times, start_parameters, telescop Mask used to select pixels and times for the fit """ - x_cm = start_parameters['x_cm'] - y_cm = start_parameters['y_cm'] length = start_parameters['length'] - width = start_parameters['wl'] * length - psi = start_parameters['psi'] - - dx = pix_x - x_cm - dy = pix_y - y_cm - - lon = dx * np.cos(psi) + dy * np.sin(psi) - lat = dx * np.sin(psi) - dy * np.cos(psi) - - mask_pixel = ((lon / (length + pix_radius)) ** 2 + ( - lat / (width + pix_radius)) ** 2) < self.sigma_space ** 2 + if self.spatial_selection == 'hillas': + x_cm = start_parameters['x_cm'] + y_cm = start_parameters['y_cm'] + width = start_parameters['wl'] * length + psi = start_parameters['psi'] + + dx = pix_x - x_cm + dy = pix_y - y_cm + + lon = dx * np.cos(psi) + dy * np.sin(psi) + lat = dx * np.sin(psi) - dy * np.cos(psi) + + mask_pixel = ((lon / (length + pix_radius)) ** 2 + ( + lat / (width + pix_radius)) ** 2) < self.sigma_space ** 2 + elif self.spatial_selection == 'dvr': + mask_pixel = (image > self.dvr_pic_threshold) + # we add-up (sum) the selected-pixel-wise map of neighbors, to find + # those who appear at least once (>0). Those should be added: + additional_pixels = (np.sum(geometry.neighbor_matrix[mask_pixel], axis=0) > 0) + mask_pixel |= additional_pixels + # if more than min_npixels_for_full_event were selected, keep whole camera: + if mask_pixel.sum() > self.dvr_pix_for_full_image: + mask_pixel = np.array(geometry.n_pixels * [True]) v = start_parameters['v'] t_start = (start_parameters['t_cm'] @@ -395,6 +457,7 @@ def fit(self, fit_params): Parameters used to compute the likelihood but not fitted """ + def f(*args): return -2 * self.log_likelihood(*args, fit_params=fit_params) @@ -433,10 +496,10 @@ def predict(self, unit_cam, fit_params): self.fit(fit_params) container.lhfit_TS = self.fcn - container.lhfit_x = (self.end_parameters['x_cm']*unit_cam).to(u.m) - container.lhfit_x_uncertainty = (self.error_parameters['x_cm']*unit_cam).to(u.m) - container.lhfit_y = (self.end_parameters['y_cm']*unit_cam).to(u.m) - container.lhfit_y_uncertainty = (self.error_parameters['y_cm']*unit_cam).to(u.m) + container.lhfit_x = (self.end_parameters['x_cm'] * unit_cam).to(u.m) + container.lhfit_x_uncertainty = (self.error_parameters['x_cm'] * unit_cam).to(u.m) + container.lhfit_y = (self.end_parameters['y_cm'] * unit_cam).to(u.m) + container.lhfit_y_uncertainty = (self.error_parameters['y_cm'] * unit_cam).to(u.m) container.lhfit_r = np.sqrt(container.lhfit_x ** 2 + container.lhfit_y ** 2) container.lhfit_phi = np.arctan2(container.lhfit_y, container.lhfit_x) if self.end_parameters['psi'] > np.pi: diff --git a/lstchain/tests/test_lstchain.py b/lstchain/tests/test_lstchain.py index c742853a0b..d88bab4b76 100644 --- a/lstchain/tests/test_lstchain.py +++ b/lstchain/tests/test_lstchain.py @@ -95,6 +95,9 @@ def test_r0_to_dl1_lhfit_mc(tmp_path, mc_gamma_testfile): ["type", "*", 0.0], ["type", "LST_LST_LSTCam", 0.0] ], + "spatial_selection": "hillas", + "dvr_pic_threshold": 8, + "dvr_pix_for_full_image": 500, "sigma_space": 3, "sigma_time": 4, "time_before_shower": [ @@ -122,6 +125,11 @@ def test_r0_to_dl1_lhfit_mc(tmp_path, mc_gamma_testfile): config['lh_fit_config']["no_asymmetry"] = True config['lh_fit_config']["verbose"] = 0 r0_to_dl1(mc_gamma_testfile, custom_config=config, output_filename=tmp_path / "tmp.h5") + os.remove(tmp_path / "tmp.h5") + config['lh_fit_config']["spatial_selection"] = 'dvr' + config['lh_fit_config']["use_interleaved"] = True + config['waveform_nsb_tuning']['nsb_tuning'] = True + r0_to_dl1(mc_gamma_testfile, custom_config=config, output_filename=tmp_path / "tmp.h5") @pytest.mark.private_data @@ -139,6 +147,9 @@ def test_r0_to_dl1_lhfit_observed(tmp_path): ["type", "*", 0.0], ["type", "LST_LST_LSTCam", 0.0] ], + "spatial_selection": "hillas", + "dvr_pic_threshold": 8, + "dvr_pix_for_full_image": 500, "sigma_space": 3, "sigma_time": 4, "time_before_shower": [ @@ -156,6 +167,9 @@ def test_r0_to_dl1_lhfit_observed(tmp_path): "verbose": 0 } r0_to_dl1(test_r0_path, custom_config=config, output_filename=tmp_path / "tmp2.h5") + os.remove(tmp_path / "tmp2.h5") + config['lh_fit_config']["spatial_selection"] = 'dvr' + r0_to_dl1(test_r0_path, custom_config=config, output_filename=tmp_path / "tmp2.h5") def test_content_dl1(simulated_dl1_file):