Coverage for lst_auto_rta/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 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 

28from gammapy.data import DataStore, EventList 

29from SourceAnalyse import event 

30 

31warnings.filterwarnings("ignore") 

32warnings.simplefilter("ignore") 

33# import toml 

34import glob 

35import time 

36 

37# from PIL import Image 

38 

39 

40def CreatePlots(directory, run_id, RA, DEC, RA_map, DEC_map): 

41 datastore = DataStore.from_dir(directory) 

42 

43 Analysis = event(directory, run_id) 

44 print("Create On Region...") 

45 Analysis.create_on_region(RA, DEC) 

46 print("Create Exclusion Mask") 

47 Analysis.create_exclusion_mask() 

48 print("create reduction chain") 

49 Analysis.reduction_chain() 

50 print("create data maker") 

51 Analysis.data_maker() 

52 print("data run") 

53 # Analysis.data_run() 

54 # Analysis.signal_info() 

55 # Analysis.fit_spectrum() 

56 

57 # Analysis.plot_timee() 

58 # excess_reflected,excess_reflected_uncertainty,c,d,sig_reflected,sig_livetime = Analysis.signal_info() 

59 # Analysis.flux_plot(plot=False,) 

60 

61 print("Compute Acceptance model") 

62 Analysis.calculate_acceptance_model(plot=False) 

63 print("Add acceptance map") 

64 Analysis.add_acceptance_map() 

65 print("Compute map") 

66 Analysis.map_geometry(RA_map, DEC_map, plot=False, binsize=0.05, npix_x=150, npix_y=150) 

67 ring_excess, ring_sqrt_ts, ring_bg, uncertainty_excess_ring = Analysis.ring_data_estimation() 

68 sources = Analysis.excess_significance_plot(directory=directory, plot=True) 

69 return (sources[0][3], sources[0][4]) 

70 

71 

72def main(): 

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

74 

75 # Required arguments 

76 parser.add_argument( 

77 "--input-filter", 

78 "-f", 

79 type=str, 

80 dest="input_data", 

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

82 default=None, 

83 required=True, 

84 ) 

85 

86 parser.add_argument( 

87 "--input-directory", 

88 "-d", 

89 type=str, 

90 dest="input_dir", 

91 help="path to input DL3 files directory", 

92 default=None, 

93 required=True, 

94 ) 

95 

96 parser.add_argument( 

97 "--output-dir", 

98 "-o", 

99 type=str, 

100 dest="output_fits_dir", 

101 help="path to output files", 

102 default=None, 

103 required=True, 

104 ) 

105 parser.add_argument( 

106 "--run_id", "-r", type=int, dest="run_id", help="path to output files", default=None, required=True 

107 ) 

108 args = parser.parse_args() 

109 

110 file_filter = args.input_data # "dl3_LST-1.Run04190.*.fits" 

111 directory = args.input_dir # "/home/sami.caroff/cta-lstchain-enrique/cta-lstchain/lstchain/scripts/" 

112 print(directory + file_filter) 

113 filelist = glob.glob(directory + file_filter) 

114 

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

116 for file in filelist[1:]: 

117 hdu2 = fits.open(file) 

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

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

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

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

122 print(file) 

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

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

125 ) / (nev_hdu1 + nev_hdu2) 

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

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

128 ) / (nev_hdu1 + nev_hdu2) 

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

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

131 ) / (nev_hdu1 + nev_hdu2) 

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

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

134 ) / (nev_hdu1 + nev_hdu2) 

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

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

137 ) / (nev_hdu1 + nev_hdu2) 

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

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

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

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

142 

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

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

145 

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

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

148 ] 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

178 

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

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

181 

182 hdu1.writeto( 

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

184 overwrite=True, 

185 ) 

186 print("create directory ") 

187 print(directory[:-10] + "plots") 

188 if not os.path.exists(directory[:-10] + str(int(args.run_id)) + "/plots"): 

189 os.mkdir(directory[:-10] + str(int(args.run_id)) + "/plots") 

190 

191 for i in range(150): 

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

193 os.system("cp -f " + directory[:-10] + str(int(args.run_id) - i) + "/DL3/dl3_LST* " + directory + "/.") 

194 

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

196 print("runlist is : ") 

197 runlist_to_analyse = [] 

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

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

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

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

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

203 ) ** 2 < 9: 

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

205 

206 print("lstchain_create_dl3_index_files --file-pattern='dl3_LST*' --input-dl3-dir='" + directory + "' --overwrite ") 

207 os.system( 

208 "lstchain_create_dl3_index_files --file-pattern='dl3_LST*' --input-dl3-dir='" + directory + "' --overwrite " 

209 ) 

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

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

212 print(runlist_to_analyse) 

213 # RA,DEC = CreatePlots(directory,[runlist_to_analyse[-1]],hdu1['EVENTS'].header['RA_PNT'],hdu1['EVENTS'].header['DEC_PNT'],hdu1['EVENTS'].header['RA_PNT'],hdu1['EVENTS'].header['DEC_PNT']) 

214 # RA,DEC = CreatePlots(directory,[runlist_to_analyse[-1]],RA,DEC,hdu1['EVENTS'].header['RA_PNT'],hdu1['EVENTS'].header['DEC_PNT']) 

215 # RA,DEC = CreatePlots(directory,runlist_to_analyse,hdu1['EVENTS'].header['RA_PNT'],hdu1['EVENTS'].header['DEC_PNT'],hdu1['EVENTS'].header['RA_PNT'],hdu1['EVENTS'].header['DEC_PNT']) 

216 # CreatePlots(directory,runlist_to_analyse,RA,DEC,hdu1['EVENTS'].header['RA_PNT'],hdu1['EVENTS'].header['DEC_PNT']) 

217 

218 

219if __name__ == "__main__": 

220 main()