Reconstruction

grand.analysis reconstructs arrival direction, the distance to Xmax and an electromagnetic-energy proxy from recorded times and amplitudes. It needs the optional iminuit dependency (pip install -e ".[analysis]"; it is in the conda environment). Notebook 11 (see Tutorial notebooks) runs the chain step by step, on showers with a known answer and on ten GP13 candidates.

Note

Each fit’s tests show that it recovers what its own forward model generates; agreement with full simulations and measured showers is not yet tested.

grand.analysis.fitting.plane_wave	Plane-wave fit of the arrival direction from antenna times.
grand.analysis.fitting.spherical	Spherical-wave fit of the arrival direction and source distance.
grand.analysis.fitting.adf	Angular distribution function (ADF) fit of the antenna amplitudes.
grand.analysis.energy_reco.voltage	Electromagnetic-energy proxy from the fitted radio amplitude.
grand.analysis.signals.extraction	Peak amplitude and peak time extraction from antenna traces.
grand.analysis.cramer_rao_bounds.cramer_rao	Cramér-Rao lower bounds on the parameters of the reconstruction fits.
grand.analysis.physics.cherenkov_angle	Cherenkov angle from a two-emission-point model around the source.
grand.analysis.physics.atmosphere	Refractive index of the atmosphere, locally and averaged along a path.
grand.analysis.coords.array_shower	Shower direction vectors and the rotation into the shower frame.
grand.analysis.geom.angles	Angles between the shower axis, the antennas and the geomagnetic field.
grand.analysis.geom.footprint	Shower core and Cherenkov footprint on the ground.
grand.analysis.constants	Physical constants and site defaults used by the reconstruction.

grand.analysis.fitting.plane_wave

Plane-wave fit of the arrival direction from antenna times.

grand.analysis.fitting.plane_wave.PWF_semianalytical(Xants, tants, verbose=False, c=299792458.0, n=1.000136, sigma=None)

Solve the minimization problem using a semi-analytical approach.

(see section 2.1.2)

Parameters:

• (ndarray) (tants)

• (ndarray)

• (bool) (verbose)

• (float) (c)

• ndarray) (n (float or)

• ndarray (sigma (float or) – seconds; weights the antennas.

• optional) (Timing uncertainty of each antenna, in) – seconds; weights the antennas.

Returns:

• ndarray (Theta in [0, pi] and phi in [0, 2*pi), in radians: the direction)

• the shower comes from.

See also

grand.analysis.fitting.spherical.recons_swf

The spherical-wave fit, which also gives the distance to the source.

grand.analysis.fitting.adf.recons_ADF

The amplitude fit.

Examples

Recover the direction of a plane wave from the times it reaches 30 antennas. Positions are in meters, times in seconds, angles in radians:

import numpy as np
from grand.analysis.fitting.plane_wave import PWF_model, PWF_semianalytical

rng = np.random.default_rng(0)
xants = np.column_stack([rng.uniform(-2000, 2000, 30),
                         rng.uniform(-2000, 2000, 30),
                         np.full(30, 1264.0)])
tants = PWF_model((np.radians(70.0), np.radians(30.0)), xants)

theta, phi = PWF_semianalytical(xants, tants)
print("zenith %.2f deg, azimuth %.2f deg" % (np.degrees(theta), np.degrees(phi)))

zenith 70.00 deg, azimuth 30.00 deg

grand.analysis.fitting.plane_wave.mean(X: ndarray, sigma=None)

Return the mean of X along its first axis, weighted by sigma.

Parameters:

• X (ndarray) – Values, one row per antenna.

• sigma (ndarray, optional) – Uncertainties: a vector weights each row by 1/sigma, a covariance matrix by the column sums of its inverse. Unweighted if omitted.

Returns:

The (weighted) mean.

Return type:

ndarray

grand.analysis.fitting.plane_wave.PWF_loss(params, Xants, tants, verbose=False, c=299792458.0, n=1.000136, sigma=None)

Define Chi2 by summing model residuals over individual antennas.

After maximizing likelihood over reference time.

grand.analysis.fitting.plane_wave.PWF_residuals(params, Xants, tants, verbose=False, c=299792458.0, n=1.000136)

Compute timing residuals for each antenna using the plane wave model.

Note that this is defined at up to an additive constant, that when minimizing the loss over it, amounts to centering the residuals.

grand.analysis.fitting.plane_wave.PWF_model(params, Xants, c=299792458.0, n=1.000136, groundAltitude=1231)

Generate plane wavefront timings.

Parameters:

• params (sequence) – (theta, phi) in radians.

• Xants (ndarray, shape (N, 3)) – Antenna positions: x North, y West, z the height above sea level (m).

• c (float, optional) – Speed of light and refractive index.

• n (float, optional) – Speed of light and refractive index.

• groundAltitude (float, optional) – Height above sea level of the frame’s origin, where the source distance is measured from (meters). The default, 1231 m, is the GP13 site; for simulation files pass the ground altitude grand.analysis.geom.antenna_positions_from_run() returns. For the plane wave it only shifts all times by one constant.

Returns:

Arrival times in seconds, relative to the origin’s.

Return type:

ndarray, shape (N,)

grand.analysis.fitting.spherical

Spherical-wave fit of the arrival direction and source distance.

grand.analysis.fitting.spherical.SWF_loss(theta, phi, r_xmax, t_s, Xants, tants, sigma=None, cr=299792458.0)

Define Chi2 by summing model residuals over antennas (i).

The loss is:

loss = sum_i ( cr (tants[i] - t_s) - n_i |Xants[i] - X_s| )**2
X_s  = r_xmax * (sin(theta) cos(phi), sin(theta) sin(phi), cos(theta))
       + (0, 0, groundAltitude)

where Xants are the antenna positions (shape (N, 3)), tants the trigger times (shape (N,)) and n_i the mean refractive index along the path. The source lies at -r_xmax * K, on the side the shower comes from (K is the propagation direction).

Parameters:

• theta (float) – Polar angle of the shower direction (radians).

• phi (float) – Azimuth angle of the shower direction (radians).

• r_xmax (float) – Distance from the ground to the emission point (meters).

• t_s (float) – Emission time of the source (seconds).

• Xants (np.ndarray) – Antenna positions in the detector reference frame, shape (N, 3).

• tants (np.ndarray) – Measured signal times at antennas (seconds), shape (N,).

• sigma (float, ndarray of shape (N,) or (N, N), optional) – Timing uncertainty (seconds): one for all antennas, one per antenna, or their covariance matrix (seconds squared). If provided, each residual is divided by its uncertainty.

• cr (float, optional) – Propagation speed of the signal (default: speed of light).

Returns:

Chi-square value (normalized if sigma is provided).

Return type:

float

grand.analysis.fitting.spherical.recons_swf(theta_pwf, phi_pwf, tants, Xants, sigma=None, maxiter=1000, seed=42)

Perform a SWF reconstruction using differential evolution minimization.

The minimization is performed around PWF direction to improve convergence.

Parameters:

• theta_pwf (float) – Initial zenith angle from Plane Wave Fit (radians).

• phi_pwf (float) – Initial azimuth angle from Plane Wave Fit (radians).

• tants (np.ndarray) – Measured antenna times (seconds).

• Xants (np.ndarray) – Antenna positions, shape (N, 3), in meters: x North, y West, z above sea level.

• sigma (float, ndarray of shape (N,) or (N, N), optional) – Timing uncertainty (seconds), as for SWF_loss(). A single value scales the chi2 without moving its minimum; per-antenna values weight the antennas.

• maxiter (int, optional) – Maximum number of iterations for differential evolution.

• seed (int, optional) – Random seed for reproducibility.

Returns:

(theta_swf, phi_swf, r_xmax_swf, t_s_swf): the direction the shower comes from, in radians; the distance from the source to (0, 0, groundAltitude), in meters; and the emission time, in seconds.

Return type:

tuple

grand.analysis.fitting.spherical.compute_Xsource_cartesian_coords(theta_swf, phi_swf, r_xmax, groundAltitude=1231)

Compute the Cartesian coordinates of the emission point (Xsource).

From spherical coordinates.

Parameters:

• theta_swf (float) – Polar angle (radians).

• phi_swf (float) – Azimuth angle (radians).

• r_xmax (float) – Distance to the source (meters).

• groundAltitude (float, optional) – Height above sea level of the frame’s origin, where the source distance is measured from (meters). The default, 1231 m, is the GP13 site; for simulation files pass the ground altitude grand.analysis.geom.antenna_positions_from_run() returns.

Returns:

Cartesian coordinates of the source in GRAND detector frame, shape (1, 3).

Return type:

np.ndarray

grand.analysis.fitting.spherical.SWF_model(theta, phi, r_xsource, t_s, Xants, groundAltitude=1231, cr=299792458.0)

Compute the expected arrival times at each antenna based on the Spherical Wave Front (SWF) model.

The emission point (Xsource) is located at a distance r_xsource along the direction opposite to the shower propagation, referenced from the ground altitude. Arrival times are computed assuming straight-line propagation with an average refraction index.

Parameters:

• theta (float) – Polar angle of the shower direction (radians).

• phi (float) – Azimuth angle of the shower direction (radians).

• r_xsource (float) – Distance from the shower maximum (emission point) to the reference point at the ground (meters).

• t_s (float) – Emission time of the source (seconds).

• Xants (np.ndarray) – Antenna positions, shape (N, 3): x North, y West, z the height above sea level, in meters.

• groundAltitude (float, optional) – Height above sea level of the frame’s origin, where the source distance is measured from (meters). The default, 1231 m, is the GP13 site; for simulation files pass the ground altitude grand.analysis.geom.antenna_positions_from_run() returns.

• cr (float, optional) – Propagation speed of the signal (default: speed of light).

Returns:

Expected arrival times at each antenna, shape (N,).

Return type:

np.ndarray

grand.analysis.fitting.adf

Angular distribution function (ADF) fit of the antenna amplitudes.

grand.analysis.fitting.adf.ADF_parameters(theta, phi, delta_omega, amplitude, Xants, Xsource, groundAltitude=1231, Bvec=None)

Compute all geometric parameters for the ADF function.

Inputs:

theta, phi : shower direction angles (rad) (from ADF_recons, best fit values) delta_omega : ADF shape parameter (output from ADF_recons, best fit values) amplitude : ADF amplitude (output from ADF_recons, best fit values) Xants : (N,3) positions of antennas: x North, y West, z above sea level (m) Xsource : (3,) position of Xsource (from SWF) Bvec : (3,) magnetic field groundAltitude : height above sea level of the frame’s origin (m): 1231 m (GP13) by default; for simulation files, the ground altitude grand.analysis.geom.antenna_positions_from_run returns

Returns:

• eta ((N,) azimuthal angle in shower plane)

• omega ((N,) angle wrt shower axis)

• omega_cr ((N,) Cherenkov angle for each antenna)

• l_ant ((N,) distance from Xsource to antenna)

• adf ((N,) ADF amplitude for each antenna)

grand.analysis.fitting.adf.ADF_loss(params, Aants, Xants, Xsource, uncertainty=0.075)

Compute chi² for the ADF function.

Parameters:

• params (sequence of float) – [theta, phi, delta_omega, amplitude].

• Aants (ndarray) – Measured peak amplitudes, shape (N,).

• Xants (ndarray) – Antenna positions, shape (N, 3).

• Xsource (ndarray) – Shower source position, shape (3,) (from SWF).

• uncertainty (float, optional) – Relative uncertainty on amplitudes (default: 7.5%).

Returns:

chi² value.

Return type:

float

grand.analysis.fitting.adf.recons_ADF(theta_pwf, phi_pwf, Aants, Xants, Xsource)

Fit the ADF parameters to measured antenna peak amplitudes.

This function performs a chi² minimization using the ADF_loss function to reconstruct the best-fit shower parameters: direction (theta, phi), ADF width (delta_omega), and amplitude.

Inputs:

theta_pwf : initial guess for shower zenith angle (rad) phi_pwf : initial guess for shower azimuth angle (rad) Aants : measured peak amplitudes at antennas (N,) Xants : antenna positions (N,3) Xsource : Xsource position (3,)

Returns:

• theta_adf (reconstructed zenith angle (rad))

• phi_adf (reconstructed azimuth angle (rad))

• delta_omega (best-fit ADF shape parameter)

• amplitude (best-fit ADF amplitude)

grand.analysis.fitting.adf.ADF_fun(l_ant, amplitude, omega_cr, delta_omega)

Compute a simple model of the ADF for a shower on all omega values.

Don’t consider geomagnetic effect (very low) and average cherenkov angle along all eta values.

Parameters:

• l_ant (float) – Mean distance from the antennas to the shower source (in meters).

• amplitude (float) – Scaling factor of the ADF.

• omega_cr (float) – Mean Cherenkov opening angle (in radians).

• delta_omega (float) – Width parameter of the ADF.

Returns:

• omega (np.ndarray) – Array of angles (in radians) over which the ADF is evaluated.

• f_adf (np.ndarray) – Corresponding ADF amplitudes at each angle.

grand.analysis.energy_reco.voltage

Electromagnetic-energy proxy from the fitted radio amplitude.

grand.analysis.energy_reco.voltage.recons_energy_from_voltage(amplitude, sin_alpha, a=19600000.0, b=7900000.0)

Reconstruct the electromagnetic energy from the fitted radio amplitude.

The amplitude is in ADC counts. Just a first proxy.

Parameters:

• amplitude (float or array-like) – Scaling factor from the best-fit ADF (in ADC units).

• sin_alpha (float or array-like) – Sine of the geomagnetic angle between the shower axis and the geomagnetic field.

• a (float) – Slope parameter obtained from the fit on simulations (default: 1.96e7).

• b (float) – Offset parameter obtained from the fit on simulations (default: 7.90e6).

• arXiv (See)

Returns:

energy – Reconstructed electromagnetic energy in eV. Events yielding negative reconstructed energies are set to zero. Such cases may correspond to genuine cosmic-ray events for which the voltage-based energy proxy fails, as this reconstruction is not very robust and should only be interpreted as a first-order estimator.

Return type:

float or array-like

grand.analysis.signals.extraction

Peak amplitude and peak time extraction from antenna traces.

grand.analysis.signals.extraction.get_peak_amplitude(trace, channels, return_envelope=False)

Compute the peak amplitude of a trace (in ADC counts or uV or uV/m).

The envelope is the Euclidean norm of the Hilbert envelopes of the selected components, sqrt(sum_i |hilbert(ch_i)|**2). (The Hilbert envelope of the norm would be biased by -4% to +6%, depending on the carrier frequency.)

Parameters:

• trace (np.ndarray) – 2D array of shape (n_channels, n_samples).

• channels (list or array-like) – Indices of the channels to use for the amplitude computation.

• return_envelope (bool, optional) – If True, also return the Hilbert envelope of the signal.

Returns:

• peak_amp (float) – Maximum amplitude of the Hilbert envelope.

• hilbert_amp (np.ndarray, optional) – Hilbert envelope of the signal (only if return_envelope=True).

grand.analysis.signals.extraction.compute_t0(t_object)

Compute the initial time (t0) for each trace in nanoseconds.

Parameters:

t_object (object) – Object containing ‘du_seconds’ and ‘du_nanoseconds’ attributes.

Returns:

Array of t0 values for each trace in nanoseconds.

Return type:

np.ndarray

grand.analysis.signals.extraction.get_peak_time(trace, t0, channels, dt_ns=2)

Compute the time of the peak signal amplitude.

The peak time is obtained from the maximum of the Hilbert envelope and converted from sample index to physical time.

Parameters:

• trace (np.ndarray) – 2D array of shape (n_channels, n_samples).

• t0 (float or np.ndarray) – Initial time offset(s) in nanoseconds.

• channels (list or array-like) – Indices of the channels to use.

• dt_ns (float, optional) – Sampling time in nanoseconds (default is 2 ns).

Returns:

Time of the peak amplitude in seconds.

Return type:

float or np.ndarray

grand.analysis.signals.extraction.convert_voltage_to_ADC(trace, channels, adc_full_scale=8192, voltage_ref=0.9)

Convert voltage traces to ADC counts.

As grand.ADC.process() does, without its noise: the counts are truncated toward zero and saturate at +-adc_full_scale.

Parameters:

• trace (np.ndarray) – 2D array of voltage traces in microvolts (µV).

• channels (int, list, slice or boolean mask) – The channels to convert, selected as NumPy indexing does.

• adc_full_scale (int, optional) – ADC full scale value (default: 8192).

• voltage_ref (float, optional) – Reference voltage in volts (default: 0.9 V).

Returns:

Trace converted to ADC counts, as floats holding whole numbers; the channels not selected are returned unchanged.

Return type:

np.ndarray

Raises:

ValueError – If a selected channel holds NaN or infinity.

grand.analysis.cramer_rao_bounds.cramer_rao

Cramér-Rao lower bounds on the parameters of the reconstruction fits.

Each bound is the square root of a diagonal element of the inverse Fisher information matrix, built from numerical derivatives of the fit models at the fitted parameters and the assumed measurement uncertainties. It is the smallest standard deviation an unbiased fit can reach, so it gives the order of magnitude of the reconstruction errors.

grand.analysis.cramer_rao_bounds.cramer_rao.CRB_ADF_SWF(theta_swf: float, phi_swf: float, r_xsource: float, t_s: float, theta_adf: float, phi_adf: float, delta_omega: float, scaling_factor: float, Xants: ndarray, uncertainty_amplitude=0.075, uncertainty_time=5e-09) → ndarray

Cramer-Rao Bound (CRB) calculation for the joint ADF (amplitudes) and SWF (times) model.

Parameters:

• theta_swf (float) – Zenith angle from SWF reconstruction (radians).

• phi_swf (float) – Azimuth angle from SWF reconstruction (radians).

• r_xsource (float) – Distance to the source (meters).

• t_s (float) – Emission time of the source (seconds).

• theta_adf (float) – Zenith angle from ADF reconstruction (radians).

• phi_adf (float) – Azimuth angle from ADF reconstruction (radians).

• delta_omega (float) – Width parameter of the ADF, in radians (not an uncertainty).

• scaling_factor (float) – Scaling factor for the amplitude model.

• Xants (np.ndarray) – Antenna positions, shape (N, 3).

• uncertainty_amplitude (float, optional) – Relative uncertainty in amplitude measurements (current 7.5%).

• uncertainty_time (float, optional) – Uncertainty in time measurements (seconds).

Returns:

The smallest standard deviations (square roots of the diagonal of the inverse Fisher matrix) an unbiased fit can reach, for the parameters in the order of the arguments: theta_swf, phi_swf (radians), r_xsource (m), t_s (s), theta_adf, phi_adf, delta_omega (radians), scaling_factor. All NaN if the Fisher matrix cannot be inverted. They are bounds for the joint time-and-amplitude model, so lower than a time-only spherical fit can reach.

Return type:

np.ndarray, shape (8,)

grand.analysis.cramer_rao_bounds.cramer_rao.CRB_PWF(theta_pwf: float, phi_pwf: float, Xants: ndarray, uncertainty_time=5e-09) → ndarray

Cramer-Rao Bound (CRB) calculation for the Plane Wave Fit (PWF) model.

Parameters:

• theta_pwf (float) – Polar angle from PWF reconstruction (radians).

• phi_pwf (float) – Azimuth angle from PWF reconstruction (radians).

• Xants (np.ndarray) – Antenna positions, shape (N, 3).

• uncertainty_time (float, optional) – Uncertainty in time measurements (seconds).

Returns:

The smallest standard deviations of theta_pwf and phi_pwf an unbiased fit can reach, in radians (square roots of the diagonal of the inverse Fisher matrix). At zenith 0 or pi the azimuth is undefined: its bound is inf, with a warning, and the zenith bound is the one for a known azimuth.

Return type:

np.ndarray, shape (2,)

grand.analysis.physics.cherenkov_angle

Cherenkov angle from a two-emission-point model around the source.

grand.analysis.physics.cherenkov_angle.compute_Cerenkov(Xant, K, xsourceDist, Xsource, delta)

Compute Cherenkov angle by minimizing the time delay between light rays from shower points and the observer.

Inputs: - Xant : np.array, shape (3,) -> antenna position in meters - K : np.array, shape (3,) -> shower direction (unit vector) - xsourceDist : float -> distance from Xsource to shower core along K - Xsource : np.array, shape (3,) -> position of Xsource - delta : float -> distance along shower axis to points before/after Xmax

Returns:

- omega_cr

Return type:

float -> Cherenkov angle in radians

grand.analysis.physics.cherenkov_angle.compute_delay(omega, Xmax, Xa, Xb, Xant, U, K, alpha, delta, xmaxDist)

Compute residual time delay for a given Cherenkov angle guess.

Inputs: - omega : float -> guessed Cherenkov angle (radians) - Xmax : np.array, shape (3,) -> shower maximum position - Xa, Xb : np.array, shape (3,) -> points before/after Xmax along shower axis - Xant : np.array, shape (3,) -> antenna position - U : np.array, shape (3,) -> unit vector from shower core to antenna - K : np.array, shape (3,) -> shower direction (unit vector) - alpha : float -> angle between shower direction and vector to antenna - delta : float -> distance along shower axis from Xmax - xmaxDist : float -> distance from Xmax to shower core along K

Returns:

- res

Return type:

float -> residual (should be zero for correct Cherenkov angle)

grand.analysis.physics.cherenkov_angle.minor_equation(omega, n2, n1, alpha, delta, xmaxDist)

Compute time delay (in m).

Compute [c*delta(t)]^2

grand.analysis.physics.cherenkov_angle.compute_observer_position(omega, Xmax, Xant, U, K, xmaxDist, alpha)

Compute observer (antenna) position given Cherenkov angle.

Inputs: - omega : float -> Cherenkov angle - Xmax : np.array, shape (3,) -> shower maximum position - Xant : np.array, shape (3,) -> antenna position - U : np.array, shape (3,) -> unit vector from shower core to antenna - K : np.array, shape (3,) -> shower direction - xmaxDist : float -> distance from Xmax to shower core - alpha : float -> angle between shower direction and vector to antenna

Returns:

- X

Return type:

np.array, shape (3,) -> computed observer position

grand.analysis.physics.cherenkov_angle.rotation(angle, axis)

Compute 3x3 rotation matrix around axis using Rodrigues formula.

Inputs: - angle : float -> rotation angle in radians - axis : np.array, shape (3,) -> unit rotation axis

Returns:

- mat

Return type:

np.array, shape (3,3) -> rotation matrix

grand.analysis.physics.cherenkov_angle.der(func, x, args=[], eps=1e-07)

Forward estimate of derivative.

grand.analysis.physics.cherenkov_angle.newton(func, x0, tol=1e-07, nstep_max=100, args=[], verbose=False)

Newton-Raphson zero-finding method with numerical derivative.

Inputs: - func : callable -> function for which zero is sought - x0 : float -> initial guess - tol : float -> relative tolerance - nstep_max : int -> maximum iterations - args : list -> extra arguments to func - verbose : bool -> print iteration info

Returns:

- x

Return type:

float -> estimated zero of func

grand.analysis.physics.atmosphere

Refractive index of the atmosphere, locally and averaged along a path.

grand.analysis.physics.atmosphere.RefractionIndexAtPosition(X)

Return the refractive index at position X.

Exponential refractivity model, with the altitude taken over a spherical Earth.

Parameters:

X (array-like) – Position in meters, shape (3,).

Returns:

Refractive index.

Return type:

float

grand.analysis.physics.atmosphere.ZHSEffectiveRefractionIndex(X0, Xa)

Return the refractive index averaged between X0 and Xa.

As in ZHS: the refractivity is integrated along the straight path from the emission point to the antenna.

Parameters:

• X0 (array-like) – Emission point in meters, shape (3,).

• Xa (array-like) – Antenna position in meters, shape (3,).

Returns:

Effective refractive index along the path.

Return type:

float

grand.analysis.coords.array_shower

Shower direction vectors and the rotation into the shower frame.

grand.analysis.coords.array_shower.shower_direction_vector(theta, phi)

Returns the unit vector the shower propagates along.

It points away from where the shower comes from: the opposite of grand.dataio.xmax_frame.arrival_direction(), which points toward the source.

Parameters:

• theta (float) – Zenith angle of the direction the shower comes from, in radians.

• phi (float) – Azimuth of that direction, in radians, from north toward west.

Returns:

K, in the frame x North, y West, z Up; K[2] <= 0.

Return type:

ndarray, shape (3,)

grand.analysis.coords.array_shower.transformation_matrix(theta, phi, Bvec)

Builds the rotation matrix to the shower frame.

theta : zenith angle (rad) phi : azimuth angle (rad) Bvec : global magnetic field vector (3,)

grand.analysis.coords.array_shower.to_shower_frame(theta, phi, Bvec, Xant, Xsource)

Transforms global coordinates -> shower frame.

Xant : (3,) or (N,3) antenna positions Xsource : (3,) core/Xmax position :returns: coordinates in the shower frame :rtype: X_shower

grand.analysis.geom.angles

Angles between the shower axis, the antennas and the geomagnetic field.

grand.analysis.geom.angles.eta(theta, phi, Bvec, Xants, Xsource)

Computes the angle eta (azimuth angle in the shower plane).

theta : shower zenith angle (rad) phi : shower azimuth angle (rad) Bvec : magnetic field vector (3,) Xants : antenna positions (3,) or (N,3) Xsource: Xsource position in GRAND detector frame (3,)

grand.analysis.geom.angles.distance_source_antenna(Xants, Xsource)

Computes the distance(s) between source and antenna(s).

Xants : antenna positions (3,) or (N,3) Xsource: Xsource position in GRAND detector frame (3,)

grand.analysis.geom.angles.omega(theta, phi, Xants, Xsource)

Computes the angle omega between shower direction and antenna vector(s).

theta : shower zenith angle (rad) phi : shower azimuth angle (rad) Xants : antenna positions (3,) or (N,3) Xsource: Xsource position in GRAND detector frame (3,)

grand.analysis.geom.angles.sin_geomag_angle(theta, phi, B=None)

Compute the sine of the geomagnetic angle (alpha).

The angle is between the shower axis and the geomagnetic field.

Parameters:

• theta (float or array-like) – Zenith angle(s) of the shower in radians.

• phi (float or array-like) – Azimuth angle(s) of the shower in radians.

• B (array-like, shape (3,), optional) – Geomagnetic field vector (default: constants.Bn, read when called).

Returns:

sin_alpha – Sine of the geomagnetic angle

Return type:

float or ndarray

grand.analysis.geom.footprint

Shower core and Cherenkov footprint on the ground.

grand.analysis.geom.footprint.compute_core(k, Xsource, groundAltitude=1231)

Compute the shower core position.

The core is the intersection of the shower axis with the ground plane.

Warning

This function assumes a flat ground surface and does not take local topography into account. TODO: extend the implementation to include realistic ground elevation models.

Parameters:

• k (ndarray, shape (3,)) – Shower direction unit vector.

• Xsource (ndarray, shape (3,)) – Emission point coordinates [m].

• groundAltitude (float) – Height above sea level of the frame’s origin, where the source distance is measured from (meters). The default, 1231 m, is the GP13 site; for simulation files pass the ground altitude grand.analysis.geom.antenna_positions_from_run() returns.

Returns:

xc – Core position at ground [m].

Return type:

ndarray, shape (3,)

grand.analysis.geom.footprint.antenna_positions_from_run(trun)

Returns a run’s antenna positions in the fits’ frame, and its ground altitude.

The fits take x North, y West and z the height above sea level, with the source measured from (0, 0, groundAltitude). TRun.du_xyz gives z relative to the run’s origin, origin_geoid, so feeding it to the fits with the default groundAltitude (1231 m, GP13) placed the source that far too high: arrival times moved by up to 250 ns.

Parameters:

trun (grand.dataio.TRun) – The run tree, at the entry of the wanted run.

Returns:

• Xants (ndarray, shape (N, 3)) – Positions in meters, z above sea level, in trun.du_id order.

• ground_altitude (float) – origin_geoid’s height, to pass as groundAltitude.

grand.analysis.geom.footprint.generate_cone_surface_vectors(k, omega, n=10)

Generate vectors uniformly distributed on a cone surface.

This cone represents the Cherenkov emission cone. The cone axis is aligned with the shower direction, and the opening angle corresponds to the Cherenkov angle. The generated vectors sample the surface of this cone and are used to compute the projected Cherenkov footprint on the ground.

Parameters:

• k (ndarray, shape (3,)) – Unit vector defining the cone axis (shower propagation direction).

• omega (float) – Cone opening angle in radians (Cherenkov angle).

• n (int, optional) – Number of vectors generated uniformly around the cone (default: 10).

Returns:

• vectors (ndarray, shape (n, 3)) – Unit vectors lying on the surface of the Cherenkov cone, aligned

• with the shower axis.

grand.analysis.constants

Physical constants and site defaults used by the reconstruction.

The geomagnetic field at the GP13 DAQ position (B, and Bvec and Bn with x along magnetic north) is evaluated on first use, not at import. Evaluating it needs the compiled geometry core and the geomagnetic-model data, so doing it at import made import grand.analysis need both – and broke the documentation build, which stands a mock in for the core.

grand.analysis.constants.groundAltitude = 1231

Height above sea level, in meters, of the reconstruction frame’s origin: the GP13 site. The fits take antenna heights above sea level; the simulation files store du_xyz relative to origin_geoid instead (1264 m for the committed Xiaodushan samples), so use grand.analysis.geom.antenna_positions_from_run() for those (#252).