Source code for ctapipe.image.muon.intensity_fitter

```"""
Class for performing a HESS style 2D fit of muon images

To do:
- Deal with astropy untis better, currently stripped and no checks made
- unit tests
- create container class for output

"""
from functools import lru_cache
from math import erf

import numpy as np
from astropy import units as u
from astropy.coordinates import SkyCoord
from iminuit import Minuit
from numba import double, vectorize
from scipy.constants import alpha
from scipy.ndimage import correlate1d

from ctapipe.image.pixel_likelihood import neg_log_likelihood_approx

from ...containers import MuonEfficiencyContainer
from ...coordinates import CameraFrame, TelescopeFrame
from ...core import TelescopeComponent
from ...core.traits import FloatTelescopeParameter, IntTelescopeParameter

__all__ = ["MuonIntensityFitter"]

# ratio of the areas of the unit circle and a square of side lengths 2
CIRCLE_SQUARE_AREA_RATIO = np.pi / 4

# Sqrt of 2, as it is needed multiple times
SQRT2 = np.sqrt(2)

@vectorize([double(double, double, double)])
"""
Function for integrating the length of a chord across a circle

Parameters
----------
rho: float or ndarray
fractional distance of impact point from circle center
phi: float or ndarray in radians
rotation angles to calculate length

Returns
-------
float or ndarray:
chord length
"""
chord = 1 - (rho**2 * np.sin(phi) ** 2)
valid = chord >= 0

if not valid:
return 0

if rho <= 1.0:
# muon has hit the mirror
chord = radius * (np.sqrt(chord) + rho * np.cos(phi))
else:
# muon did not hit the mirror
chord = 2 * radius * np.sqrt(chord)

return chord

"""Perform line integration along a given axis in the mirror frame
given an impact point on the mirror

Parameters
----------
angle: ndarray
Angle along which to integrate mirror

Returns
--------
float: length from impact point to mirror edge

"""

return mirror_length

return mirror_length - hole_length

"""Calculate number of pixels of diameter ``pixel_diameter`` on the circumference
"""
circumference = 2 * np.pi * radius
n_pixels = u.Quantity(circumference / pixel_diameter)
return int(n_pixels.to_value(u.dimensionless_unscaled))

@lru_cache(maxsize=1000)
def linspace_two_pi(n_points):
return np.linspace(-np.pi, np.pi, n_points)

def create_profile(
impact_parameter,
phi,
pixel_diameter,
oversampling=3,
):
"""
Perform intersection over all angles and return length

Parameters
----------
impact_parameter: float
Impact distance from mirror center
ang: ndarray
Angles over which to integrate
phi: float
Rotation angle of muon image

Returns
-------
ndarray
Chord length for each angle
"""
circumference = 2 * np.pi * radius
pixels_on_circle = int(circumference / pixel_diameter)

ang = phi + linspace_two_pi(pixels_on_circle * oversampling)

length = correlate1d(length, np.ones(oversampling), mode="wrap", axis=0)
length /= oversampling

return ang, length

def image_prediction(
impact_parameter,
phi,
center_x,
center_y,
ring_width,
pixel_x,
pixel_y,
pixel_diameter,
oversampling=3,
min_lambda=300 * u.nm,
max_lambda=600 * u.nm,
):
"""
Parameters
----------
impact_parameter: quantity[length]
Impact distance of muon
center_x: quantity[angle]
Muon ring center in telescope frame
center_y: quantity[angle]
Muon ring center in telescope frame
Radius of muon ring in telescope frame
ring_width: quantity[angle]
Gaussian width of muon ring
pixel_x: quantity[angle]
Pixel x coordinate in telescope
pixel_y: quantity[angle]
Pixel y coordinate in telescope

Returns
-------
ndarray:
Predicted signal
"""
return image_prediction_no_units(
impact_parameter.to_value(u.m),
oversampling=oversampling,
min_lambda_m=min_lambda.to_value(u.m),
max_lambda_m=max_lambda.to_value(u.m),
)

@vectorize([double(double, double, double)])
def gaussian_cdf(x, mu, sig):
"""
Function to compute values of a given gaussians
cumulative distribution function (cdf)

Parameters
----------
x: float or ndarray
point, at which the cdf should be evaluated
mu: float or ndarray
gaussian mean
sig: float or ndarray
gaussian standard deviation

Returns
-------
float or ndarray: cdf-value at x
"""
return 0.5 * (1 + erf((x - mu) / (SQRT2 * sig)))

def image_prediction_no_units(
impact_parameter_m,
oversampling=3,
min_lambda_m=300e-9,
max_lambda_m=600e-9,
):
"""Function for producing the expected image for a given set of trial
muon parameters without using astropy units but expecting the input to
be in the correct ones.

See [chalmecalvet2013]_
"""

# First produce angular position of each pixel w.r.t muon center
ang = np.arctan2(dy, dx)

# Produce smoothed muon profile
ang_prof, profile = create_profile(
impact_parameter_m,
oversampling=oversampling,
)

# Produce gaussian weight for each pixel given ring width
# The weight is the integral of the ring's radial gaussian profile inside the
# ring's width
cdfs = gaussian_cdf(
)
gauss = cdfs[0] - cdfs[1]

# interpolate profile to find prediction for each pixel
pred = np.interp(ang, ang_prof, profile)

# Multiply by integrated emissivity between 300 and 600 nm, and rest of factors to
# get total number of photons per pixel
# ^ would be per radian, but no need to put it here, would anyway cancel out below

pred *= alpha * (min_lambda_m**-1 - max_lambda_m**-1)
# multiply by angle (in radians) subtended by pixel width as seen from ring center

# multiply by gaussian weight, to account for "fraction of muon ring" which falls
# within the pixel
pred *= gauss

# Now it would be the total light in an area S delimited by: two radii of the
# ring, tangent to the sides of the pixel in question, and two circles concentric
# with the ring, also tangent to the sides of the pixel.
# A rough correction, assuming pixel is round, is introduced here:
# [pi*(pixel_diameter/2)**2]/ S. Actually, for the large rings (relative to pixel
# size) we are concerned with, a good enough approximation is the ratio between a
# circle's area and that of the square whose side is equal to the circle's
# diameter. In any case, since in the end we do a data-MC comparison of the muon
# ring analysis outputs, it is not critical that this value is exact.

pred *= CIRCLE_SQUARE_AREA_RATIO

return pred

def build_negative_log_likelihood(
image,
telescope_description,
oversampling,
min_lambda,
max_lambda,
spe_width,
pedestal,
):
"""Create an efficient negative log_likelihood function that does
not rely on astropy units internally by defining needed values as closures
in this function.

The likelihood is the gaussian approximation,
i.e. the unnumbered equation on page 22 between (24) and (25), from [denaurois2009]_

The logarithm of the likelihood is calculated analytically as far as possible
and terms constant under differentation are discarded.
"""

# get all the neeed values and transform them into appropriate units
optics = telescope_description.optics
mirror_area = optics.mirror_area.to_value(u.m**2)

focal_length = optics.equivalent_focal_length

cam = telescope_description.camera.geometry
camera_frame = CameraFrame(focal_length=focal_length, rotation=cam.cam_rotation)
cam_coords = SkyCoord(x=cam.pix_x, y=cam.pix_y, frame=camera_frame)
tel_coords = cam_coords.transform_to(TelescopeFrame())

# Use only a subset of pixels, indicated by mask:

pixel_diameter = 2 * (
np.sqrt(cam.pix_area[0] / np.pi) / focal_length * u.rad

min_lambda = min_lambda.to_value(u.m)
max_lambda = max_lambda.to_value(u.m)

def negative_log_likelihood(
impact_parameter,
phi,
center_x,
center_y,
ring_width,
optical_efficiency_muon,
):
"""
Likelihood function to be called by minimizer

Parameters
----------
impact_parameter: float
Impact distance from telescope center
center_x: float
center of muon ring in the telescope frame
center_y: float
center of muon ring in the telescope frame
ring_width: float
Gaussian width of muon ring
optical_efficiency_muon: float
Efficiency of the optical system

Returns
-------
float: Likelihood that model matches data
"""

# Generate model prediction
prediction = image_prediction_no_units(
impact_parameter_m=impact_parameter,
oversampling=oversampling,
min_lambda_m=min_lambda,
max_lambda_m=max_lambda,
)

# scale prediction by optical efficiency of the telescope
prediction *= optical_efficiency_muon

return neg_log_likelihood_approx(image, prediction, spe_width, pedestal)

return negative_log_likelihood

geometry = telescope_description.camera.geometry
optics = telescope_description.optics

focal_length = optics.equivalent_focal_length.to_value(u.m)
pixel_area = geometry.pix_area[0].to_value(u.m**2)
pixel_radius = np.sqrt(pixel_area / np.pi) / focal_length

initial_guess = {}
initial_guess["phi"] = 0
initial_guess["optical_efficiency_muon"] = 0.1

return initial_guess

[docs]class MuonIntensityFitter(TelescopeComponent):
spe_width = FloatTelescopeParameter(
help="Width of a single photo electron distribution", default_value=0.5
).tag(config=True)

min_lambda_m = FloatTelescopeParameter(
help="Minimum wavelength for Cherenkov light in m", default_value=300e-9
).tag(config=True)

max_lambda_m = FloatTelescopeParameter(
help="Minimum wavelength for Cherenkov light in m", default_value=600e-9
).tag(config=True)

help="Hole radius of the reflector in m",
default_value=[
("type", "LST_*", 0.308),
("type", "MST_*", 0.244),
("type", "SST_1M_*", 0.130),
],
).tag(config=True)

oversampling = IntTelescopeParameter(
help="Oversampling for the line integration", default_value=3
).tag(config=True)

"""

Parameters
----------
center_x: Angle quantity
Initial guess for muon ring center in telescope frame
center_y: Angle quantity
Initial guess for muon ring center in telescope frame
Radius of muon ring from circle fitting
pixel_x: ndarray
X position of pixels in image from circle fitting
pixel_y: ndarray
Y position of pixel in image from circle fitting
image: ndarray
Amplitude of image pixels
mask marking the pixels to be used in the likelihood fit

Returns
-------
MuonEfficiencyContainer
"""
telescope = self.subarray.tel[tel_id]
if telescope.optics.n_mirrors != 1:
raise NotImplementedError(
"Currently only single mirror telescopes"
f" are supported in {self.__class__.__name__}"
)

negative_log_likelihood = build_negative_log_likelihood(
image,
telescope,
oversampling=self.oversampling.tel[tel_id],
min_lambda=self.min_lambda_m.tel[tel_id] * u.m,
max_lambda=self.max_lambda_m.tel[tel_id] * u.m,
spe_width=self.spe_width.tel[tel_id],
pedestal=pedestal,
)
negative_log_likelihood.errordef = Minuit.LIKELIHOOD

initial_guess = create_initial_guess(center_x, center_y, radius, telescope)

# Create Minuit object with first guesses at parameters
# strip away the units as Minuit doesnt like them

minuit = Minuit(negative_log_likelihood, **initial_guess)

minuit.print_level = 0

minuit.errors["impact_parameter"] = 0.5
minuit.errors["optical_efficiency_muon"] = 0.05

minuit.limits["impact_parameter"] = (0, None)
minuit.limits["phi"] = (-np.pi, np.pi)
minuit.limits["ring_width"] = (0.0, None)
minuit.limits["optical_efficiency_muon"] = (0.0, None)

minuit.fixed["center_x"] = True
minuit.fixed["center_y"] = True

# Perform minimisation