import astropy.units as u
import numpy as np
from rta_reconstruction.utils.coordinates import pol2cart, sky_to_camera
[docs]
def disp(cog_x, cog_y, src_x, src_y, hillas_psi):
"""
Compute the disp parameters
Parameters
----------
cog_x: `numpy.ndarray` or float
cog_y: `numpy.ndarray` or float
src_x: `numpy.ndarray` or float
src_y: `numpy.ndarray` or float
hillas_psi: `numpy.ndarray` or float
Returns
-------
(disp_dx, disp_dy, disp_norm, disp_angle, disp_sign):
disp_dx: 'astropy.units.m`
disp_dy: 'astropy.units.m`
disp_norm: 'astropy.units.m`
disp_angle: 'astropy.units.rad`
disp_sign: `numpy.ndarray`
"""
disp_dx = src_x - cog_x
disp_dy = src_y - cog_y
disp_norm = disp_dx * np.cos(hillas_psi) + disp_dy * np.sin(hillas_psi)
disp_sign = np.sign(disp_norm)
disp_norm = np.abs(disp_norm)
# disp_sign : indicates in which direction, "positive" or "negative", we must move along the
# reconstructed image axis (with direction defined by the versor cos(hillas_psi), sin(hillas_psi))
# we must move from cog_x, cog_y to get closest to the true direction (src_x, src_y)
if hasattr(disp_dx, "__len__"):
disp_angle = np.arctan(disp_dy / disp_dx)
disp_angle[disp_dx == 0] = np.pi / 2.0 * np.sign(disp_dy[disp_dx == 0])
else:
if disp_dx == 0:
disp_angle = np.pi / 2.0 * np.sign(disp_dy)
else:
disp_angle = np.arctan(disp_dy / disp_dx)
return disp_dx, disp_dy, disp_norm, disp_angle, disp_sign
[docs]
def miss(disp_dx, disp_dy, hillas_psi):
"""
Compute miss
Parameters
----------
disp_dx: `numpy.ndarray` or float
disp_dy: `numpy.ndarray` or float
hillas_psi: `numpy.ndarray` or float
Returns
-------
"""
return np.abs(np.sin(hillas_psi) * disp_dx - np.cos(hillas_psi) * disp_dy)
[docs]
def disp_vector_pol2cart(disp_norm, disp_angle, disp_sign):
"""
Compute `disp_norm.dx` and `disp_norm.dy` vector from `disp_norm.norm`,
`disp_norm.angle` and `disp_norm.sign`
Parameters
----------
disp_norm: float
disp_angle: float
disp_sign: float
Returns
-------
disp_dx, disp_dy
"""
return np.transpose(pol2cart(disp_norm, disp_angle, disp_sign))
[docs]
def disp_to_pos(disp_dx, disp_dy, cog_x, cog_y):
"""
Calculates source position in camera coordinates(x,y)
from the reconstructed disp.
Parameters
----------
disp: DispContainer
cog_x: float
Coordinate x of the center of gravity of Hillas ellipse
cog_y: float
Coordinate y of the center of gravity of Hillas ellipse
Returns
-------
(source_pos_x, source_pos_y)
"""
source_pos_x = cog_x + disp_dx
source_pos_y = cog_y + disp_dy
return source_pos_x, source_pos_y
[docs]
def update_disp_with_effective_focal_length(data, effective_focal_length=29.30565 * u.m):
"""Update disp parameters using effective focal length
Parameters
----------
data: Pandas DataFrame
config: dictionnary containing configuration
"""
source_pos_in_camera = sky_to_camera(
u.Quantity(data["mc_alt"].values, u.rad, copy=False),
u.Quantity(data["mc_az"].values, u.rad, copy=False),
effective_focal_length,
u.Quantity(data["mc_alt_tel"].values, u.rad, copy=False),
u.Quantity(data["mc_az_tel"].values, u.rad, copy=False),
)
expected_src_pos_x_m = source_pos_in_camera.x.to_value(u.m)
expected_src_pos_y_m = source_pos_in_camera.y.to_value(u.m)
data["src_x"] = expected_src_pos_x_m
data["src_y"] = expected_src_pos_y_m
disp_dx, disp_dy, disp_norm, disp_angle, disp_sign = disp.disp(
data["x"].values, data["y"].values, expected_src_pos_x_m, expected_src_pos_y_m, data["psi"].values
)
data["disp_dx"] = disp_dx
data["disp_dy"] = disp_dy
data["disp_norm"] = disp_norm
data["disp_angle"] = disp_angle
data["disp_sign"] = disp_sign