Coverage for lst_auto_rta/merge_DL3_no_analysis.py: 0%
112 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-22 14:47 +0000
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-22 14:47 +0000
1#!/usr/bin/env python
2# S. Caroff
4"""
5Script to merge DL3 and Analyse them
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"""
15import argparse
16import glob
17import os
19# import matplotlib.colors as mcolors
20import warnings
21from pathlib import Path
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
31warnings.filterwarnings("ignore")
32warnings.simplefilter("ignore")
33# import toml
34import glob
35import time
37# from PIL import Image
40def CreatePlots(directory, run_id, RA, DEC, RA_map, DEC_map):
41 datastore = DataStore.from_dir(directory)
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()
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,)
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])
72def main():
73 parser = argparse.ArgumentParser(description="DL3 merger")
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 )
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 )
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()
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)
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)
143 # print(len(hdu1['EVENTS'].data))
144 hdu1["EVENTS"].data = hdu1["EVENTS"].data[np.argsort(hdu1["EVENTS"].data["TIME"])]
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"
179 hdu1["GTI"].data[0][0] = hdu1["EVENTS"].data[0][1]
180 hdu1["GTI"].data[0][1] = hdu1["EVENTS"].data[-1][1]
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")
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 + "/.")
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]))
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'])
219if __name__ == "__main__":
220 main()