Coverage for src/pyhiperta/R0_DL1.py: 11%

1# Copyright 2026 CNRS 

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

3 

4"""HiPeRTA file processing pipeline in python.""" 

5 

6import argparse 

7 

8import numpy as np 

9import tables 

10 

11from pyhiperta.calibration import calibrate, select_channel 

12from pyhiperta.integration import integrate_local_peak, windowed_integration_correction 

13from pyhiperta.R0Reader import R0HDF5Dataset 

14from pyhiperta.waveform_indexing import waveform1Dto2D, waveform2Dto1D, waveform_1D_to_2D_maps 

15 

16 

17def main(): 

18 

19 import matplotlib.pyplot as plt 

20 from matplotlib.backends.backend_pdf import PdfPages 

21 

22 def plotImage(pdf, tabX, tabY, signal, imageIndex, eventId, eventType, tabFixPixelOrder=None): 

23 strInfo = f"image {imageIndex}, id = {eventId}, type = {eventType},\n Signal mean = {signal.mean()} pe, min = {signal.min()}, max = {signal.max()}, std = {signal.std()}" 

24 print(strInfo) 

25 fig = plt.figure(figsize=(16, 10)) 

26 fig.patch.set_alpha(1.0) 

27 if tabFixPixelOrder is not None: 

28 plt.scatter(tabX[tabFixPixelOrder], tabY[tabFixPixelOrder], c=signal, s=120) 

29 else: 

30 plt.scatter(tabX, tabY, c=signal, s=120) 

31 

32 plt.axis("equal") 

33 plt.xlabel("x") 

34 plt.ylabel("y") 

35 plt.colorbar() 

36 plt.text(0.05, 0.95, strInfo, transform=fig.transFigure, size=14) 

37 pdf.savefig() # saves the current figure into a pdf page 

38 plt.close() 

39 

40 parser = argparse.ArgumentParser() 

41 parser.add_argument("--R0_files", "-r0", type=str, required=True, help="input file, or folder.") 

42 parser.add_argument("--output_file", "-o", type=str, required=True, help="output path of produced pdf.") 

43 parser.add_argument("--first_image", "-f", type=int, default=1, help="Index of first image to plot.") 

44 parser.add_argument("--last_image", "-l", type=int, default=100, help="Index of last image to plot.") 

45 parser.add_argument( 

46 "--indexing_step", "-idx", type=int, default=1, help="Step between the images indices selected for plotting." 

47 ) 

48 

49 args = parser.parse_args() 

50 

51 dataset = R0HDF5Dataset(args.R0_files, 1) 

52 gains = dataset.read_gains() 

53 pedestals = dataset.read_pedestals() 

54 pedestals /= dataset.events_n_frames() # pedestal per frame option of the poor 

55 camera_geometry = dataset.read_camera_geometry() 

56 map_1D_to_2D = waveform_1D_to_2D_maps(dataset.read_camera_geometry()) 

57 

58 with tables.open_file(args.R0_files, "r") as hfile: 

59 tabEventId = hfile.root.r0.event.telescope.waveform.tel_001.col("event_id") 

60 tabEventType = hfile.root.r0.event.subarray.trigger.col("event_type") 

61 # tabPixelOrder = hfile.root.configuration.instrument.telescope.camera.pixel_order.tel_001.read() 

62 

63 ref_pulse_sample_times, ref_pulse_shape = dataset.read_reference_pulse() 

64 gain_correction = windowed_integration_correction( 

65 ref_pulse_sample_times, ref_pulse_shape, 1, lambda x: integrate_local_peak(x, 3, 3) 

66 ) 

67 print("gain correction: ", gain_correction) 

68 

69 with PdfPages(args.output_file) as pdf: 

70 for im_idx in range(args.first_image, args.last_image, args.indexing_step): 

71 R0_waveform_high, R0_waveform_low = dataset[im_idx] 

72 print( 

73 "Waveform high: ", 

74 R0_waveform_high.shape, 

75 R0_waveform_high.mean(), 

76 R0_waveform_high.min(), 

77 R0_waveform_high.max(), 

78 ) 

79 print( 

80 "Waveform low: ", 

81 R0_waveform_low.shape, 

82 R0_waveform_low.mean(), 

83 R0_waveform_low.min(), 

84 R0_waveform_low.max(), 

85 ) 

86 print("gains high: ", gains[0].shape, gains[0].mean(), gains[0].min(), gains[0].max()) 

87 print("gains low: ", gains[1].shape, gains[1].mean(), gains[1].min(), gains[1].max()) 

88 print("pedestals high: ", pedestals[0].shape, pedestals[0].mean(), pedestals[0].min(), pedestals[0].max()) 

89 print("pedestals low: ", pedestals[1].shape, pedestals[1].mean(), pedestals[1].min(), pedestals[1].max()) 

90 R0_waveform, selected_gains, selected_pedestals, selected_correction = select_channel( 

91 R0_waveform_high, R0_waveform_low, gains, pedestals, gain_correction, 4000 

92 ) 

93 R0_calibrated = calibrate(R0_waveform, selected_gains, selected_pedestals) 

94 R0_integrated = integrate_local_peak(R0_calibrated, 3, 3) 

95 R0_integrated *= selected_correction 

96 

97 R0_integrated = R0_integrated.squeeze() 

98 

99 plotImage( 

100 pdf, 

101 camera_geometry[0], 

102 camera_geometry[1], 

103 R0_integrated, 

104 im_idx, 

105 tabEventId[im_idx], 

106 tabEventType[im_idx], 

107 ) 

108 

109 fft_image = waveform2Dto1D(map_1D_to_2D, np.fft.fft2(waveform1Dto2D(map_1D_to_2D, R0_integrated))) 

110 print("fft_image: ", fft_image.shape, fft_image.mean(), fft_image.min(), fft_image.max()) 

111 plotImage( 

112 pdf, 

113 camera_geometry[0], 

114 camera_geometry[1], 

115 np.abs(fft_image), 

116 im_idx, 

117 tabEventId[im_idx], 

118 tabEventType[im_idx], 

119 ) 

120 

121 fig = plt.figure(figsize=(16, 10)) 

122 fig.patch.set_alpha(1.0) 

123 plt.imshow(waveform1Dto2D(map_1D_to_2D, R0_integrated)) 

124 plt.axis("equal") 

125 plt.xlabel("x") 

126 plt.ylabel("y") 

127 plt.colorbar() 

128 pdf.savefig() # saves the current figure into a pdf page 

129 plt.close() 

130 

131 fig = plt.figure(figsize=(16, 10)) 

132 fig.patch.set_alpha(1.0) 

133 plt.imshow(np.abs(np.fft.fft2(waveform1Dto2D(map_1D_to_2D, R0_integrated)))) 

134 plt.axis("equal") 

135 plt.xlabel("x") 

136 plt.ylabel("y") 

137 plt.colorbar() 

138 pdf.savefig() # saves the current figure into a pdf page 

139 plt.close() 

140 

141 

142if __name__ == "__main__": 142 ↛ 143line 142 didn't jump to line 143 because the condition on line 142 was never true

143 main()