diff --git a/magicctapipe/image/__init__.py b/magicctapipe/image/__init__.py index 8e867fac3..320f66376 100644 --- a/magicctapipe/image/__init__.py +++ b/magicctapipe/image/__init__.py @@ -3,17 +3,12 @@ PixelTreatment, get_num_islands_MAGIC, clean_image_params, - apply_dynamic_cleaning, ) from .leakage import ( get_leakage, ) -from .modifier import ( - add_noise_in_pixels, - random_psf_smearer, - set_numba_seed, -) + __all__ = [ "MAGICClean", @@ -21,8 +16,4 @@ "get_num_islands_MAGIC", "clean_image_params", "get_leakage", - "apply_dynamic_cleaning", - "add_noise_in_pixels", - "random_psf_smearer", - "set_numba_seed", ] diff --git a/magicctapipe/image/cleaning.py b/magicctapipe/image/cleaning.py index 35c449d63..004ae139c 100644 --- a/magicctapipe/image/cleaning.py +++ b/magicctapipe/image/cleaning.py @@ -14,7 +14,6 @@ "PixelTreatment", "get_num_islands_MAGIC", "clean_image_params", - "apply_dynamic_cleaning", ] @@ -720,35 +719,3 @@ def clean_image_params(geom, image, clean, peakpos): # time_grad[tel_id] = 1 return hillas_p, leakage_p, timing_p - - -def apply_dynamic_cleaning(image, signal_pixels, threshold, fraction): - """ - Application of the dynamic cleaning - Parameters - ---------- - image: `np.ndarray` - Pixel charges - signal_pixels - threshold: `float` - Minimum average charge in the 3 brightest pixels to apply - the dynamic cleaning (else nothing is done) - fraction: `float` - Pixels below fraction * (average charge in the 3 brightest pixels) - will be removed from the cleaned image - Returns - ------- - mask_dynamic_cleaning: `np.ndarray` - Mask with the selected pixels after the dynamic cleaning - """ - - max_3_value_index = np.argsort(image)[-3:] - mean_3_max_signal = np.mean(image[max_3_value_index]) - - if mean_3_max_signal < threshold: - return signal_pixels - - dynamic_threshold = fraction * mean_3_max_signal - mask_dynamic_cleaning = (image >= dynamic_threshold) & signal_pixels - - return mask_dynamic_cleaning diff --git a/magicctapipe/image/modifier.py b/magicctapipe/image/modifier.py deleted file mode 100644 index a9ead51d3..000000000 --- a/magicctapipe/image/modifier.py +++ /dev/null @@ -1,132 +0,0 @@ -import logging - -import numpy as np -from numba import njit - - -__all__ = [ - "add_noise_in_pixels", - "random_psf_smearer", - "set_numba_seed", -] - -log = logging.getLogger(__name__) - - -# number of neighbors of completely surrounded pixels of hexagonal cameras: -N_PIXEL_NEIGHBORS = 6 -SMEAR_PROBABILITIES = np.full(N_PIXEL_NEIGHBORS, 1 / N_PIXEL_NEIGHBORS) - - -def add_noise_in_pixels( - rng, - image, - extra_noise_in_dim_pixels, - extra_bias_in_dim_pixels, - transition_charge, - extra_noise_in_bright_pixels, -): - """ - Addition of Poissonian noise to the pixels - Parameters - ---------- - rng : `numpy.random.default_rng` - Random number generator - image: `np.ndarray` - Charges (p.e.) in the camera - extra_noise_in_dim_pixels: `float` - Mean additional number of p.e. to be added (Poisson noise) to - pixels with charge below transition_charge. To be tuned by - comparing the starting MC and data - extra_bias_in_dim_pixels: `float` - Mean bias (w.r.t. original charge) of the new charge in pixels. - Should be 0 for non-peak-search pulse integrators. To be tuned by - comparing the starting MC and data - transition_charge: `float` - Border between "dim" and "bright" pixels. To be tuned by - comparing the starting MC and data - extra_noise_in_bright_pixels: `float` - Mean additional number of p.e. to be added (Poisson noise) to - pixels with charge above transition_charge. This is unbiased, - i.e. Poisson noise is introduced, and its average subtracted, - so that the mean charge in bright pixels remains unaltered. - This is because we assume that above transition_charge the - integration window is determined by the Cherenkov light, and - would not be modified by the additional NSB noise (presumably - small compared to the C-light). To be tuned by - comparing the starting MC and data - Returns - ------- - image: `np.ndarray` - Modified (noisier) image - """ - - bright_pixels = image > transition_charge - noise = np.where( - bright_pixels, extra_noise_in_bright_pixels, extra_noise_in_dim_pixels - ) - bias = np.where( - bright_pixels, - -extra_noise_in_bright_pixels, - extra_bias_in_dim_pixels - extra_noise_in_dim_pixels, - ) - - image = image + rng.poisson(noise) + bias - - return image - - -@njit(cache=True) -def set_numba_seed(seed): - np.random.seed(seed) - - -@njit(cache=True) -def random_psf_smearer(image, fraction, indices, indptr): - """ - Random PSF smearer - Parameters - ---------- - image: `np.ndarray` - Charges (p.e.) in the camera - indices : `camera_geometry.neighbor_matrix_sparse.indices` - Pixel indices. - indptr : camera_geometry.neighbor_matrix_sparse.indptr - fraction: `float` - Fraction of the light in a pixel that will be distributed among its - immediate surroundings, i.e. immediate neighboring pixels, according - to Poisson statistics. Some light is lost for pixels which are at - the camera edge and hence don't have all possible neighbors - Returns - ------- - new_image: `np.ndarray` - Modified (smeared) image - """ - - new_image = image.copy() - for pixel in range(len(image)): - if image[pixel] <= 0: - continue - - to_smear = np.random.poisson(image[pixel] * fraction) - - if to_smear == 0: - continue - - # remove light from current pixel - new_image[pixel] -= to_smear - - # add light to neighbor pixels - neighbors = indices[indptr[pixel] : indptr[pixel + 1]] - n_neighbors = len(neighbors) - - # all neighbors are equally likely to receive the charge - # we always distribute the charge into 6 neighbors, so that charge - # on the edges of the camera is lost - neighbor_charges = np.random.multinomial(to_smear, SMEAR_PROBABILITIES) - - for n in range(n_neighbors): - neighbor = neighbors[n] - new_image[neighbor] += neighbor_charges[n] - - return new_image