<jupyter_start><jupyter_text>Fantastic coordinates and how to transform themLST Analysis Bootcamp - LNL, Legnaro, Padova - 26.11.2018<NAME> (@thomasgas)INFN, UNITS, UGA, LAPP  Intro- Short presentation on the different coordinates frames that we have in ctapipe at the moment.- More a **user point-of-view** presentation rather then a developer one.- Changes foreseen in ctapipe before the version 1.0, but these reference frames should be **ALMOST** ok. - Development of frames and transformation might change (faster computation, improved , etc..) **CREED_VTK** (CTA Rendering Event-by-Event Display with VTK)I developed a library for the 3D rendering of the telescopes. Follow the instructions on https://github.com/thomasgas/CREED_VTK to install it inside your ctapipe environment and to use this notebook.  inside the ctapipe environment- git clone https://github.com/thomasgas/CREED_VTK.git- cd CREED_VTK- python setup.py install **Still under development...feel free to add comments (issues) and suggestions for funtionalities!**  Some imports...<jupyter_code>import astropy.units as u
import copy
import numpy as np
import matplotlib.pyplot as plt

from ctapipe.io import event_source
from ctapipe.calib import CameraCalibrator
from ctapipe.utils import get_dataset, get_dataset_path

from ctapipe.reco import HillasReconstructor
from ctapipe.reco.hillas_intersection import HillasIntersection

from ctapipe.image import tailcuts_clean, hillas_parameters, HillasParameterizationError
from ctapipe.image.timing_parameters import timing_parameters

from ctapipe.visualization import ArrayDisplay, CameraDisplay
from ctapipe.coordinates import HorizonFrame, GroundFrame, TiltedGroundFrame
from ctapipe.coordinates import NominalFrame, CameraFrame, TelescopeFrame<jupyter_output><empty_output><jupyter_text>Open some data. Default is the test dataset in ctapipe<jupyter_code>try:
    DATA_PATH = "/home/thomas/Programs/astro/CTAPIPE_DAN/"
    filename = 'gamma_20deg_0deg_run100___cta-prod3_desert-2150m-Paranal-merged.simtel.gz'
    filename = DATA_PATH + filename
    source = event_source(filename)
    # select layout
    layout = np.loadtxt("utils_coords/CTA.prod3Sb.3HB9-NG.lis", usecols=0, dtype=int)
    layout = set(layout)
    source.allowed_tels = layout
except:
    filename = get_dataset_path("gamma_test.simtel.gz")
    source = event_source(filename)
    
source.max_events = 10
events = [copy.deepcopy(event) for event in source]
event = events[0]
# Calibration and find "great event"
cal = CameraCalibrator(None, None, r1_product='HESSIOR1Calibrator', extractor_product='NeighbourPeakIntegrator')
for event in events:
    cal.calibrate(event)
    
# Find "big" event
events_amplitude = []
for event in events:
    event_amplitude = 0
    for tel_id in event.r0.tels_with_data:
        if event.dl1.tel[tel_id].image is not None:
            event_amplitude += event.dl1.tel[tel_id].image[0].sum()
    events_amplitude.append(event_amplitude)
events_amplitude = np.array(events_amplitude)

mm = events_amplitude.argmax()
event = events[mm]
event.r0.tels_with_data
tel_id = 4
event.inst.subarray.tel[tel_id]<jupyter_output><empty_output><jupyter_text>HorizonFrameSpherical system to describe the position of an object in terms of altitude and azimuth.Pointing direction or source position in the sky is described as an HorizonFrame.<jupyter_code>array_pointing = HorizonFrame(alt=event.mcheader.run_array_direction[1],
                              az=event.mcheader.run_array_direction[0])
print(array_pointing)<jupyter_output><HorizonFrame Coordinate (pointing_direction=None, array_direction=None): (az, alt) in deg
    (0., 69.99999967)><jupyter_text>Optional (**present implementation**): pass "pointing direction" and "array direction" to transform directly from CameraFrame and HorizonFrame.  CameraFrameCamera coordinate frame. Simple physical cartesian frame describing the 2D dimensional position of objects in the focal plane of the telescope.**Typical usage**: Position of pixels in the focal plane.<jupyter_code>geom = event.inst.subarray.tel[tel_id].camera
pix_x, pix_y = geom.pix_x, geom.pix_y
fl = event.inst.subarray.tel[tel_id].optics.equivalent_focal_length

camera_coord = CameraFrame(x=pix_x, y=pix_y,
                           focal_length=fl, rotation = 0*u.deg)
print(camera_coord)
%matplotlib qt
plt.scatter(camera_coord.x, camera_coord.y)
plt.axis('square');<jupyter_output><empty_output><jupyter_text>TelescopeFrameTelescope coordinate frame. Cartesian system to describe the angular offset of a given position in reference to pointing direction of a given telescope. Pointing corrections should applied to the transformation between this frame and the camera frame.Not really used now...<jupyter_code>telescope_coords = camera_coord.transform_to(TelescopeFrame())
print(telescope_coords)
%matplotlib qt
plt.scatter(telescope_coords.x, telescope_coords.y)
plt.axis('square');<jupyter_output><empty_output><jupyter_text>NominalFrameNominal coordinate frame. Cartesian system to describe the angular offset of a given position in reference to the pointing direction of a nominal array pointing. This is identical, for most of the cases, to the TelescopeFrame, a part for **divergent pointing**.- 2D reconstruction (HillasIntersector) is performed in this frame - 3D reconstruction (HillasReconstructor) doesn't need this frame Let's play a bit with 3 divergent LSTs<jupyter_code>array_pointing = HorizonFrame(alt=45 * u.deg, az=0 * u.deg)

point_direction_1 = HorizonFrame(alt=45 * u.deg, az = -5 * u.deg)
point_direction_2 = HorizonFrame(alt=45 * u.deg, az = 5 * u.deg)
point_direction_3 = HorizonFrame(alt=50 * u.deg, az = 0 * u.deg)

nom_frame_1 = NominalFrame(array_direction=array_pointing, pointing_direction=point_direction_1)
nom_frame_2 = NominalFrame(array_direction=array_pointing, pointing_direction=point_direction_2)
nom_frame_3 = NominalFrame(array_direction=array_pointing, pointing_direction=point_direction_3)

nom_1 = camera_coord.transform_to(nom_frame_1)
nom_2 = camera_coord.transform_to(nom_frame_2)
nom_3 = camera_coord.transform_to(nom_frame_3)
%matplotlib qt
plt.scatter(y=nom_1.x, x=nom_1.y, c = "red", s = 10)
plt.scatter(y=nom_2.x, x=nom_2.y, c = "blue", s = 10)
plt.scatter(y=nom_3.x, x=nom_3.y, c = "green", s = 10)
plt.axis('square')
plt.show()
array_pointing = HorizonFrame(alt=45 * u.deg, az=0 * u.deg)

point_direction_1 = HorizonFrame(alt=70 * u.deg, az = -5 * u.deg)
point_direction_2 = HorizonFrame(alt=70 * u.deg, az = 5 * u.deg)
point_direction_3 = HorizonFrame(alt=72 * u.deg, az = 0 * u.deg)

nom_frame_1 = NominalFrame(array_direction=array_pointing, pointing_direction=point_direction_1)
nom_frame_2 = NominalFrame(array_direction=array_pointing, pointing_direction=point_direction_2)
nom_frame_3 = NominalFrame(array_direction=array_pointing, pointing_direction=point_direction_3)

nom_1 = camera_coord.transform_to(nom_frame_1).transform_to(array_pointing)
nom_2 = camera_coord.transform_to(nom_frame_2).transform_to(array_pointing)
nom_3 = camera_coord.transform_to(nom_frame_3).transform_to(array_pointing)
%matplotlib qt
from utils_coords.plot_on_sky import sky_in_box
az_list = [nom_1.az.value, nom_2.az.value, nom_3.az.value]
zen_list = [90-nom_1.alt.value, 90-nom_2.alt.value, 90-nom_3.alt.value]
labels = ["1", "2", "3"]
fig3 = plt.figure(figsize=(6, 6))
sky_in_box(fig3, az_deg=az_list, zen_deg=zen_list, label=labels, center=array_pointing)
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>GroundFrameGround coordinate frame. This is a simple cartesian frame describing the 3D position of objects compared to the array ground level in relation to the conter of the array. **Typical usage**: positions of telescopes on the ground (x, y, z)<jupyter_code>%matplotlib qt
event.inst.subarray.peek()
%matplotlib qt
event.inst.subarray.select_subarray("Prod3b layout", layout).peek()<jupyter_output>/home/thomas/Programs/astro/ctapipe/ctapipe/instrument/subarray.py:118: UserWarning: SubarrayDescription.pos_x is deprecated. Use tel_coords.x
  warnings.warn("SubarrayDescription.pos_x is deprecated. Use "
/home/thomas/Programs/astro/ctapipe/ctapipe/instrument/subarray.py:125: UserWarning: SubarrayDescription.pos_y is deprecated. Use tel_coords.y
  warnings.warn("SubarrayDescription.pos_y is deprecated. Use "<jupyter_text>Telescopes position in 3D<jupyter_code>from CREED_VTK import CREED_VTK 
try:
    render = CREED_VTK(event, telescopes_ids=list(layout))
except:
    render = CREED_VTK(event, telescopes_ids= event.inst.subarray.tel_ids.tolist())
render.add_gnd_tels()
render.add_gnd_frame(size=2000)
render.tel_labels()

render.camera_view(elev=20)
render.show(width= 1000, height=800)<jupyter_output><empty_output><jupyter_text>TiltedGroundFrameTilted ground coordinate frame. Cartesian system describing the 2D dimensional projected positions of objects in a tilted plane described by the pointing direction (HorizonFrame). **Tipical usage**: shower core position reconstruction.<jupyter_code>array_pointing = HorizonFrame(alt=event.mcheader.run_array_direction[1],
                              az=event.mcheader.run_array_direction[0])
try:
    render = CREED_VTK(event, telescopes_ids=list(layout))
except:
    render = CREED_VTK(event, telescopes_ids= event.inst.subarray.tel_ids.tolist())
render.add_gnd_tels()
render.add_gnd_frame(size=2000)
render.tel_labels() 

render.add_tilted_tels()
render.add_tilted_frame(size=2000)

render.camera_view(elev=20)
render.show(width= 1000, height=800)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09438508275717578, 'fraction_numerical': 0.019170333123821496, 'mean_word_length': 6.1674174174174174, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 37, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '118357', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3333, 'n_tokens_neox': 3125, 'n_words': 755, 'path': 'notebooks/coordinates_ctapipe_Bootcamp_2018.ipynb', 'repo_name': 'Mike7477/cta-lstchain', 'size': '9557'}