Signals and traces

Traces, signal processing and detector layouts: grand.basis.

grand.basis.traces_event	Handling a set of 3D traces
grand.basis.signal
grand.basis.type_trace	The two kinds of trace a detector unit carries.
grand.basis.du_network	Handling Detector Unit (DU) network, footprint plots
grand.basis.pipeline	Functions used in the pipline execution.

grand.basis.traces_event

Handling a set of 3D traces

grand.basis.traces_event.get_psd(trace, f_samp_mhz, nperseg=0)

Reference estimation of power spectrum density by Welch’s method

Parameters:

• trace (ndarray) – Time trace to transform.

• f_samp_mhz (float) – Sampling frequency, in MHz.

• nperseg (int, optional) – Segment length for the periodogram.

Returns:

Frequency axis in MHz, and power spectral density.

Return type:

tuple of ndarray

class grand.basis.traces_event.Handling3dTraces(name='NotDefined')

Bases: object

Handling a set of traces associated to one event observed on Detector Unit network

Initialisation, with two methods:

• init_traces()

• optionally with init_network()

Features:

• some plots : trace , power spectrum, footprint, …

• compute time where the trace is maximun

• compute time for each sample of trace

• compute a common time for trace

Public attributes:

• name str: name of the set of trace

• traces float(nb_du, 3, nb_sample): trace 3D

• idx2idt int(nb_du): array of identifier of DU

• t_start_ns float(nb_du): [ns] time of first sample of trace

• t_samples float(nb_du, nb_dim, nb_sample): [ns]

• f_samp_mhz float(nb_du,): [MHz] frequency sampling

• idt2idx dict: for each identifier return the index in array

• unit_trace str: string unit of trace

• network object: content position network

init_traces(traces, du_id=None, t_start_ns=None, f_samp_mhz=2000)

Parameters:

• traces (ndarray, shape (n_du, 3, n_samples)) – The traces.

• du_id (sequence, optional) – Identifier of each unit.

• t_start_ns (ndarray, optional) – Start time of each trace, in nanoseconds.

• f_samp_mhz (float or ndarray, optional) – Sampling frequency, in MHz.

init_network(du_pos)

Parameters:

du_pos (ndarray, shape (n_du, 3)) – Position of each unit, in the array frame.

set_unit_axis(str_unit='TBD', axis_name='idx', type_tr='Trace')

Parameters:

• str_unit (str, optional) – Unit of the trace values, for axis labels.

• axis_name (str, optional) – Naming of the three components.

• type_tr (str, optional) – Kind of trace, used in plot titles.

set_periodogram(size)

Parameters:

size (int) – Segment length used when computing spectra.

apply_bandpass(fr_min, fr_max, causal=True)

band filter with butterfly window

Parameters:

• fr_min (float) – Lower edge, in MHz.

• fr_max (float) – Upper edge, in MHz.

• causal (bool, optional) – Apply a causal filter rather than a zero-phase one.

keep_only_trace_with_ident(l_idt)

Keep trace with identifier defined in list <l_idt>

Parameters:

l_idt (sequence) – Identifiers of the units to keep.

keep_only_trace_with_index(l_idx)

Keep trace at index defined in list <l_idx>

Parameters:

l_idx (sequence of int) – Indices of the traces to keep.

reduce_nb_trace(new_nb_du)

reduces the number of traces to the first <new_nb_du>

Feature to reduce computation, for debugging

Parameters:

new_nb_du (int) – Number of traces to keep, taken from the start.

downsize_sampling(fact)

Downsampling with scipy decimate function

Parameters:

fact (int) – Decimation factor.

remove_trace_low_signal(threshold, norm_traces=None)

Remove trace where <norm_traces> is lower than <threshold>

Parameters:

• threshold (float) – Minimum amplitude to keep.

• norm_traces (ndarray, optional) – Precomputed norms, to avoid recomputing them.

Returns:

Indices of the traces that were kept.

Return type:

ndarray

copy(new_traces=None, deepcopy=True)

Return a copy of current object where traces can be modify

The type of copy is copy with reference, not a deepcopy https://stackoverflow.com/questions/3975376/why-updating-shallow-copy-dictionary-doesnt-update-original-dictionary/3975388#3975388

if new_traces is :

• None : object with same value

• 0 : the return object has a traces with same shape but set to 0

• np.array : the return object has new_traces as traces

Parameters:

• new_traces (ndarray, optional) – Traces for the copy; the existing ones by default.

• deepcopy (bool, optional) – Copy the arrays rather than referencing them.

Returns:

The copy.

Return type:

Handling3dTraces

get_delta_t_ns()

Return sampling rate in ns

Returns:

Sampling interval of each trace, in nanoseconds.

Return type:

ndarray

get_max_abs()

Find absolute maximal value in trace for each detector

Returns:

Largest absolute value of each trace.

Return type:

ndarray

get_max_norm()

Return array of maximal of 3D norm in trace for each detector

Returns:

Largest three-component norm of each trace.

Return type:

ndarray

get_norm()

Return norm of traces for each time sample

Returns:

Three-component norm at every sample.

Return type:

ndarray, shape (n_du, n_samples)

get_tmax_vmax(hilbert=True, interpol='parab')

Return time where norm of the amplitude of the Hilbert tranform is max

Parameters:

• hilbert (bool, optional) – Use the Hilbert envelope rather than the raw trace.

• interpol (str, optional) – Interpolation around the peak: "parab" for a three-point parabola, or "auto".

Returns:

Time of maximum in nanoseconds, and the value there.

Return type:

tuple of ndarray

get_min_max_t_start()

Return time interval of time start of trace

Returns:

Earliest and latest trace start time, in nanoseconds.

Return type:

tuple of float

get_nb_trace()

Return the number of 3d traces

Returns:

Number of traces held.

Return type:

int

get_size_trace()

Return the number of sample in trace

Returns:

Number of samples per trace.

Return type:

int

get_snr_and_noise()

Return a crude estimation of SNR and noise level

Crude estimation because:

• noise is estimated at the end of trace

• max value has different estimator.

Returns:

• snr (ndarray) – Signal-to-noise ratio of each trace.

• v_max (ndarray) – Peak of each trace’s Hilbert envelope, in the traces’ unit.

• noise (ndarray) – Noise level of each trace (largest standard deviation over the three axes of the trace’s last samples), in the same unit.

get_extended_traces()

Return extended traces to time interval of event

Compute and return traces extended to the entire duration of the event with common time

Returns:

Traces padded onto a common time axis, and that axis.

Return type:

tuple of ndarray

plot_trace_idx(idx, to_draw='012')

Draw 3 traces associated to DU with index idx

Parameters:

• idx (int) – Index of the trace to plot.

• to_draw (str, optional) – Which components to draw.

plot_trace_du(du_id, to_draw='012')

Draw 3 traces associated to DU idx2idt

Parameters:

• du_id (int) – Identifier of the unit to plot.

• to_draw (str, optional) – Which components to draw.

plot_psd_trace_idx(idx, to_draw='012')

Draw power spectrum for 3 traces associated to DU at index idx

Parameters:

• idx (int) – Index of the trace.

• to_draw (str, optional) – Which components to draw.

plot_psd_trace_du(du_id, to_draw='012')

Draw power spectrum for 3 traces associated to DU idx2idt

Parameters:

• du_id (int) – Identifier of the unit.

• to_draw (str, optional) – Which components to draw.

plot_all_traces_as_image()

Interactive image double click open traces associated

plot_histo_t_start()

Histogram of time start

plot_footprint_4d_max()

Plot time max and max value by component

plot_footprint_val_max()

Plot footprint max value

plot_footprint_time_max()

Plot footprint time associated to max value

plot_footprint_time_slider()

Plot footprint max value

grand.basis.signal

• This module contains several signal processing functionalities to be applied to simulation/data

• operations are meant to be on the signal traces for individual antennas, suitable to be used both in Grandlib format/ read from hdf5 files

• expects signal traces to be of the size (3,lengthoftrace)

grand.basis.signal.find_max_with_parabola_interp_3pt(x_trace, y_trace, idx_max)

Returns the interpolated maximum of a trace, from a parabola through 3 points.

Fits \(ax^2 + bx + c\) through the largest sample and its two neighbors, offset to \((x_0, y_0)\), and returns the vertex:

\[x_m = x_0 - \frac{b}{2a}, \qquad y_m = y_0 - \frac{b^2}{4a}\]

The trace is assumed non-negative, and idx_max is assumed to be a genuine local maximum with a sample on either side.

Parameters:

• x_trace (ndarray) – Sample positions.

• y_trace (ndarray) – Sample values.

• idx_max (int) – Index of the largest sample.

Returns:

Position and value of the interpolated maximum, from a parabola through three points.

Return type:

tuple of float

Examples

A peak between two samples is recovered more precisely than the sample grid allows.

import numpy as np
from grand.basis.signal import find_max_with_parabola_interp_3pt

t = np.arange(512) * 0.5
true_peak = 100.25                              # deliberately between samples
y = np.exp(-((t - true_peak) ** 2) / (2 * 6.0 ** 2))

idx = int(np.argmax(y))
t_interp, v = find_max_with_parabola_interp_3pt(t, y, idx)
print("nearest sample : %.4f ns" % t[idx])
print("interpolated   : %.4f ns" % t_interp)

nearest sample : 100.0000 ns
interpolated   : 100.2500 ns

grand.basis.signal.find_max_with_parabola_interp(x_trace, y_trace, idx_max, factor_hill=0.96)

Parabolic interpolation of the maximum with more than 3 points

trace : all values >= 0

algo:

1. find begin idx, ie trace[–idx_max] > v_max*factor_hill

2. find end idx, ie trace[idx_max++] > v_max*factor_hill

3. if nb idx <= 2 : mode pic else mode hill

4. Mode pic : 3 values and the middle one is max 4.1 offset of (x0, v0) 4.2 solve coef a, b => x_m = offset - b/2a ; v_m=offset - b^2/4a

5. Mode hill: 5.0 offset of (x, y) of first sample 5.1 solve overdetermined linear system with a, b, c 5.2 x_m =offset - b/2a ; v_m=offset - b^2/4a + c

Parameters:

• x_trace (ndarray) – Sample positions.

• y_trace (ndarray) – Sample values.

• idx_max (int) – Index of the largest sample.

• factor_hill (float, optional) – Fraction of the peak defining how much of the hill to fit.

Returns:

Position and value of the interpolated maximum.

Return type:

tuple of float

grand.basis.signal.get_filter(time, trace, fr_min, fr_max)

Filter signal in given bandwidth

@note

At present Butterworth filter only is implemented, others: what is close to hardware filter?

Parameters:

• time (ndarray) – Time axis, in nanoseconds.

• trace (ndarray) – Trace to filter.

• fr_min (float) – Lower band edge, in Hz. Note the unit: most frequency arguments in this package are in MHz, and this one is not – passing 50 rather than 50e6 silently returns an all-zero trace.

• fr_max (float) – Upper band edge, in Hz.

Returns:

The band-passed trace.

Return type:

ndarray

Examples

The band edges are in hertz. Passing megahertz does not raise; it silently returns an empty trace, which is the failure most easily mistaken for a quiet event.

import numpy as np
from grand.basis.signal import get_filter

t = np.arange(1024) * 0.5                       # ns
tone = np.sin(2 * np.pi * 100.0 * t * 1e-3)     # a 100 MHz tone

kept = get_filter(t, tone, 50e6, 200e6)         # edges in Hz
print("in band, kept: %.3f" % np.std(kept))

# The same numbers read as MHz pass nothing at all -- no error, just zeros.
print("edges given in MHz: %.3f" % np.std(get_filter(t, tone, 50.0, 200.0)))

in band, kept: 0.701
edges given in MHz: 0.000

grand.basis.signal.get_peakamptime_norm_hilbert(a2_time, a3_trace)

Get peak Hilbert amplitude norm of trace (v_max) and its time t_max without interpolation

Parameters:

• a2_time (ndarray) – Time axis per trace.

• a3_trace (ndarray, shape (n_du, 3, n_samples)) – Traces.

Returns:

Peak time, peak amplitude, the norm, and its Hilbert envelope.

Return type:

tuple of ndarray

grand.basis.signal.get_fastest_size_fft(sig_size, f_samp_mhz, padding_fact=1)

Returns an FFT-friendly transform length and its frequency axis.

The same as grand.sim.efield2voltage.get_fastest_size_fft(), whose padding argument is named padding_factor.

Parameters:

• sig_size (int) – Length of the trace, in samples.

• f_samp_mhz (float or ndarray) – Sampling frequency, in MHz.

• padding_fact (float, optional) – Zero-padding factor; at least 1.

Returns:

Transform length, and the frequency axis in MHz.

Return type:

tuple

grand.basis.signal.interpol_at_new_x(a_x, a_y, new_x)

#RK: This function is copied to galaxy.py and rf_chain.py where it is used. Remove it from here if it is not used anywhere else.

Interpolation of discreet function F defined by set of point F(a_x)=a_y for new_x value and set to zero outside interval definition a_x

Parameters:

• a_x (ndarray) – Sample positions.

• a_y (ndarray) – Sample values.

• new_x (ndarray) – Positions to interpolate onto.

Returns:

Interpolated values, zero outside the range of a_x.

Return type:

ndarray

Examples

Values outside the measured range come back as zero rather than extrapolated: the tabulated antenna and RF-chain data have no meaning beyond 30-250 MHz.

import numpy as np
from grand.basis.signal import interpol_at_new_x

a_x = np.linspace(30.0, 250.0, 100)             # the measured band, in MHz
a_y = np.ones_like(a_x)

print(interpol_at_new_x(a_x, a_y, np.array([10.0, 140.0, 400.0])))

[0. 1. 0.]

grand.basis.type_trace

The two kinds of trace a detector unit carries.

ElectricField is the field at the antenna, before the instrument touches it; Voltage is what comes out of it. Both are dataclasses over a time axis, and both cache their own transform.

They are not symmetric, which is worth knowing before reaching for one: ElectricField carries the antenna position and the coordinate frame its three components are expressed in, and Voltage carries neither – only t and V. A voltage trace therefore cannot say what frame it is in, and the caller has to keep track.

class grand.basis.type_trace.ElectricField(a_time: ndarray, e_xyz: CartesianRepresentation, pos_xyz: CartesianRepresentation | None = None, frame: LTP | GRANDCS | None = None)

Bases: object

A three-component electric field trace, with its time axis and frame.

Holds the samples for one detection unit: three components against a shared time axis, together with the coordinate frame they are expressed in.

Parameters:

• a_time (np.ndarray) – Time of each sample, in s (the voltage chain passes ns times 1e-9)

• e_xyz (CartesianRepresentation) – The three components of the field, shape (3, n_samples), in µV/m

• pos_xyz (Union[CartesianRepresentation, None]) – Position of the detection unit, in m, in frame

• frame (Union[LTP, GRANDCS, None]) – Frame in which e_xyz and pos_xyz are expressed

a_time: ndarray

Time of each sample, in s (the voltage chain passes ns times 1e-9)

e_xyz: CartesianRepresentation

The three components of the field, shape (3, n_samples), in µV/m

pos_xyz: CartesianRepresentation | None = None

Position of the detection unit, in m, in frame

frame: LTP | GRANDCS | None = None

Frame in which e_xyz and pos_xyz are expressed

get_fft(size_sig_pad)

Parameters:

size_sig_pad (int) – Padded transform length.

Returns:

The one-sided spectrum of the three components.

Return type:

ndarray, shape (3, n_freq)

get_delta_time_s()

Returns the sampling interval, in seconds.

Returns:

Time between consecutive samples. Note the unit: the traces themselves are conventionally described in nanoseconds.

Return type:

float

get_nb_sample()

Returns the number of samples in the trace.

Returns:

Length of the time axis.

Return type:

int

class grand.basis.type_trace.Voltage(t: ndarray, V: ndarray)

Bases: object

A three-component voltage trace and its time axis.

Holds the samples for one detection unit: V against t. Unlike ElectricField it carries no position and no coordinate frame – the fields are t and V and nothing else – so what frame the components are in is the caller’s to remember.

Parameters:

• t (np.ndarray) – Time of each sample

• V (np.ndarray) – The three components, shape (3, n_samples), in the unit of the caller’s input (µV in the voltage chain)

t: ndarray

Time of each sample

V: ndarray

The three components, shape (3, n_samples), in the unit of the caller’s input (µV in the voltage chain)

get_fft(size_sig_pad)

Parameters:

size_sig_pad (int) – Padded transform length.

Returns:

The one-sided spectrum of the three components.

Return type:

ndarray, shape (3, n_freq)

get_delta_time_s()

Returns the sampling interval, in seconds.

Returns:

Time between consecutive samples. Note the unit: the traces themselves are conventionally described in nanoseconds.

Return type:

float

get_nb_sample()

Returns the number of samples in the trace.

Returns:

Length of the time axis.

Return type:

int

grand.basis.du_network

Handling Detector Unit (DU) network, footprint plots

grand.basis.du_network.closest_node(node, nodes)

Simple computation of distance mouse DU

Parameters:

• node (array_like) – Reference position.

• nodes (ndarray) – Positions to search.

Returns:

Index of the nearest of nodes.

Return type:

int

class grand.basis.du_network.DetectorUnitNetwork(name='NotDefined')

Bases: object

Handling a DU network

Public attributs:

• name str: name of the set of trace

• du_pos float(nb_du, 3): position of DU

• idx2idt int(nb_du): array of identifier of DU

init_pos_id(du_pos, du_id=None)

Init object with array DU position and identifier

Parameters:

• du_pos (ndarray, shape (n_du, 3)) – Position of each unit, in meters.

• du_id (sequence, optional) – Identifier of each unit.

keep_only_du_with_index(l_idx)

Keep DU at index defined in list <l_idx>

Parameters:

l_idx (sequence of int) – Indices to keep.

reduce_nb_du(new_nb_du)

Feature to debug and reduce computation

Parameters:

new_nb_du (int) – Number of units to keep, taken from the start.

get_sub_network(l_id)

Reduce networh to DU in list <l_id>

Parameters:

l_id (sequence) – Identifiers to keep.

Returns:

A network holding only those units.

Return type:

DetectorUnitNetwork

get_nb_du()

Return the number of DU

Returns:

Number of detection units.

Return type:

int

get_surface()

Return suface in km2

Returns:

Area covered by the array, in square kilometres.

Return type:

float

get_max_dist_du()

Return the largest distance between two DUs of the network.

Returns:

The distance, in kilometres; 0 for fewer than two DUs.

Return type:

float

plot_du_pos()

Plot DU position

plot_footprint_1d(a_values, title='', traces=None, scale='log', unit='')

Interactive footprint double click on DU draw trace associated and power spectrum

Parameters:

• a_values (ndarray) – Per-unit value to color by.

• title (str, optional) – Plot title.

• traces (Handling3dTraces, optional) – Traces to show when a unit is clicked.

• scale (str, optional) – "log" or "lin".

• unit (str, optional) – Unit of a_values, for the color bar.

plot_footprint_4d(o_tr, v_plot, title='', same_scale=True, unit='')

Plot footprint of time max by DU and value max by component

Parameters:

• o_tr (Handling3dTraces) – Traces to animate.

• v_plot (str, optional) – Quantity to plot.

• title (str, optional) – Plot title.

• same_scale (bool, optional) – Use one color scale for all components.

• unit (str, optional) – Unit for the color bar.

Returns:

Draws the figure.

Return type:

None

plot_footprint_time(a_time, a3_values, title='')

Interactive plot, footprint max value for time defined by user with slider widget

Parameters:

• a_time (ndarray) – Time axis, in nanoseconds.

• a3_values (ndarray) – Per-unit, per-time values.

• title (str, optional) – Plot title.

Returns:

Draws the figure.

Return type:

None

grand.basis.pipeline

Functions used in the pipline execution.

RK TODO: implement techniques used in LP’s execute_pipeline where output_trees

and time-traces are saved on event-by-event basis.

@GRAND Collaboration, 2023

grand.basis.pipeline.logger = <Logger grand.basis.pipeline (WARNING)>

Final usage: from grand.io.pipeline import Pipeline

pipeline = Pipeline() pipeline.Add(“reader”,f_input=”f.root”) pipeline.Add(“coreas2root”, …) pipeline.Add(“efield2voltage”, …..) pipeline.Add(“some_reco”, ….) pipeline.Add(“something_else”, …) pipeline.Add(“writer”, f_output=”out.root”)

class grand.basis.pipeline.Pipeline

Bases: object

A sketch of a processing pipeline; not finished.

The intent is that stages are added by name and events pushed through them. It is incomplete and has drifted from grand.dataio; see the module docstring.

trees_dict = {}

command_dict = {}

Add(name, **kwargs)

Adds a named stage to the pipeline.

Parameters:

• name (str) – Stage name, such as "reader" or "efield2voltage".

• **kwargs – Stage-specific options.