"""HiPeRTA file processing pipeline in python."""
import numpy as np
import tables
from pyhiperta.calibration import calibrate, select_channel
from pyhiperta.integration import integrate_local_peak, windowed_integration_correction
from pyhiperta.R0Reader import R0HDF5Dataset
from pyhiperta.waveform_indexing import waveform1Dto2D, waveform2Dto1D, waveform_1D_to_2D_maps
[docs]
def main():
import argparse
import matplotlib.pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
def plotImage(pdf, tabX, tabY, signal, imageIndex, eventId, eventType, tabFixPixelOrder=None):
strInfo = f"image {imageIndex}, id = {eventId}, type = {eventType},\n Signal mean = {signal.mean()} pe, min = {signal.min()}, max = {signal.max()}, std = {signal.std()}"
print(strInfo)
fig = plt.figure(figsize=(16, 10))
fig.patch.set_alpha(1.0)
if tabFixPixelOrder is not None:
plt.scatter(tabX[tabFixPixelOrder], tabY[tabFixPixelOrder], c=signal, s=120)
else:
plt.scatter(tabX, tabY, c=signal, s=120)
plt.axis("equal")
plt.xlabel("x")
plt.ylabel("y")
plt.colorbar()
plt.text(0.05, 0.95, strInfo, transform=fig.transFigure, size=14)
pdf.savefig() # saves the current figure into a pdf page
plt.close()
parser = argparse.ArgumentParser()
parser.add_argument("--R0_files", "-r0", type=str, required=True, help="input file, or folder.")
parser.add_argument("--output_file", "-o", type=str, required=True, help="output path of produced pdf.")
parser.add_argument("--first_image", "-f", type=int, default=1, help="Index of first image to plot.")
parser.add_argument("--last_image", "-l", type=int, default=100, help="Index of last image to plot.")
parser.add_argument(
"--indexing_step", "-idx", type=int, default=1, help="Step between the images indices selected for plotting."
)
args = parser.parse_args()
dataset = R0HDF5Dataset(args.R0_files, 1)
gains = dataset.read_gains()
pedestals = dataset.read_pedestals()
pedestals /= dataset.events_n_frames() # pedestal per frame option of the poor
camera_geometry = dataset.read_camera_geometry()
map_1D_to_2D = waveform_1D_to_2D_maps(dataset.read_camera_geometry())
with tables.open_file(args.R0_files, "r") as hfile:
tabEventId = hfile.root.r0.event.telescope.waveform.tel_001.col("event_id")
tabEventType = hfile.root.r0.event.subarray.trigger.col("event_type")
# tabPixelOrder = hfile.root.configuration.instrument.telescope.camera.pixel_order.tel_001.read()
ref_pulse_sample_times, ref_pulse_shape = dataset.read_reference_pulse()
gain_correction = windowed_integration_correction(
ref_pulse_sample_times, ref_pulse_shape, 1, lambda x: integrate_local_peak(x, 3, 3)
)
print("gain correction: ", gain_correction)
with PdfPages(args.output_file) as pdf:
for im_idx in range(args.first_image, args.last_image, args.indexing_step):
R0_waveform_high, R0_waveform_low = dataset[im_idx]
print(
"Waveform high: ",
R0_waveform_high.shape,
R0_waveform_high.mean(),
R0_waveform_high.min(),
R0_waveform_high.max(),
)
print(
"Waveform low: ",
R0_waveform_low.shape,
R0_waveform_low.mean(),
R0_waveform_low.min(),
R0_waveform_low.max(),
)
print("gains high: ", gains[0].shape, gains[0].mean(), gains[0].min(), gains[0].max())
print("gains low: ", gains[1].shape, gains[1].mean(), gains[1].min(), gains[1].max())
print("pedestals high: ", pedestals[0].shape, pedestals[0].mean(), pedestals[0].min(), pedestals[0].max())
print("pedestals low: ", pedestals[1].shape, pedestals[1].mean(), pedestals[1].min(), pedestals[1].max())
R0_waveform, selected_gains, selected_pedestals, selected_correction = select_channel(
R0_waveform_high, R0_waveform_low, gains, pedestals, gain_correction, 4000
)
R0_calibrated = calibrate(R0_waveform, selected_gains, selected_pedestals)
R0_integrated = integrate_local_peak(R0_calibrated, 3, 3)
R0_integrated *= selected_correction
R0_integrated = R0_integrated.squeeze()
plotImage(
pdf,
camera_geometry[0],
camera_geometry[1],
R0_integrated,
im_idx,
tabEventId[im_idx],
tabEventType[im_idx],
)
fft_image = waveform2Dto1D(map_1D_to_2D, np.fft.fft2(waveform1Dto2D(map_1D_to_2D, R0_integrated)))
print("fft_image: ", fft_image.shape, fft_image.mean(), fft_image.min(), fft_image.max())
plotImage(
pdf,
camera_geometry[0],
camera_geometry[1],
np.abs(fft_image),
im_idx,
tabEventId[im_idx],
tabEventType[im_idx],
)
fig = plt.figure(figsize=(16, 10))
fig.patch.set_alpha(1.0)
plt.imshow(waveform1Dto2D(map_1D_to_2D, R0_integrated))
plt.axis("equal")
plt.xlabel("x")
plt.ylabel("y")
plt.colorbar()
pdf.savefig() # saves the current figure into a pdf page
plt.close()
fig = plt.figure(figsize=(16, 10))
fig.patch.set_alpha(1.0)
plt.imshow(np.abs(np.fft.fft2(waveform1Dto2D(map_1D_to_2D, R0_integrated))))
plt.axis("equal")
plt.xlabel("x")
plt.ylabel("y")
plt.colorbar()
pdf.savefig() # saves the current figure into a pdf page
plt.close()
if __name__ == "__main__":
main()