Coverage for lst_auto_rta/merge_DL3_ACADA.py: 0%

1#!/usr/bin/env python 

2# S. Caroff 

3 

4""" 

5Script to merge DL3 and Analyse them 

6 

7Usage: 

8$> python merge_DL3.py 

9--input-filter "dl3_LST-1.Run04190.*.fits" 

10--input-directory "./DL3/" 

11--output-dir "./DL3/" 

12--run_id 9864 

13""" 

14 

15import argparse 

16import glob 

17import os 

18 

19# import matplotlib.colors as mcolors 

20import warnings 

21from pathlib import Path 

22 

23import astropy.units as u 

24import matplotlib.pyplot as plt 

25import numpy as np 

26from astropy.io import fits 

27from astropy.table import Column, QTable, Table, vstack 

28 

29# from SourceAnalyse import event 

30from gammapy.data import DataStore, EventList 

31 

32warnings.filterwarnings("ignore") 

33warnings.simplefilter("ignore") 

34# import toml 

35import glob 

36import time 

37 

38# from PIL import Image 

39 

40 

41def main(): 

42 parser = argparse.ArgumentParser( 

43 description="Automatic Script for the DL3 check", formatter_class=argparse.ArgumentDefaultsHelpFormatter 

44 ) 

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

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

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

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

49 parser.add_argument("-rf", "--run-fake", default="", help="local run id ") 

50 args = parser.parse_args() 

51 config = vars(args) 

52 

53 mypath = config["directory"] + config["date"] + "/" + config["run_id"] + "/" + config["add_string"] + "/dl3/" 

54 mypath_rundir = config["directory"] + config["date"] + "/" + config["run_id"] + "/" + config["add_string"] 

55 

56 file_filter = "dl3_v06_*.fits" # "dl3_LST-1.Run04190.*.fits" 

57 print(mypath + file_filter) 

58 filelist = glob.glob(mypath + file_filter) 

59 print(filelist) 

60 

61 hdu1 = fits.open(filelist[0]) 

62 for file in filelist[:]: 

63 hdu2 = fits.open(file) 

64 nev_hdu1 = len(hdu1["EVENTS"].data) 

65 nev_hdu2 = len(hdu2["EVENTS"].data) 

66 if (hdu1["EVENTS"].header["RA_PNT"] - hdu2["EVENTS"].header["RA_PNT"]) ** 2 < 2**2: 

67 print(hdu1["EVENTS"].header["RA_PNT"], " ", hdu2["EVENTS"].header["RA_PNT"]) 

68 print(file) 

69 hdu1["EVENTS"].header["RA_PNT"] = ( 

70 hdu1["EVENTS"].header["RA_PNT"] * nev_hdu1 + hdu2["EVENTS"].header["RA_PNT"] * nev_hdu2 

71 ) / (nev_hdu1 + nev_hdu2) 

72 hdu1["EVENTS"].header["DEC_PNT"] = ( 

73 hdu1["EVENTS"].header["DEC_PNT"] * nev_hdu1 + hdu2["EVENTS"].header["DEC_PNT"] * nev_hdu2 

74 ) / (nev_hdu1 + nev_hdu2) 

75 hdu1["EVENTS"].header["ALT_PNT"] = ( 

76 hdu1["EVENTS"].header["ALT_PNT"] * nev_hdu1 + hdu2["EVENTS"].header["ALT_PNT"] * nev_hdu2 

77 ) / (nev_hdu1 + nev_hdu2) 

78 hdu1["EVENTS"].header["AZ_PNT"] = ( 

79 hdu1["EVENTS"].header["AZ_PNT"] * nev_hdu1 + hdu2["EVENTS"].header["AZ_PNT"] * nev_hdu2 

80 ) / (nev_hdu1 + nev_hdu2) 

81 hdu1["EVENTS"].header["DEADC"] = ( 

82 hdu1["EVENTS"].header["DEADC"] * nev_hdu1 + hdu2["EVENTS"].header["DEADC"] * nev_hdu2 

83 ) / (nev_hdu1 + nev_hdu2) 

84 if not (hdu1["EVENTS"].header["OBS_ID"] == hdu2["EVENTS"].header["OBS_ID"]): 

85 print("BEWARE ! You are merging DL3 event with different OBS_ID... You should not do that...") 

86 # print(len(hdu1['EVENTS'].data)) 

87 hdu1["EVENTS"].data = np.append(hdu1["EVENTS"].data, hdu2["EVENTS"].data) 

88 

89 # print(len(hdu1['EVENTS'].data)) 

90 hdu1["EVENTS"].data = hdu1["EVENTS"].data[np.argsort(hdu1["EVENTS"].data["TIME"])] 

91 

92 hdu1["EVENTS"].data = hdu1["EVENTS"].data[ 

93 (hdu1["EVENTS"].data["TIME"] - hdu1["EVENTS"].data["TIME"][0]) < 8 * 3600 

94 ] 

95 # hdu1['EVENTS'].header['RA_PNT'] = np.mean(hdu1['EVENTS'].data['RA']) 

96 # hdu1['EVENTS'].header['DEC_PNT'] = np.mean(hdu1['EVENTS'].data['DEC']) 

97 # hdu1['EVENTS'].data['TIME'] = hdu1['EVENTS'].data['TIME'] 

98 # hdu1['EVENTS'].data['RA'] = hdu1['EVENTS'].data['RA'] 

99 # hdu1['EVENTS'].data['DEC'] = hdu1['EVENTS'].data['DEC'] 

100 # hdu1['EVENTS'].data['ENERGY'] = hdu1['EVENTS'].data['ENERGY'] 

101 # hdu1['EVENTS'].data[:][] = hdu1['EVENTS'].data[:][1]*u.s 

102 hdu1["EVENTS"].header["OBS_ID"] = config["run_fake"] 

103 hdu1["EVENTS"].header["N_TELS"] = 1 

104 hdu1["EVENTS"].header["TELLIST"] = "LST-1" 

105 hdu1["EVENTS"].header["TSTART"] = hdu1["EVENTS"].data["TIME"][0] 

106 hdu1["EVENTS"].header["TSTOP"] = hdu1["EVENTS"].data["TIME"][-1] 

107 hdu1["EVENTS"].header["ONTIME"] = hdu1["EVENTS"].data["TIME"][-1] - hdu1["EVENTS"].data["TIME"][0] 

108 hdu1["EVENTS"].header["LIVETIME"] = hdu1["EVENTS"].header["ONTIME"] * hdu1["EVENTS"].header["DEADC"] 

109 hdu1["EVENTS"].header["INSTRUME"] = "LST" 

110 hdu1["EVENTS"].header["DATE-OBS"] = 0 

111 hdu1["EVENTS"].header["TIME-OBS"] = 0 

112 hdu1["EVENTS"].header["DATE-END"] = 0 

113 hdu1["EVENTS"].header["TIME-END"] = 0 

114 hdu1["EVENTS"].header["TELAPSE"] = 0 

115 hdu1["POINTING"].header["ALT_PNT"] = hdu1["EVENTS"].header["ALT_PNT"] 

116 hdu1["POINTING"].header["AZ_PNT"] = hdu1["EVENTS"].header["AZ_PNT"] 

117 hdu1["POINTING"].header["TIME"] = hdu1["EVENTS"].header["TSTART"] 

118 prefix = (5 - len(str(hdu1["EVENTS"].header["OBS_ID"]))) * "0" 

119 # print(args.output_fits_dir+'dl3_LST-1.Run'+prefix+str(hdu1['EVENTS'].header['OBS_ID'])+'.fits') 

120 hdu1["EVENTS"].columns["TIME"].unit = "s" 

121 hdu1["EVENTS"].columns["ENERGY"].unit = "TeV" 

122 hdu1["EVENTS"].columns["RA"].unit = "deg" 

123 hdu1["EVENTS"].columns["DEC"].unit = "deg" 

124 

125 hdu1["GTI"].data[0][0] = hdu1["EVENTS"].data[0][1] 

126 hdu1["GTI"].data[0][1] = hdu1["EVENTS"].data[-1][1] 

127 

128 if not os.path.exists(mypath_rundir + "plots/"): 

129 os.mkdir(mypath_rundir + "plots/") 

130 

131 hdu1.writeto( 

132 mypath_rundir + "plots/dl3_LST-1.Run" + prefix + str(hdu1["EVENTS"].header["OBS_ID"]) + ".fits", overwrite=True 

133 ) 

134 print("create directory ") 

135 print(mypath_rundir + "plots/") 

136 # if not os.path.exists(mypath_rundir+'plots/'): 

137 # os.mkdir(mypath_rundir+'plots/') 

138 

139 files_and_directories = os.listdir(config["directory"] + config["date"] + "/") 

140 for item in files_and_directories: 

141 print( 

142 "cp -f " 

143 + config["directory"] 

144 + config["date"] 

145 + "/" 

146 + item 

147 + "/" 

148 + config["add_string"] 

149 + "/plots/dl3_LST* " 

150 + mypath_rundir 

151 + "plots/." 

152 ) 

153 os.system( 

154 "cp -f " 

155 + config["directory"] 

156 + config["date"] 

157 + "/" 

158 + item 

159 + "/" 

160 + config["add_string"] 

161 + "/plots/dl3_LST* " 

162 + mypath_rundir 

163 + "plots/." 

164 ) 

165 

166 filelist_runs = glob.glob( 

167 config["directory"] + config["date"] + "/" + config["run_id"] + "/" + config["add_string"] + "/plots/dl3_LST*" 

168 ) 

169 print("runlist is : ") 

170 runlist_to_analyse = [] 

171 for i in range(len(filelist_runs)): 

172 print(filelist_runs[i][-10:-5]) 

173 hdutemp = fits.open(filelist_runs[i]) 

174 if (hdutemp["EVENTS"].header["RA_PNT"] - hdu1["EVENTS"].header["RA_PNT"]) ** 2 < 9 and ( 

175 hdutemp["EVENTS"].header["DEC_PNT"] - hdu1["EVENTS"].header["DEC_PNT"] 

176 ) ** 2 < 9: 

177 runlist_to_analyse.append(int(filelist_runs[i][-10:-5])) 

178 

179 print( 

180 "lstchain_create_dl3_index_files --file-pattern='dl3_LST*' --input-dl3-dir='" 

181 + config["directory"] 

182 + config["date"] 

183 + "/" 

184 + config["run_id"] 

185 + "/" 

186 + config["add_string"] 

187 + "/plots/' --overwrite " 

188 ) 

189 os.system( 

190 "lstchain_create_dl3_index_files --file-pattern='dl3_LST*' --input-dl3-dir='" 

191 + config["directory"] 

192 + config["date"] 

193 + "/" 

194 + config["run_id"] 

195 + "/" 

196 + config["add_string"] 

197 + "/plots/' --overwrite" 

198 ) 

199 print(hdu1["EVENTS"].header["RA_PNT"]) 

200 print(" final runlist for analysis is : ") 

201 print(runlist_to_analyse) 

202 

203 

204if __name__ == "__main__": 

205 main()