Coverage for lst_auto_rta/Auto_Check_DL1.py: 0%

1#!/usr/bin/env python3 

2import argparse 

3import os 

4from os import listdir 

5from os.path import isfile, join 

6 

7import matplotlib.pyplot as plt 

8import numpy as np 

9from astropy.table import Table, vstack 

10from ctapipe.containers import EventType 

11from ctapipe.io import read_table 

12from matplotlib.backends.backend_pdf import PdfPages 

13 

14parser = argparse.ArgumentParser( 

15 description="Automatic Script for the DL1 check", formatter_class=argparse.ArgumentDefaultsHelpFormatter 

16) 

17parser.add_argument("-d", "--directory", default="/fefs/onsite/pipeline/rta/data/", help="Directory for data") 

18parser.add_argument("-da", "--date", default="20230705", help="Date of the run to check") 

19parser.add_argument("-r", "--run-id", default="13600", help="run id to check") 

20#parser.add_argument("-add", "--add-string", default="reco/", help="add a string to the path") 

21parser.add_argument("-add", "--add-string", default="" , help="add a string to the path") 

22 

23args = parser.parse_args() 

24config = vars(args) 

25 

26 

27def plot_variable(table, pdf, name): 

28 print(name) 

29 mask = table["is_good_event"] == 1 

30 table = table[mask] 

31 plt.hist(table[name], bins=100) 

32 plt.xlabel(name) 

33 pdf.savefig() 

34 plt.close() 

35 # plt.show() 

36 

37 

38def main(): 

39 mypath_run_dir = config["directory"] + config["date"] + "/" + config["run_id"] + "/" + config["add_string"] 

40 if not os.path.exists(mypath_run_dir + "plots"): 

41 os.mkdir(mypath_run_dir + "plots") 

42 mypath = mypath_run_dir + "DL1/" 

43 filename = [f for f in listdir(mypath) if isfile(join(mypath, f)) if "dl1_v06_" in f] 

44 params = [] 

45 for i in range(len(filename)): 

46 tablename = "/dl1/event/telescope/parameters/LST_LSTCam" 

47 params.append(read_table(mypath + filename[i], tablename)) 

48 params = vstack(params) 

49 max_trigger_time = 3487407909 

50 mask = params["trigger_time"] < max_trigger_time 

51 

52 start_time = np.sort(params["trigger_time"][mask])[0] 

53 end_time = np.sort(params["trigger_time"][mask])[-1] 

54 times = end_time - start_time 

55 print(times) 

56 nbins = int(times) 

57 # filter_ = np.diff(params['event_id'][mask]) < 10000 

58 with PdfPages(mypath_run_dir + "plots/output_DL1.pdf") as pdf: 

59 plt.plot((np.sort(params["trigger_time"][mask][:-1])) - start_time, np.diff(np.sort(params["event_id"][mask]))) 

60 plt.xlabel("Time (s)") 

61 plt.ylabel("Delta event id") 

62 pdf.savefig() 

63 plt.close() 

64 print("Delta event id") 

65 plt.hist(np.diff(np.sort(params["event_id"][mask])), bins=100) 

66 plt.yscale("log") 

67 plt.xlabel("Delta event id") 

68 pdf.savefig() 

69 plt.close() 

70 print("Delta event id Histo") 

71 

72 plt.hist(params["trigger_time"][mask] - start_time, bins=nbins, alpha=0.5, label="All") 

73 plt.xlabel("Time (s)") 

74 plt.ylabel("Rate (Hz)") 

75 

76 plt.legend() 

77 pdf.savefig() 

78 plt.close() 

79 print("Rate versus time") 

80 

81 plt.hist(params["is_good_event"], bins=100) 

82 plt.xlabel("Is Good Event") 

83 pdf.savefig() 

84 plt.close() 

85 print("Is Good Event") 

86 

87 plt.hist(params["event_type"], bins=32) 

88 plt.xlabel("event_type") 

89 pdf.savefig() 

90 plt.close() 

91 print("event_type") 

92 

93 plt.hist(np.log10(params["total_intensity"][params["event_quality"] == 0]), alpha=0.5, bins=100) 

94 plt.yscale("log") 

95 plt.xlabel("total_intensity") 

96 pdf.savefig() 

97 plt.close() 

98 print("total_intensity") 

99 

100 plt.hist(np.log10(params["intensity"][params["event_quality"] == 0]), alpha=0.5, bins=100) 

101 plt.yscale("log") 

102 plt.xlabel("intensity") 

103 pdf.savefig() 

104 plt.close() 

105 print("intensity") 

106 

107 plt.hist(params["event_quality"], bins=8) 

108 plt.xlabel("Event Quality") 

109 pdf.savefig() 

110 plt.close() 

111 print("event_quality") 

112 

113 plt.plot(params["trigger_time"][mask] - start_time, np.rad2deg(params["az_tel"][mask])) # [mask2])) 

114 plt.xlabel("Time (s)") 

115 plt.ylabel("az_tel (deg)") 

116 pdf.savefig() 

117 plt.close() 

118 print("az_tel") 

119 

120 plt.plot(params["trigger_time"][mask] - start_time, np.rad2deg(params["alt_tel"][mask])) # [mask2])) 

121 plt.xlabel("Time (s)") 

122 plt.ylabel("alt_tel (deg)") 

123 pdf.savefig() 

124 plt.close() 

125 print("alt_tel") 

126 

127 plot_variable(params, pdf, "x") 

128 plot_variable(params, pdf, "y") 

129 plot_variable(params, pdf, "phi") 

130 plot_variable(params, pdf, "width") 

131 plot_variable(params, pdf, "length") 

132 plot_variable(params, pdf, "intensity") 

133 plot_variable(params, pdf, "n_pixels") 

134 plot_variable(params, pdf, "skewness") 

135 plot_variable(params, pdf, "r") 

136 plot_variable(params, pdf, "kurtosis") 

137 plot_variable(params, pdf, "psi") 

138 plot_variable(params, pdf, "time_gradient") 

139 plot_variable(params, pdf, "intercept") 

140 plot_variable(params, pdf, "leakage_pixels_width_1") 

141 plot_variable(params, pdf, "leakage_pixels_width_2") 

142 plot_variable(params, pdf, "leakage_intensity_width_1") 

143 plot_variable(params, pdf, "leakage_intensity_width_2") 

144 plot_variable(params, pdf, "alt_tel") 

145 plot_variable(params, pdf, "az_tel") 

146 plot_variable(params, pdf, "trigger_time") 

147 # plot_variable(params,pdf,"log_intensity") 

148 plot_variable(params, pdf, "wl") 

149 

150 

151if __name__ == "__main__": 

152 main()