Coverage for src/pyhiperta/integration.py: 100%

1# Copyright 2026 CNRS 

2# This software is distributed under the terms of the CeCILL-C free software license. 

3 

4"""Integration of waveforms videos into a single (charge, time of max) pair of image.""" 

5 

6from typing import Protocol 

7 

8import numpy as np 

9 

10 

11def integrate_local_peak( 

12 waveform: np.ndarray, nb_frames_before_max: int, nb_frames_after_max: int, time_bin_duration_ns: float = 1.0 

13) -> tuple[np.ndarray, np.ndarray]: 

14 """Integrate (sum the number of photoelectrons) in the waveform pulses in a window around the pulse maximum. 

15 

16 The method does the following: find the maximum in time for each pixel independently, then sum the frames 

17 from `max - nb_slice_before_peak` to `max + nb_slice_after_peak` (inclusive on both end). 

18 If the window wouldn't fit entirely in the available frames, because the maximum position is too 

19 close to one of the edges, the window is shifted towards the center until it includes the same number 

20 of frames. 

21 

22 Parameters 

23 ---------- 

24 waveform : np.ndarray 

25 Batch or single video samples of the shower events. Shape: (N_batch, N_frames, N_pixels) 

26 nb_frames_before_max : int 

27 Number of frames before the maximum to include in the integration window. 

28 nb_frames_after_max : int 

29 Number of frames after the maximum to include in the integration window. 

30 time_bin_duration_ns : float 

31 Amount of time between 2 frames in the waveforms. It is used to compute the maximum signal time. 

32 

33 Returns 

34 ------- 

35 Tuple[np.ndarray, np.ndarray] 

36 Signal charge (integrated waveform signal) and peak maximum time. 

37 

38 Examples 

39 -------- 

40 >>> waveforms = 10.0 * np.random.random((12, 40, 1855)) # "random shower samples" 

41 >>> images = integrate_local_peak(waveforms, 3, 3) # shape (12, 1855)1 

42 

43 Raises 

44 ------ 

45 ValueError 

46 If the number of frames to include before or after the maximum is strictly negative. 

47 """ 

48 if nb_frames_before_max < 0 or nb_frames_after_max < 0: 

49 raise ValueError( 

50 "Number of frames before or after maximum must be greater or equal to 0. " 

51 f"Got nb_frames_before_max {nb_frames_before_max} and nb_frames_after_max {nb_frames_after_max}" 

52 ) 

53 # Find maximum position (argmax). Window center position will be the maximum position if it allows the entire 

54 # window to be in the array, otherwise the window is shifted to be in the array (clip). 

55 peaks_idx = np.clip( 

56 np.argmax(waveform, axis=-2, keepdims=True), 

57 nb_frames_before_max, 

58 waveform.shape[1] - nb_frames_after_max - 1, 

59 dtype=np.int32, 

60 ) 

61 # Take the windows in the array (different slice for each pixel) 

62 # See example with argmax in `take_along_axis` documentation. 

63 # +1 in arange to include the last slice 

64 window_idx = ( 

65 peaks_idx + np.arange(-nb_frames_before_max, nb_frames_after_max + 1, dtype=np.int32)[..., np.newaxis] 

66 ) 

67 waveform_windows = np.ascontiguousarray(np.take_along_axis(waveform, window_idx, axis=-2)) 

68 # charge is just the sum 

69 charge = waveform_windows.sum(axis=-2) 

70 # max time is sum(waveform * t) / sum(waveform). 

71 max_times = np.float32(time_bin_duration_ns) * ( 

72 (waveform_windows * window_idx).sum(axis=-2, dtype=np.float32) / charge 

73 ) 

74 return charge, max_times 

75 

76 

77class IntegrationFunctionType(Protocol): 

78 """Describes the signature of a function that can be used for windowed integration and associated correction.""" 

79 

80 def __call__( 

81 self, waveform: np.ndarray, nb_frames_before_max: int, nb_frames_after_max: int, *args, **kwargs 

82 ) -> tuple[np.ndarray, np.ndarray]: ... 

83 

84 

85def windowed_integration_correction( 

86 ref_pulse_sample_times_ns: np.ndarray, 

87 ref_pulse_waveform: np.ndarray, 

88 integration_function: IntegrationFunctionType, 

89 nb_frames_before_max_real_data: int, 

90 nb_frames_after_max_real_data: int, 

91 real_data_sampling_period_ns: float, 

92) -> np.ndarray: 

93 """Calculate the correction to apply to compensate the bias introduced by the windowed integration of real data. 

94 

95 The windowed integration results are biased with respect to an integration on the total time range, because the 

96 integration window is typically shorter than the time range of the full signal (this is on purpose: the edges of 

97 the signals are cut-off because they are closer to noise level than the pulse center). 

98 

99 The bias introduced depends on the integration window parameters: position around maximum and range. To be able to 

100 compare integrated charge computed with different integration windows, we need to be able to correct the biases. 

101 The correction is such that 

102 correction * windowed_integration(noise_less_signal, n_slice_before, n_slice_after) = integration(noise_less_signal) 

103 

104 The noise less signal is called the "reference pulse" and typically provided by camera teams, it is the response of 

105 a pixel to a single photo-electron. 

106 

107 The waveform signals (reference pulse and real data) come as as number of p.e. at a certain time, ie it is binned 

108 number of p.e. wrt time. To compute the correction, the reference pulse is re-binned (typically it has a different 

109 bin size than real data) to real data bins which introduces a small error: the ratio between real data bins and the 

110 reference bins is not an integer, so the bins on the edge of the integration window have a signal that is a bit 

111 mixed with the neighbor's bin (outside of the window). 

112 

113 Parameters 

114 ---------- 

115 ref_pulse_sample_times_ns : np.ndarray 

116 Array containing the timestamps at which the number of photoelectrons are measured, for the reference pulse. 

117 Shape: (N_samples,) 

118 ref_pulse_waveform : np.ndarray 

119 Array containing the number of photoelectrons at each timestamp, for the reference pulse. ref_pulse_waveform[0] 

120 is the reference pulse of the channel 0 (low gain) and ref_pulse_waveform[1] is the reference pulse for 

121 channel 1 (high gain) 

122 Shape: (2, N_samples) 

123 integration_function : IntegrationFunctionType 

124 Function performing a window integration. It must accept the "waveform", "nb_frames_before_max" and 

125 "nb_frames_after_max" arguments. 

126 nb_frames_before_max_real_data : int 

127 Number of frames before the maximum to include in the integration window, when used with real data. 

128 The number of frames of the reference pulse to integrate will be computed based on the reference pulse sampling 

129 period and `real_data_sampling_period_ns`. 

130 nb_frames_after_max_real_data : int 

131 Number of frames after the maximum to include in the integration window, when used with real data. 

132 The number of frames of the reference pulse to integrate will be computed based on the reference pulse sampling 

133 period and `real_data_sampling_period_ns`. 

134 real_data_sampling_period_ns : float 

135 Sampling period (amount of time between 2 samples) of the real data, in nanoseconds. 

136 

137 Returns 

138 ------- 

139 np.ndarray 

140 Correction value for each channel such that 

141 correction * windowed_integration(noiseless_real_data) = full_integration(noiselesss_real_data) 

142 Shape: (2,) 

143 

144 Raises 

145 ------ 

146 TypeError 

147 If `integration_function` does not accept the expected arguments. 

148 """ 

149 ref_pulse_real_data_bins_times = np.arange( 

150 ref_pulse_sample_times_ns.min(), 

151 ref_pulse_sample_times_ns.max(), 

152 real_data_sampling_period_ns, 

153 dtype=np.float32, 

154 ) 

155 # for each gain: re-bin the reference pulse with bins the size of the real data bins 

156 # (preserving the total number of p.e. but getting bins as wide as the real data bins) 

157 # Note: keep only the 1st return value of histogram since we only want the counts, we know the bin edges 

158 real_data_bin_ref_pulse = np.stack( 

159 [ 

160 np.histogram(ref_pulse_sample_times_ns, bins=ref_pulse_real_data_bins_times, weights=pulse_waveform)[0] 

161 for pulse_waveform in ref_pulse_waveform 

162 ], 

163 axis=0, 

164 ) 

165 # Compute the windowed integration for each gain 

166 try: 

167 window_integrated_ref_pulse = integration_function( 

168 waveform=real_data_bin_ref_pulse[..., np.newaxis], # add dimension to "simulate" nb_pixel dimension 

169 nb_frames_before_max=nb_frames_before_max_real_data, 

170 nb_frames_after_max=nb_frames_after_max_real_data, 

171 )[0].squeeze() # and remove this dimension (after selecting integrated charge and discarding max times) 

172 except TypeError as e: 

173 raise TypeError( 

174 "Got a type error when evaluating the integration_function. " 

175 "Make sure that the integration function has the expected arguments names: " 

176 "waveform, nb_frames_before_max, nb_frames_after_max." 

177 ) from e 

178 # correction is the ratio of total integration and windowed integration 

179 return real_data_bin_ref_pulse.sum(axis=-1) / window_integrated_ref_pulse