Coverage for src / lstautorta / merge_DL3.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 logging 

18import os 

19import sys 

20import warnings 

21from pathlib import Path 

22 

23import numpy as np 

24from astropy.io import fits 

25 

26warnings.filterwarnings("ignore") 

27warnings.simplefilter("ignore") 

28 

29 

30def setup_logging(output_dir, run_id): 

31 """ 

32 Configure logging with one log file per execution. 

33 

34 If output_dir does not exist, write the log in a fallback location 

35 and raise FileNotFoundError so the script crashes immediately. 

36 """ 

37 logger = logging.getLogger("merge_DL3") 

38 logger.setLevel(logging.INFO) 

39 logger.handlers.clear() 

40 

41 formatter = logging.Formatter( 

42 fmt="%(asctime)s | %(levelname)s | %(message)s", 

43 datefmt="%Y-%m-%d %H:%M:%S", 

44 ) 

45 

46 # Always keep stdout logging 

47 stream_handler = logging.StreamHandler(sys.stdout) 

48 stream_handler.setLevel(logging.INFO) 

49 stream_handler.setFormatter(formatter) 

50 logger.addHandler(stream_handler) 

51 

52 output_path = Path(output_dir) 

53 

54 if not output_path.exists(): 

55 fallback_log = Path.cwd() / f"merge_DL3_run{run_id}.log" 

56 

57 file_handler = logging.FileHandler(fallback_log) 

58 file_handler.setLevel(logging.INFO) 

59 file_handler.setFormatter(formatter) 

60 logger.addHandler(file_handler) 

61 

62 logger.info("============================================================") 

63 logger.error("Output directory does not exist: %s", output_path) 

64 logger.error("Fallback log file: %s", fallback_log) 

65 logger.error("Crashing because output directory is missing.") 

66 logger.info("============================================================") 

67 

68 raise FileNotFoundError(f"Output directory does not exist: {output_path}") 

69 

70 log_dir = output_path / "logs" 

71 log_dir.mkdir(exist_ok=True) 

72 

73 log_file = log_dir / f"merge_DL3_run{run_id}.log" 

74 

75 file_handler = logging.FileHandler(log_file) 

76 file_handler.setLevel(logging.INFO) 

77 file_handler.setFormatter(formatter) 

78 logger.addHandler(file_handler) 

79 

80 logger.info("============================================================") 

81 logger.info("Starting merge_DL3.py") 

82 logger.info("Log file: %s", log_file) 

83 

84 return logger, log_file 

85 

86 

87def main(): 

88 parser = argparse.ArgumentParser(description="DL3 merger") 

89 

90 parser.add_argument( 

91 "--input-filter", 

92 "-f", 

93 type=str, 

94 dest="input_data", 

95 help='DL3 filter selection, example: "dl3_LST-1.Run04190.*.fits"', 

96 default=None, 

97 required=True, 

98 ) 

99 

100 parser.add_argument( 

101 "--input-directory", 

102 "-d", 

103 type=str, 

104 dest="input_dir", 

105 help="Path to input DL3 files directory", 

106 default=None, 

107 required=True, 

108 ) 

109 

110 parser.add_argument( 

111 "--output-dir", 

112 "-o", 

113 type=str, 

114 dest="output_fits_dir", 

115 help="Path to output files", 

116 default=None, 

117 required=True, 

118 ) 

119 

120 parser.add_argument( 

121 "--run_id", 

122 "-r", 

123 type=int, 

124 dest="run_id", 

125 help="Run ID to assign to merged file", 

126 default=None, 

127 required=True, 

128 ) 

129 

130 args = parser.parse_args() 

131 

132 logger, log_file = setup_logging(args.output_fits_dir, args.run_id) 

133 

134 try: 

135 logger.info("Arguments:") 

136 logger.info(" input_filter = %s", args.input_data) 

137 logger.info(" input_dir = %s", args.input_dir) 

138 logger.info(" output_fits_dir= %s", args.output_fits_dir) 

139 logger.info(" run_id = %s", args.run_id) 

140 

141 file_filter = args.input_data 

142 directory = args.input_dir 

143 

144 search_pattern = directory + file_filter 

145 logger.info("Searching files with pattern: %s", search_pattern) 

146 filelist = glob.glob(search_pattern) 

147 filelist = sorted(filelist) 

148 

149 logger.info("Found %d input files.", len(filelist)) 

150 for f in filelist: 

151 logger.info(" input file: %s", f) 

152 

153 if len(filelist) == 0: 

154 logger.error("No input files found. Exiting.") 

155 sys.exit(1) 

156 

157 logger.info("Opening first FITS file: %s", filelist[0]) 

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

159 

160 for file in filelist[1:]: 

161 logger.info("Processing file: %s", file) 

162 hdu2 = fits.open(file) 

163 

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

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

166 

167 logger.info( 

168 "Current merge status: nev_hdu1=%d, nev_hdu2=%d, RA1=%.6f, RA2=%.6f", 

169 nev_hdu1, 

170 nev_hdu2, 

171 hdu1["EVENTS"].header["RA_PNT"], 

172 hdu2["EVENTS"].header["RA_PNT"], 

173 ) 

174 

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

176 logger.info("File accepted for merge: %s", file) 

177 

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

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

180 ) / (nev_hdu1 + nev_hdu2) 

181 

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

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

184 ) / (nev_hdu1 + nev_hdu2) 

185 

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

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

188 ) / (nev_hdu1 + nev_hdu2) 

189 

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

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

192 ) / (nev_hdu1 + nev_hdu2) 

193 

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

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

196 ) / (nev_hdu1 + nev_hdu2) 

197 

198 if hdu1["EVENTS"].header["OBS_ID"] != hdu2["EVENTS"].header["OBS_ID"]: 

199 logger.warning( 

200 "Merging DL3 events with different OBS_ID values: %s vs %s", 

201 hdu1["EVENTS"].header["OBS_ID"], 

202 hdu2["EVENTS"].header["OBS_ID"], 

203 ) 

204 

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

206 logger.info("Merge done. New event count: %d", len(hdu1["EVENTS"].data)) 

207 else: 

208 logger.info("File rejected by RA_PNT criterion: %s", file) 

209 

210 hdu2.close() 

211 

212 logger.info("Sorting merged EVENTS by TIME...") 

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

214 

215 logger.info("Applying 8-hour time window selection...") 

216 initial_count = len(hdu1["EVENTS"].data) 

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

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

219 ] 

220 logger.info( 

221 "Time selection kept %d / %d events.", 

222 len(hdu1["EVENTS"].data), 

223 initial_count, 

224 ) 

225 

226 logger.info("Updating FITS headers...") 

227 hdu1["EVENTS"].header["OBS_ID"] = args.run_id 

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

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

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

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

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

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

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

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

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

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

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

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

240 

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

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

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

244 

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

246 

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

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

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

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

251 

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

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

254 

255 output_file = args.output_fits_dir + "dl3_LST-1.Run" + prefix + str(hdu1["EVENTS"].header["OBS_ID"]) + ".fits" 

256 

257 logger.info("Writing merged FITS file to: %s", output_file) 

258 hdu1.writeto(output_file, overwrite=True) 

259 logger.info("Merged FITS file written successfully.") 

260 

261 plots_dir = directory[:-10] + str(int(args.run_id)) + "/plots" 

262 logger.info("Checking plots directory: %s", plots_dir) 

263 if not os.path.exists(plots_dir): 

264 os.mkdir(plots_dir) 

265 logger.info("Created plots directory: %s", plots_dir) 

266 else: 

267 logger.info("Plots directory already exists.") 

268 

269 logger.info("Copying previous DL3 files from the last 150 runs...") 

270 for i in range(150): 

271 src_cmd = "cp -f " + directory[:-10] + str(int(args.run_id) - i) + "/DL3/dl3_LST* " + directory + "/." 

272 logger.info("Executing command: %s", src_cmd) 

273 ret = os.system(src_cmd) 

274 logger.info("Command exit code: %s", ret) 

275 

276 filelist_runs = glob.glob(directory + "dl3_LST*") 

277 logger.info("Found %d DL3 files in working directory after copy.", len(filelist_runs)) 

278 

279 runlist_to_analyse = [] 

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

281 logger.info("Inspecting candidate run file: %s", filelist_runs[i]) 

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

283 

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

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

286 ) ** 2 < 9: 

287 run_id_candidate = int(filelist_runs[i][-10:-5]) 

288 runlist_to_analyse.append(run_id_candidate) 

289 logger.info("Run selected for analysis: %s", run_id_candidate) 

290 

291 hdutemp.close() 

292 

293 cmd_index = ( 

294 "lstchain_create_dl3_index_files --file-pattern='dl3_LST*' " 

295 "--input-dl3-dir='" + directory + "' --overwrite " 

296 ) 

297 logger.info("Executing index command: %s", cmd_index) 

298 ret = os.system(cmd_index) 

299 logger.info("Index command exit code: %s", ret) 

300 

301 logger.info("Final RA_PNT: %s", hdu1["EVENTS"].header["RA_PNT"]) 

302 logger.info("Final run list for analysis: %s", runlist_to_analyse) 

303 

304 # Uncomment if needed 

305 # RA, DEC = CreatePlots( 

306 # directory, 

307 # [runlist_to_analyse[-1]], 

308 # hdu1["EVENTS"].header["RA_PNT"], 

309 # hdu1["EVENTS"].header["DEC_PNT"], 

310 # hdu1["EVENTS"].header["RA_PNT"], 

311 # hdu1["EVENTS"].header["DEC_PNT"], 

312 # logger=logger, 

313 # ) 

314 

315 hdu1.close() 

316 

317 logger.info("merge_DL3.py finished successfully.") 

318 logger.info("============================================================") 

319 

320 except Exception as e: 

321 logger.exception("Fatal error while running merge_DL3.py: %s", e) 

322 logger.info("============================================================") 

323 raise 

324 

325 

326if __name__ == "__main__": 

327 main()