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

21 

22import numpy as np 

23from astropy.io import fits 

24 

25# from SourceAnalyse import event 

26 

27warnings.filterwarnings("ignore") 

28warnings.simplefilter("ignore") 

29# import toml 

30 

31# from PIL import Image 

32 

33 

34def main(): 

35 parser = argparse.ArgumentParser( 

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

37 ) 

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

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

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

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

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

43 args = parser.parse_args() 

44 config = vars(args) 

45 

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

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

48 

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

50 print(mypath + file_filter) 

51 filelist = glob.glob(mypath + file_filter) 

52 print(filelist) 

53 

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

55 for file in filelist[:]: 

56 hdu2 = fits.open(file) 

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

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

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

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

61 print(file) 

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

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

64 ) / (nev_hdu1 + nev_hdu2) 

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

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

67 ) / (nev_hdu1 + nev_hdu2) 

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

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

70 ) / (nev_hdu1 + nev_hdu2) 

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

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

73 ) / (nev_hdu1 + nev_hdu2) 

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

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

76 ) / (nev_hdu1 + nev_hdu2) 

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

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

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

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

81 

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

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

84 

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

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

87 ] 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

117 

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

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

120 

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

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

123 

124 hdu1.writeto( 

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

126 overwrite=True, 

127 ) 

128 print("create directory ") 

129 print(mypath_rundir + "plots/") 

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

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

132 

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

134 for item in files_and_directories: 

135 print( 

136 "cp -f " 

137 + config["directory"] 

138 + config["date"] 

139 + "/" 

140 + item 

141 + "/" 

142 + config["add_string"] 

143 + "/plots/dl3_LST* " 

144 + mypath_rundir 

145 + "plots/." 

146 ) 

147 os.system( 

148 "cp -f " 

149 + config["directory"] 

150 + config["date"] 

151 + "/" 

152 + item 

153 + "/" 

154 + config["add_string"] 

155 + "/plots/dl3_LST* " 

156 + mypath_rundir 

157 + "plots/." 

158 ) 

159 

160 filelist_runs = glob.glob( 

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

162 ) 

163 print("runlist is : ") 

164 runlist_to_analyse = [] 

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

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

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

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

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

170 ) ** 2 < 9: 

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

172 

173 print( 

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

175 + config["directory"] 

176 + config["date"] 

177 + "/" 

178 + config["run_id"] 

179 + "/" 

180 + config["add_string"] 

181 + "/plots/' --overwrite " 

182 ) 

183 os.system( 

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

185 + config["directory"] 

186 + config["date"] 

187 + "/" 

188 + config["run_id"] 

189 + "/" 

190 + config["add_string"] 

191 + "/plots/' --overwrite" 

192 ) 

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

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

195 print(runlist_to_analyse) 

196 

197 

198if __name__ == "__main__": 

199 main()