Reconstruction of EAS direction

Direction reconstruction

This module contains two classes that can be used to reconstruct HiSPARC events and coincidences. The classes know how to extract the relevant information from the station and event or cluster and coincidence. Various algorithms which do the reconstruction are also defined here. The algorithms require positions and arrival times to do the reconstruction.

Each algorithm has a reconstruct_common() method which always requires arrival times, x, and y positions and optionally z positions and previous reconstruction results. The data is then prepared for the algorithm and passed to the reconstruct() method which returns the reconstructed theta and phi coordinates.

class sapphire.analysis.direction_reconstruction.EventDirectionReconstruction(station)

Reconstruct direction for station events

This class is aware of ‘events’ and ‘stations’. Initialize this class with a ‘station’ and you can reconstruct events using reconstruct_event(). To use other algorithms overwrite the direct and fit attributes.

Parameters:stationsapphire.clusters.Station object.
reconstruct_event(event, detector_ids=None, offsets=[0.0, 0.0, 0.0, 0.0], initial=None)

Reconstruct a single event

Parameters:
  • event – an event (e.g. from an events table), or any dictionary-like object containing the keys necessary for reconstructing the direction of a shower (e.g. arrival times).
  • detector_ids – list of the detectors to use for reconstruction. The detector ids are 0-based, unlike the column names in the esd data.
  • offsets – time offsets for each detector or a Station object.
  • initial – dictionary with already fitted shower parameters.
Returns:

theta, phi, and detector ids.

reconstruct_events(events, detector_ids=None, offsets=[0.0, 0.0, 0.0, 0.0], progress=True, initials=None)

Reconstruct events

Parameters:
  • events – the events table for the station from an ESD data file.
  • detector_ids – detectors to use for the reconstructions.
  • offsets – time offsets for each detector or a Station object.
  • progress – if True show a progress bar while reconstructing.
  • initials – list of dictionaries with already reconstructed shower parameters.
Returns:

list of theta, phi, and detector ids.

class sapphire.analysis.direction_reconstruction.CoincidenceDirectionReconstruction(cluster)

Reconstruct direction for coincidences

This class is aware of ‘coincidences’ and ‘clusters’. Initialize this class with a ‘cluster’ and you can reconstruct a coincidence using reconstruct_coincidence(). To use other algorithms overwrite the direct,``fit``, and curved attributes.

Parameters:clusterBaseCluster object.
reconstruct_coincidence(coincidence_events, station_numbers=None, offsets=None, initial=None)

Reconstruct a single coincidence

Parameters:
  • coincidence_events – a coincidence list consisting of three or more (station_number, event) tuples.
  • station_numbers – list of station numbers, to only use events from those stations.
  • offsets – a dictionary of either lists of detector timing offsets or Station objects for each station.
  • initial – dictionary with already fitted shower parameters.
Returns:

list of theta, phi, and station numbers.

reconstruct_coincidences(coincidences, station_numbers=None, offsets=None, progress=True, initials=None)

Reconstruct all coincidences

Parameters:
  • coincidences – a list of coincidence events, each consisting of three or more (station_number, event) tuples.
  • station_numbers – list of station numbers, to only use events from those stations.
  • offsets – dictionary with detector offsets for each station. These detector offsets should be relative to one detector from a specific station.
  • progress – if True show a progress bar while reconstructing.
  • initials – list of dictionaries with already reconstructed shower parameters.
Returns:

list of theta, phi, and station numbers.

get_station_offsets(coincidence_events, station_numbers, offsets, ts0)
determine_best_offsets(station_numbers, midnight_ts, offsets)

Determine best combined station and detector offsets

Check which station is best used as reference. Allow offsets via other stations, intermediate stations are used if it reduces the offset error.

Parameters:
  • station_numbers – list of stations in the coincidence or also other stations can are allow to be the reference station.
  • midnight_ts – timestamp of midnight before the coincidence.
  • offsets – a dictionary of Station objects for each station.
Returns:

combined detector and station offsets for given station, relative to the reference station.

determine_best_reference(error_matrix, station_numbers)
class sapphire.analysis.direction_reconstruction.CoincidenceDirectionReconstructionDetectors(cluster)

Reconstruct direction for coincidences using each detector

Instead of only the first arrival time per station this class uses the arrival time in each detector for the reconstruction.

reconstruct_coincidence(coincidence_events, station_numbers=None, offsets=None, initial=None)

Reconstruct a single coincidence

Parameters:
  • coincidence_events – a coincidence list consisting of one or more (station_number, event) tuples.
  • station_numbers – list of station numbers, to only use events from those stations.
  • offsets – dictionary with detector offsets for each station. These detector offsets should be relative to one detector from a specific station.
  • initial – dictionary with already fitted shower parameters.
Returns:

list of theta, phi, and station numbers.

class sapphire.analysis.direction_reconstruction.BaseDirectionAlgorithm

No actual direction reconstruction algorithm

Simply returns (nan, nan) as direction.

classmethod reconstruct_common(t, x, y, z=None, initial=None)

Reconstruct shower angles

Parameters:
  • t – detector arrival time in ns.
  • x,y – positions of detectors in m.
  • z – height of detectors in m.
  • initial – dictionary containing values from previous reconstructions.
Returns:

reconstructed theta and phi angles.

static reconstruct()

Reconstruct shower angles

Returns:reconstructed theta and phi angles.
class sapphire.analysis.direction_reconstruction.DirectAlgorithm

Reconstruct angles using direct analytical formula.

This implements the equations derived in Fokkema2012 sec 4.2. (DOI: 10.3990/1.9789036534383)

This algorithm assumes each detector is at the same altitude.

Note! The detectors are 0-based.

classmethod reconstruct_common(t, x, y, z=None, initial=None)

Reconstruct angles from 3 detections

This function converts the coordinates to be suitable for the algorithm.

Parameters:
  • t – arrival times in detector 0, 1 and 2 in ns.
  • x,y – positions of detector 0, 1 and 2 in m.
  • z – height of detectors 0, 1 and 2 is ignored.
  • initial – dictionary containing values from previous reconstructions is ignored.
classmethod reconstruct(dt1, dt2, r1, r2, phi1, phi2)

Reconstruct angles from 3 detections

Parameters:
  • dt# – arrival times in detector 1 and 2 relative to detector 0 in ns (!).
  • r#,phi# – position of detector 1 and 2 relative to detector 0 in m and radians.
Returns:

theta as given by Fokkema2012 eq 4.14, phi as given by Fokkema2012 eq 4.13.

classmethod rel_theta1_errorsq(theta, phi, phi1, phi2, r1=10, r2=10)

Fokkema2012, eq 4.23

classmethod rel_theta2_errorsq(theta, phi, phi1, phi2, r1=10, r2=10)

Fokkema2012, eq 4.23

static rel_phi_errorsq(theta, phi, phi1, phi2, r1=10, r2=10)

Fokkema2012, eq 4.22

classmethod dphi_dt0(theta, phi, phi1, phi2, r1=10, r2=10)

Fokkema2012, eq 4.19

static dphi_dt1(theta, phi, phi1, phi2, r1=10, r2=10)

Fokkema2012, eq 4.20

static dphi_dt2(theta, phi, phi1, phi2, r1=10, r2=10)

Fokkema2012, eq 4.21

class sapphire.analysis.direction_reconstruction.DirectAlgorithmCartesian

Reconstruct angles using direct analytical formula.

This implements the equations derived in Montanus2014. “Direction reconstruction of cosmic air showers with detectorstations at different altitudes”

Here the 2D version is used, assuming each detector is at the same altitude.

classmethod reconstruct_common(t, x, y, z=None, initial=None)

Reconstruct angles from 3 detections

This function converts the coordinates to be suitable for the algorithm.

Parameters:
  • t – arrival times in detector 0, 1 and 2 in ns.
  • x,y – positions of detector 0, 1 and 2 in m.
  • z – height of detectors 0, 1 and 2 is ignored.
  • initial – dictionary containing values from previous reconstructions is ignored.
static reconstruct(dt1, dt2, dx1, dx2, dy1, dy2)

Reconstruct angles from 3 detections

Parameters:
  • dt# – arrival times in detector 1 and 2 relative to detector 0 in ns.
  • dx#,dy# – position of detector 1 and 2 relative to detector 0 in m.
Returns:

theta as given by Montanus2014 eq 27, phi as given by Montanus2014 eq 26.

class sapphire.analysis.direction_reconstruction.DirectAlgorithmCartesian3D

Reconstruct angles using direct analytical formula.

This implements the equations derived in Montanus2014. “Direction reconstruction of cosmic air showers with detectorstations at different altitudes”

Here the 3D version is used, assuming each detector is at the same altitude.

classmethod reconstruct_common(t, x, y, z=None, initial=None)

Reconstruct angles from 3 detections

This function converts the coordinates to be suitable for the algorithm.

Parameters:
  • t – arrival times in detector 0, 1 and 2 in ns.
  • x,y,z – positions of detector 0, 1 and 2 in m.
  • initial – dictionary containing values from previous reconstructions is ignored.
static reconstruct(dt1, dt2, dx1, dx2, dy1, dy2, dz1=0, dz2=0)

Reconstruct angles from 3 detections

Parameters:
  • dt# – arrival times in detector 1 and 2 relative to detector 0 in ns.
  • dx#,dy#,dz# – position of detector 1 and 2 relative to detector 0 in m.
Returns:

theta as given by Montanus2014 eq 24, phi as given by Montanus2014 eq 22.

class sapphire.analysis.direction_reconstruction.SphereAlgorithm

Reconstruct the direction in equatorial coordinates

Note: currently incompatible with the other algorithms!

This class uses a different coordinate systems than the other algorithms. The location input is in ECEF coordinates and a timestamp is required to connect the direction to the equatorial coordinates.

classmethod reconstruct_equatorial(t, x, y, z, timestamp)

Reconstructs the source in the Equatorial Coordinate System.

Parameters:
  • t – An array with three arrival times in ns.
  • x,y,z – arrays with the ECEF locations of the three detectors / stations in meters.
  • timestamp – The UTC timestamp of the coincidence in s.
Returns:

declination and right ascension of the source. The apparent location of the cosmic ray source in the Equatorial Coordinate System.

static interaction_curve(x, y, z, t, t_int)

Calculates the curve of possible primary interactions

This uses the arrival times in three detectors. The algorithm is based on location calculations used for LORAN, DECCA, RACAL, GPS as described by N.G. Schultheiss 2012

Parameters:
  • x,y,z – arrays with the orthogonal coordinates of the three detection points in m.
  • t – arrival times of the detectors in ns.
  • t_int – interaction time in ns.
Returns:

parameters x_int, y_int, z_int

class sapphire.analysis.direction_reconstruction.FitAlgorithm3D
classmethod reconstruct_common(t, x, y, z=None, initial=None)

Reconstruct angles from 3 or more detections

This function converts the arguments to be suitable for the algorithm.

Parameters:
  • t – arrival times of the detectors in ns.
  • x,y,z – positions of the detectors in m. The height for all detectors will be set to 0 if not given.
  • initial – dictionary containing values from previous reconstructions is ignored.
classmethod reconstruct(t, x, y, z)

Reconstruct angles for many detections

Parameters:
  • t# – arrival times in the detectors in ns.
  • x#,y#,z# – position of the detectors in m.
Returns:

theta as given by Montanus2014 eq 21, phi as given by Montanus2014 eq 22.

static constraint_normal_vector(n)

This should be equal to zero

static best_fit(n_xyz, dt, dx, dy, dz)

The function to be minimized to find the direction

Parameters:
  • n_xyz – list containing the unit vector.
  • dt – list of relative arrival times in the detectors in ns.
  • dx,dy,dz – list of relative detector positions in m.
Returns:

least sum of squares as in Montanus2014, eq 36

class sapphire.analysis.direction_reconstruction.RegressionAlgorithm

Reconstruct angles using an analytical regression formula.

This implements the equations as for ISVHECRI (Montanus 2014). “Direction reconstruction of cosmic air showers with three or more detectorstations in a horizontal (for the moment) plane”

classmethod reconstruct_common(t, x, y, z=None, initial=None)

Reconstruct angles from 3 or more detections

This function converts the arguments to be suitable for the algorithm.

Parameters:
  • t – arrival times of the detectors in ns.
  • x,y,z – positions of the detectors in m. The height is ignored.
  • initial – dictionary containing values from previous reconstructions is ignored.
classmethod reconstruct(t, x, y)

Reconstruct angles for many detections

Parameters:
  • t – arrival times in the detectors in ns.
  • x,y – positions of the detectors in m.
Returns:

theta as derived by Montanus2014, phi as derived by Montanus2014.

class sapphire.analysis.direction_reconstruction.RegressionAlgorithm3D

Reconstruct angles by iteratively applying a regression formula.

This implements the equations as recently derived (Montanus 2014). “Direction reconstruction of cosmic air showers with three or more detectorstations at arbitrary altitudes”

MAX_ITERATIONS = 1000
classmethod reconstruct_common(t, x, y, z=None, initial=None)

Reconstruct angles from 3 or more detections

This function converts the arguments to be suitable for the algorithm.

Parameters:
  • t – arrival times of the detectors in ns.
  • x,y,z – positions of the detectors in m. The height for all detectors will be set to 0 if not given.
  • initial – dictionary containing values from previous reconstructions is ignored.
classmethod reconstruct(t, x, y, z)

Reconstruct angles for many detections

Parameters:
  • t – arrival times in the detectors in ns.
  • x,y,z – positions of the detectors in m.
Returns:

theta as derived by Montanus2014, phi as derived by Montanus2014.

class sapphire.analysis.direction_reconstruction.CurvedMixin

Provide methods to estimate the time delay due to front curvature

Given a core location, detector position, and shower angle the radial core distance can be determined, which can be used to determine the expected time delay.

time_delay(x, y, core_x, core_y, theta, phi)
classmethod radial_core_distance(x, y, core_x, core_y, theta, phi)

Determine the radial core distance

Parameters:
  • x,y,z – positions of the detectors in m.
  • core_x,core_y – core position at z = 0 in m.
  • theta,phi – reconstructed shower direction.
Returns:

radial core distance in m.

class sapphire.analysis.direction_reconstruction.CurvedRegressionAlgorithm

Reconstruct angles taking the shower front curvature into account.

Take the shower front curvature into account. Assumes knowledge about the shower core position.

MAX_ITERATIONS = 1000
reconstruct_common(t, x, y, z=None, initial=None)

Reconstruct angles from 3 or more detections

This function converts the arguments to be suitable for the algorithm.

Parameters:
  • t – arrival times of the detectors in ns.
  • x,y,z – positions of the detectors in m. The height is ignored.
  • initial – dictionary containing values from previous reconstructions, including core position.
reconstruct(t, x, y, core_x, core_y)

Reconstruct angles for many detections

Parameters:
  • t – arrival times in the detectors in ns.
  • x,y – positions of the detectors in m.
  • core_x,core_y – core position at z = 0 in m.
Returns:

theta as derived by Montanus2014, phi as derived by Montanus2014.

class sapphire.analysis.direction_reconstruction.CurvedRegressionAlgorithm3D

Reconstruct angles accounting for front curvature and detector altitudes

Take the shower front curvature and different detector heights into account. Assumes knowledge about the shower core position.

MAX_ITERATIONS = 1000
reconstruct_common(t, x, y, z=None, initial=None)

Reconstruct angles from 3 or more detections

This function converts the arguments to be suitable for the algorithm.

Parameters:
  • t – arrival times of the detectors in ns.
  • x,y,z – positions of the detectors in m. The height for all detectors will be set to 0 if not given.
  • initial – dictionary containing values from previous reconstructions, including core position.
reconstruct(t, x, y, z, core_x, core_y)

Reconstruct angles for many detections

Parameters:
  • t – arrival times in the detectors in ns.
  • x,y,z – positions of the detectors in m.
  • core_x,core_y – core position at z = 0 in m.
Returns:

theta as derived by Montanus2014, phi as derived by Montanus2014.

sapphire.analysis.direction_reconstruction.logic_checks(t, x, y, z)

Check for impossible reconstructions

Criteria:

  • No two detectors are at the same position.
  • Time difference between two detections should be less than distance / c.
  • All detectors on a line is bad.

To fix:

  • Time difference can still be to large in cases where a different distance becomes relevant.
Parameters:
  • t – arrival times in the detectors in ns.
  • x,y,z – positions of the detectors in m.
Returns:

True if the checks pass, False otherwise.

sapphire.analysis.direction_reconstruction.warning_only_three()