Coverage for lst_auto_rta/Auto_Check_DL2.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 DL2 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") 

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

21# parser.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 + "DL2/" 

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

44 params = [] 

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

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

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

48 params = vstack(params) 

49 max_trigger_time = 1894367064 

50 mask = params["trigger_time"] < max_trigger_time # *(params['trigger_time']>1651613520) 

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_DL2.pdf") as pdf: 

59 plot_variable(params, pdf, "log_reco_energy") 

60 plot_variable(params, pdf, "reco_energy") 

61 plot_variable(params, pdf, "reco_disp_dx") 

62 plot_variable(params, pdf, "reco_disp_dy") 

63 plot_variable(params, pdf, "reco_src_x") 

64 plot_variable(params, pdf, "reco_src_y") 

65 plot_variable(params, pdf, "gammaness") 

66 

67 mask = params["is_good_event"] == 1 

68 params = params[mask] 

69 mask = params["gammaness"] > 0.7 

70 

71 plt.hist2d(params["reco_src_x"][mask], params["reco_src_y"][mask], bins=100) 

72 plt.xlabel("reco_src_x") 

73 plt.ylabel("reco_src_y") 

74 pdf.savefig() 

75 plt.close() 

76 

77 

78if __name__ == "__main__": 

79 main()