Coverage for src / lstautorta / 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 vstack 

10from ctapipe.io import read_table 

11from matplotlib.backends.backend_pdf import PdfPages 

12 

13parser = argparse.ArgumentParser( 

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

15) 

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

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

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

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

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

21 

22args = parser.parse_args() 

23config = vars(args) 

24 

25 

26def plot_variable(table, pdf, name): 

27 print(name) 

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

29 table = table[mask] 

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

31 plt.xlabel(name) 

32 pdf.savefig() 

33 plt.close() 

34 # plt.show() 

35 

36 

37def main(): 

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

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

40 os.mkdir(mypath_run_dir + "plots") 

41 mypath = mypath_run_dir + "DL1/" 

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

43 params = [] 

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

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

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

47 params = vstack(params) 

48 max_trigger_time = 3487407909 

49 mask = params["trigger_time"] < max_trigger_time 

50 

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

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

53 times = end_time - start_time 

54 print(times) 

55 nbins = int(times) 

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

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

58 plt.plot( 

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

60 ) 

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

62 plt.ylabel("Delta event id") 

63 pdf.savefig() 

64 plt.close() 

65 print("Delta event id") 

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

67 plt.yscale("log") 

68 plt.xlabel("Delta event id") 

69 pdf.savefig() 

70 plt.close() 

71 print("Delta event id Histo") 

72 

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

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

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

76 

77 plt.legend() 

78 pdf.savefig() 

79 plt.close() 

80 print("Rate versus time") 

81 

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

83 plt.xlabel("Is Good Event") 

84 pdf.savefig() 

85 plt.close() 

86 print("Is Good Event") 

87 

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

89 plt.xlabel("event_type") 

90 pdf.savefig() 

91 plt.close() 

92 print("event_type") 

93 

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

95 plt.yscale("log") 

96 plt.xlabel("total_intensity") 

97 pdf.savefig() 

98 plt.close() 

99 print("total_intensity") 

100 

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

102 plt.yscale("log") 

103 plt.xlabel("intensity") 

104 pdf.savefig() 

105 plt.close() 

106 print("intensity") 

107 

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

109 plt.xlabel("Event Quality") 

110 pdf.savefig() 

111 plt.close() 

112 print("event_quality") 

113 

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

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

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

117 pdf.savefig() 

118 plt.close() 

119 print("az_tel") 

120 

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

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

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

124 pdf.savefig() 

125 plt.close() 

126 print("alt_tel") 

127 

128 plot_variable(params, pdf, "x") 

129 plot_variable(params, pdf, "y") 

130 plot_variable(params, pdf, "phi") 

131 plot_variable(params, pdf, "width") 

132 plot_variable(params, pdf, "length") 

133 plot_variable(params, pdf, "intensity") 

134 plot_variable(params, pdf, "n_pixels") 

135 plot_variable(params, pdf, "skewness") 

136 plot_variable(params, pdf, "r") 

137 plot_variable(params, pdf, "kurtosis") 

138 plot_variable(params, pdf, "psi") 

139 plot_variable(params, pdf, "time_gradient") 

140 plot_variable(params, pdf, "intercept") 

141 plot_variable(params, pdf, "leakage_pixels_width_1") 

142 plot_variable(params, pdf, "leakage_pixels_width_2") 

143 plot_variable(params, pdf, "leakage_intensity_width_1") 

144 plot_variable(params, pdf, "leakage_intensity_width_2") 

145 plot_variable(params, pdf, "alt_tel") 

146 plot_variable(params, pdf, "az_tel") 

147 plot_variable(params, pdf, "trigger_time") 

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

149 plot_variable(params, pdf, "wl") 

150 

151 

152if __name__ == "__main__": 

153 main()