Add function to get point on shower axis in altaz

[maxnoe]

This can be used to compute the shower axis in telescope / camera frame, e.g. like this:

import astropy.units as u
from astropy.coordinates import AltAz
import matplotlib.pyplot as plt

from ctapipe.calib import CameraCalibrator
from ctapipe.io import EventSource
from ctapipe.coordinates import get_point_on_shower_axis, TelescopeFrame
from ctapipe.visualization import CameraDisplay


plt.style.use("dark_background")
plt.rcParams["savefig.facecolor"] = "0.15"

path = "dataset://lst_prod3_calibration_and_mcphotons.simtel.zst"

if __name__ == "__main__":
    with EventSource(path, focal_length_choice="EQUIVALENT") as source:
        it = iter(source)
        for _ in range(2):
            e = next(it)
        subarray = source.subarray

    calib = CameraCalibrator(subarray)
    calib(e)

    shower = e.simulation.shower
    distance = 5 * u.km

    fig, axs = plt.subplots(2, 2, layout="constrained", figsize=(8, 8))

    telescope_positions = source.subarray.tel_coords
    axis_point = get_point_on_shower_axis(
        core_x=shower.core_x,
        core_y=shower.core_y,
        alt=shower.alt,
        az=shower.az,
        telescope_position=telescope_positions,
        distance=distance,
    )

    for ax, (tel_id, dl1) in zip(axs.flat, e.dl1.tel.items()):
        tel_index = source.subarray.tel_ids_to_indices(tel_id)[0]

        pointing = AltAz(
            alt=e.pointing.tel[tel_id].altitude, az=e.pointing.tel[tel_id].azimuth
        )
        tel_frame = TelescopeFrame(telescope_pointing=pointing)

        geometry = source.subarray.tel[tel_id].camera.geometry
        disp = CameraDisplay(geometry.transform_to(tel_frame), image=dl1.image, ax=ax)

        ax.set_title(f"LST-{tel_id}")
        ax.set_axis_off()
        ax.set_ylabel("")
        ax.set_xlabel("")

        source_tel_frame = AltAz(alt=shower.alt, az=shower.az).transform_to(tel_frame)
        p = axis_point[tel_index].transform_to(tel_frame)

        disp.axes.plot(
            [source_tel_frame.fov_lon.deg, p.fov_lon.deg],
            [source_tel_frame.fov_lat.deg, p.fov_lat.deg],
            color="xkcd:baby blue",
            marker=".",
        )
        disp.overlay_coordinate(source_tel_frame, color="xkcd:light yellow")

    fig.suptitle(
        f"$E={{}}${shower.energy.to(u.GeV):.0f}, Impact = ({shower.core_x:.0f}, {shower.core_y:.0f})"
    )

    plt.savefig("shower_axis.png", dpi=200)

[moralejo]

Looks good to me. Can you show an example with larger images?

[maxnoe]

Looks good to me. Can you show an example with larger images?

Here is a 5 TeV shower:

[maxnoe]

This is at the moment relatively slow (using this for Algorithm 6 of the reconstruction comparison paper takes around 400 ms per event) but most time is spent in astropy coordinate transformations, which is directly affected by a patch made to improve performance there in the next astropy release.

[moralejo]

Looks good to me. Can you show an example with larger images?

Here is a 5 TeV shower:

Thanks, the lower-E one was not so convincing because the pixelation on small images makes the effect of a misalignment, but this one looks very convincing.

[ctao-dpps-sonarqube]

Analysis Details

0 Issues

• 0 Bugs
• 0 Vulnerabilities
• 0 Code Smells

Coverage and Duplications

• 100.00% Coverage (92.80% Estimated after merge)
• 0.00% Duplicated Code (0.70% Estimated after merge)

Project ID: cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs

View in SonarQube