Use case 1: Shower tracks

import showermodel as sm

The ShowerModel package includes tools to make geometric calculations and to visualize the projection of linear shower tracks in the local coordinate system of a telescope. This notebook describes the use of the classes Track, Telescope and Projection.

Construction of a shower track

The default Track object corresponds to a vertical shower with impact point at the origin of coordinates. By default, the track is assumed to start at the top of the atmosphere at 112.83 km a.s.l and to end at ground level at 2.20 km a.s.l., where the atmosphere is discretized in 550 steps. However, the atmosphere extension and discritization as well as the track initial point can be modified when constructing a Track object.

# Default vertical track
# track = sm.Track()

# Track with impact point at 0.5 km east, 0.1 km north, zenith angle of 20 degrees and azimuth angle of 30 degrees
# track = sm.Track(x0=0.5, y0=0.1, theta=20., az=30.)  # alt=70. can be used instead of theta=20.

# Setting ground level and top of the atmosphere (km asl), and the number of discretization steps
track = sm.Track(theta=20., h0=1.8, h_top=100., N_steps=600)

# Setting track initial point. Note that theta may be greater than 90 degrees (ascending track). If so, x0, y0 can
# be given instead of xi, yi
# track = sm.Track(xi=0.5, yi=0.1, zi=14., theta=160.)

The generated Track object is a subclass of DataFrame and stores the (x, y, z) coordinates in km of each track point and the travel time in microseconds assuming the speed of light.

track

	x	y	z	t
0	0.0	0.029785	0.081833	348.282762
1	0.0	0.089355	0.245500	347.701807
2	0.0	0.148924	0.409167	347.120851
3	0.0	0.208494	0.572833	346.539896
4	0.0	0.268064	0.736500	345.958940
...	...	...	...	...
595	0.0	35.473813	97.463500	2.614299
596	0.0	35.533383	97.627167	2.033344
597	0.0	35.592953	97.790833	1.452388
598	0.0	35.652522	97.954500	0.871433
599	0.0	35.712092	98.118167	0.290478

600 rows × 4 columns

A Track object has also some useful attributes and methods.

track = sm.Track(x0=0.5, y0=0.1, theta=20., az=120., h0=1.8)

print("Unit vector parallel to shower axis (upwards):", (track.ux, track.uy, track.uz))
print("Coordinates of ground impact point:", (track.x0, track.y0, 0.))
print("Coordinates of shower point at 3.5 km a.s.l.:", track.h_to_xyz(3.5))
print("Travel time in us of shower front at 27 km above observation level:", track.z_to_t(27.))

Unit vector parallel to shower axis (upwards): (0.29619813272602386, -0.17101007166283427, 0.9396926207859084)
Coordinates of ground impact point: (0.5, 0.1, 0.0)
Coordinates of shower point at 3.5 km a.s.l.: (1.0358526974630422, -0.20937469912627193, 1.7)
Travel time in us of shower front at 27 km above observation level: 298.272204399376

Help on Track.

# track?

Calculation of the shower track projection relative to a telescope position

The default Telescope object corresponds to a telescope positioned at the origin of coordinates, pointing at zenith direction and with an angular aperture of 10 degrees. These parameters can be modified by the user (see also UC3).

# Telescope positioned at 0.1 km east and 0.2 km south
# telescope = sm.Telescope(x=0.1, y=-0.2)

# Setting the pointing direction and aperture
telescope = sm.Telescope(x=0.1, y=-0.2, z=0.1, theta=20., az=120., apert=12.)

Help on Telescope.

# sm.Telescope?

The Projection constructor calculates the coordinates (alt/az and theta/phi) of a shower track relative to a telescope position.

projection = sm.Projection(telescope, track)

Equivalent methods of Track and Telescope objects are also available.

# Equivalent method of Track objects
# projection = track.Projection(telescope)

# Equivalent method of Telescope object
# projection = telescope.Projection(track)

The generated Projection object is a subclass of DataFrame.

projection

	distance	alt	az	theta	phi	beta	time	FoV
0	0.515540	0.103984	56.887513	81.003132	306.074599	81.003132	1.447039	False
1	0.588703	20.151143	63.697352	59.876772	306.074599	59.876772	0.974510	False
2	0.720876	34.150606	69.578691	44.939486	306.074599	44.939486	0.698813	False
3	0.886033	43.201629	74.608921	35.078197	306.074599	35.078197	0.533135	False
4	1.068993	49.134256	78.897894	28.446219	306.074599	28.446219	0.426844	False
...	...	...	...	...	...	...	...	...
545	117.162399	69.891810	119.345934	0.249013	306.074599	0.249013	0.000030	True
546	117.377224	69.892011	119.347125	0.248557	306.074599	0.248557	0.000023	True
547	117.592049	69.892210	119.348311	0.248103	306.074599	0.248103	0.000017	True
548	117.806874	69.892409	119.349493	0.247651	306.074599	0.247651	0.000010	True
549	118.021699	69.892607	119.350671	0.247200	306.074599	0.247200	0.000003	True

550 rows × 8 columns

Note 1: beta is the angular distance (in degrees) of each shower point relative to the apparent position of the cosmic-ray source (at infinity).

Note 2: time is the arrival time (in microseconds) of light emitted from each shower point and reaching the telescope.

There are some useful attributes and methods of both Telescope and Projection objects.

print("Angular aperture in degrees of the telescope:", telescope.apert)
print("Unit vector parallel to the telescope pointing direction:", (telescope.ux, telescope.uy, telescope.uz))
print("Position angle (from north) of the right-hand direction from the telescope point of view:", telescope.phi_right)
print("FoV coordinates theta/phi corresponding to alt=72., az=120. relative to telescope:",
      telescope.altaz_to_thetaphi(72., 120.))  # Equivalent to projection.altaz_to_thetaphi(72., 120.)
print("theta/phi coordinates of the apparent position of the cosmic-ray source:", (projection.theta_inf, projection.phi_inf))
print("Angular distance of the ground interaction point relative to the apparent position of the cosmic-ray source:",
      projection.beta_0)

Angular aperture in degrees of the telescope: 12.0
Unit vector parallel to the telescope pointing direction: (0.29619813272602386, -0.17101007166283427, 0.9396926207859084)
Position angle (from north) of the right-hand direction from the telescope point of view: 151.51876171866053
FoV coordinates theta/phi corresponding to alt=72., az=120. relative to telescope: (1.9999999999999472, 61.51876171866053)
theta/phi coordinates of the apparent position of the cosmic-ray source: (0.0, 255.55500518658698)
Angular distance of the ground interaction point relative to the apparent position of the cosmic-ray source: 93.01201998391794

Help on Projection.

# sm.Projection?

Shower track visualtization

Projection objects have also a method to visualize the shower track in both zenith projection or the telescope field of view projection. The generated plots can be stored in PolarAxesSubplot objects for later modification.

# Just show the plots
projection.show();

# Store the plots in PolarAxesSubplot objects
# ax1, ax2 = projection.show()

# Equivalent method of Track objects
# projection, (ax1, ax2) = track.show_projection(telescope)

There are some plot options.

# Maximum offset angle in the FoV plot
projection.show(max_theta=6.);  # telescope.apert = 12.

# No axes and blue mark at slanth depth of 400 g/cm2 (it depends on atmosphere model, see UC2)
projection.show(axes=False, X_mark=400.);

Help on projection.show method.

# projection.show?

Example of track ‘search’

These tools can be used to determine the telescope pointing direction to observe different parts of a shower track.

track = sm.Track(x0=0.54, y0=-0.56, theta=120., az=47.) # Ascending shower starting at ground level (e.g., neutrino shower)
telescope = sm.Telescope(x=-0.2, y=-0.1, z=0.35)  # Telescope placed on the top of a hill
projection = sm.Projection(telescope, track)
projection.query('(alt < 20.)')  # Lower part of the shower

	distance	alt	az	theta	phi	beta	time	FoV
0	0.878923	-16.486430	133.374960	106.486430	226.625040	78.790391	0.470574	False
1	0.892676	-3.100588	156.333694	93.100588	203.666306	105.025577	1.858302	False
2	1.070003	8.213345	174.418306	81.786655	185.581694	126.317190	3.791641	False
3	1.347818	15.227347	186.702718	74.772653	173.297282	140.232964	6.060163	False
4	1.676899	19.332899	194.889401	70.667101	165.110599	149.060236	8.499688	False

# Setting the new pointing
telescope = sm.Telescope(x=-0.2, y=-0.1, z=0.35, alt=-3., az=156.)