Coverage for pyhiperta/R0_DL1.py: 0%

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

2 

3import numpy as np 

4import tables 

5 

6from pyhiperta.calibration import calibrate, select_channel 

7from pyhiperta.integration import integrate_local_peak, windowed_integration_correction 

8from pyhiperta.R0Reader import R0HDF5Dataset 

9from pyhiperta.waveform_indexing import waveform1Dto2D, waveform2Dto1D, waveform_1D_to_2D_maps 

10 

11 

12def main(): 

13 import argparse 

14 

15 import matplotlib.pyplot as plt 

16 from matplotlib.backends.backend_pdf import PdfPages 

17 

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

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

20 print(strInfo) 

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

22 fig.patch.set_alpha(1.0) 

23 if tabFixPixelOrder is not None: 

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

25 else: 

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

27 

28 plt.axis("equal") 

29 plt.xlabel("x") 

30 plt.ylabel("y") 

31 plt.colorbar() 

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

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

34 plt.close() 

35 

36 parser = argparse.ArgumentParser() 

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

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

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

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

41 parser.add_argument( 

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

43 ) 

44 

45 args = parser.parse_args() 

46 

47 dataset = R0HDF5Dataset(args.R0_files, 1) 

48 gains = dataset.read_gains() 

49 pedestals = dataset.read_pedestals() 

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

51 camera_geometry = dataset.read_camera_geometry() 

52 map_1D_to_2D = waveform_1D_to_2D_maps(dataset.read_camera_geometry()) 

53 

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

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

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

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

58 

59 ref_pulse_sample_times, ref_pulse_shape = dataset.read_reference_pulse() 

60 gain_correction = windowed_integration_correction( 

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

62 ) 

63 print("gain correction: ", gain_correction) 

64 

65 with PdfPages(args.output_file) as pdf: 

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

67 R0_waveform_high, R0_waveform_low = dataset[im_idx] 

68 print( 

69 "Waveform high: ", 

70 R0_waveform_high.shape, 

71 R0_waveform_high.mean(), 

72 R0_waveform_high.min(), 

73 R0_waveform_high.max(), 

74 ) 

75 print( 

76 "Waveform low: ", 

77 R0_waveform_low.shape, 

78 R0_waveform_low.mean(), 

79 R0_waveform_low.min(), 

80 R0_waveform_low.max(), 

81 ) 

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

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

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

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

86 R0_waveform, selected_gains, selected_pedestals, selected_correction = select_channel( 

87 R0_waveform_high, R0_waveform_low, gains, pedestals, gain_correction, 4000 

88 ) 

89 R0_calibrated = calibrate(R0_waveform, selected_gains, selected_pedestals) 

90 R0_integrated = integrate_local_peak(R0_calibrated, 3, 3) 

91 R0_integrated *= selected_correction 

92 

93 R0_integrated = R0_integrated.squeeze() 

94 

95 plotImage( 

96 pdf, 

97 camera_geometry[0], 

98 camera_geometry[1], 

99 R0_integrated, 

100 im_idx, 

101 tabEventId[im_idx], 

102 tabEventType[im_idx], 

103 ) 

104 

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

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

107 plotImage( 

108 pdf, 

109 camera_geometry[0], 

110 camera_geometry[1], 

111 np.abs(fft_image), 

112 im_idx, 

113 tabEventId[im_idx], 

114 tabEventType[im_idx], 

115 ) 

116 

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

118 fig.patch.set_alpha(1.0) 

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

120 plt.axis("equal") 

121 plt.xlabel("x") 

122 plt.ylabel("y") 

123 plt.colorbar() 

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

125 plt.close() 

126 

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

128 fig.patch.set_alpha(1.0) 

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

130 plt.axis("equal") 

131 plt.xlabel("x") 

132 plt.ylabel("y") 

133 plt.colorbar() 

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

135 plt.close() 

136 

137 

138if __name__ == "__main__": 

139 main()