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

21 

22import numpy as np 

23from astropy.io import fits 

24from gammapy.data import DataStore 

25from SourceAnalyse import event 

26 

27warnings.filterwarnings("ignore") 

28warnings.simplefilter("ignore") 

29# import toml 

30 

31# from PIL import Image 

32 

33 

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

35 datastore = DataStore.from_dir(directory) 

36 

37 Analysis = event(directory, run_id) 

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

39 Analysis.create_on_region(RA, DEC) 

40 print("Create Exclusion Mask") 

41 Analysis.create_exclusion_mask() 

42 print("create reduction chain") 

43 Analysis.reduction_chain() 

44 print("create data maker") 

45 Analysis.data_maker() 

46 print("data run") 

47 # Analysis.data_run() 

48 # Analysis.signal_info() 

49 # Analysis.fit_spectrum() 

50 

51 # Analysis.plot_timee() 

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

53 # Analysis.flux_plot(plot=False,) 

54 

55 print("Compute Acceptance model") 

56 Analysis.calculate_acceptance_model(plot=False) 

57 print("Add acceptance map") 

58 Analysis.add_acceptance_map() 

59 print("Compute map") 

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

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

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

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

64 

65 

66def main(): 

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

68 

69 # Required arguments 

70 parser.add_argument( 

71 "--input-filter", 

72 "-f", 

73 type=str, 

74 dest="input_data", 

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

76 default=None, 

77 required=True, 

78 ) 

79 

80 parser.add_argument( 

81 "--input-directory", 

82 "-d", 

83 type=str, 

84 dest="input_dir", 

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

86 default=None, 

87 required=True, 

88 ) 

89 

90 parser.add_argument( 

91 "--output-dir", 

92 "-o", 

93 type=str, 

94 dest="output_fits_dir", 

95 help="path to output files", 

96 default=None, 

97 required=True, 

98 ) 

99 parser.add_argument( 

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

101 ) 

102 args = parser.parse_args() 

103 

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

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

106 print(directory + file_filter) 

107 filelist = glob.glob(directory + file_filter) 

108 

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

110 for file in filelist[1:]: 

111 hdu2 = fits.open(file) 

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

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

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

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

116 print(file) 

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

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

119 ) / (nev_hdu1 + nev_hdu2) 

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

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

122 ) / (nev_hdu1 + nev_hdu2) 

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

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

125 ) / (nev_hdu1 + nev_hdu2) 

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

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

128 ) / (nev_hdu1 + nev_hdu2) 

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

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

131 ) / (nev_hdu1 + nev_hdu2) 

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

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

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

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

136 

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

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

139 

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

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

142 ] 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

172 

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

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

175 

176 hdu1.writeto( 

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

178 overwrite=True, 

179 ) 

180 print("create directory ") 

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

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

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

184 

185 for i in range(150): 

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

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

188 

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

190 print("runlist is : ") 

191 runlist_to_analyse = [] 

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

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

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

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

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

197 ) ** 2 < 9: 

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

199 

200 print( 

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

202 ) 

203 os.system( 

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

205 ) 

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

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

208 print(runlist_to_analyse) 

209 # 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']) 

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

211 # 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']) 

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

213 

214 

215if __name__ == "__main__": 

216 main()