Coverage for lst_auto_rta/Theta_square.py: 0%
181 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-11-03 14:47 +0000
1#!/usr/bin/env python
3import astropy
4import gammapy
5import matplotlib
6import numpy as np
7import regions
9print("gammapy:", gammapy.__version__)
10print("numpy:", np.__version__)
11print("astropy", astropy.__version__)
12print("regions", regions.__version__)
13print("matplotlib", matplotlib.__version__)
15import os
17import astropy.units as u
18import matplotlib.pyplot as plt
19import matplotlib.style as style
20import numpy as np
21from astropy.coordinates import SkyCoord
23style.use("tableau-colorblind10")
24import argparse
25from pathlib import Path
27from gammapy.data import DataStore
28from gammapy.datasets import (
29 Datasets,
30 MapDataset,
31 SpectrumDataset,
32)
33from gammapy.makers import (
34 MapDatasetMaker,
35 ReflectedRegionsBackgroundMaker,
36 SafeMaskMaker,
37 SpectrumDatasetMaker,
38)
39from gammapy.makers.utils import make_theta_squared_table
40from gammapy.maps import Map, MapAxis, RegionGeom, WcsGeom
41from gammapy.visualization import plot_spectrum_datasets_off_regions, plot_theta_squared_table
42from matplotlib.offsetbox import AnchoredText
43from regions import CircleSkyRegion
44from scipy.optimize import curve_fit
46parser = argparse.ArgumentParser(
47 description="Automatic Script for the DL1 check", formatter_class=argparse.ArgumentDefaultsHelpFormatter
48)
49parser.add_argument("-d", "--directory", default="/fefs/onsite/pipeline/rta/data/", help="Directory for data")
50parser.add_argument("-da", "--date", default="20230705", help="Date of the run to check")
51parser.add_argument("-r", "--run-id", default="13600", help="run id to check")
52parser.add_argument("-add", "--add-string", default="" , help="add a string to the path")
53#parser.add_argument("-add", "--add-string", default="reco/", help="add a string to the path")
54parser.add_argument("-RA", "--right-ascension", default="270.19042", help="right-ascension in deg")
55parser.add_argument("-DEC", "--declination", default="78.46806", help="declination in deg")
57args = parser.parse_args()
58config = vars(args)
60print(config)
62location_data = (
63 config["directory"] + config["date"] + "/" + config["run_id"] + "/" + config["add_string"] + "/DL3"
64) # path to DL3 folder
65source_name = config["run_id"] # e.g., Crab, GRB210807A
66cut_type = "standard" # e.g., loose, hard, ...
67filename_output = "{name}_{cut}".format(name=source_name, cut=cut_type)
69source_position = SkyCoord(ra=config["right_ascension"], dec=config["declination"], unit="deg", frame="icrs")
70max_offset_run = 5 * u.deg
71work_directory = location_data
72path_plot = Path(work_directory + "/../plots")
73print(work_directory + "/../plots")
74path_plot.mkdir(exist_ok=True)
77e_min = 0.05 * u.TeV
78e_max = 20.0 * u.TeV
79on_radius = 0.2 * u.deg
80exclusion_radius = 0.2 * u.deg
81max_theta_squared = (3 * exclusion_radius.to(u.deg).value) ** 2
82n_bin_per_square_deg = 120
84data_store = DataStore.from_dir(location_data)
86print(data_store.obs_table)
88obs_ids = data_store.obs_table[source_position.separation(data_store.obs_table.pointing_radec) < max_offset_run][
89 "OBS_ID"
90]
91print(len(obs_ids))
92obs_collection = data_store.get_observations(obs_ids, required_irf="point-like")
94print("OBS_ID ", obs_ids)
96for i in range(len(obs_collection)):
97 obs_collection[i].events.peek()
98 plt.savefig(os.path.join(path_plot, "{}__obs.png".format(filename_output)), dpi=300)
99 plt.close()
101on_region = CircleSkyRegion(center=source_position, radius=on_radius)
102exclude_region = CircleSkyRegion(center=source_position, radius=exclusion_radius)
104energy_axis = MapAxis.from_edges(
105 np.logspace(np.log10(e_min.to_value(u.TeV)), np.log10(e_max.to_value(u.TeV)), 2),
106 unit="TeV",
107 name="energy",
108 interp="log",
109)
110geom_on_region = RegionGeom.create(region=on_region, axes=[energy_axis])
111geom = WcsGeom.create(skydir=source_position, npix=(200, 200), binsz=0.05, frame="icrs", axes=[energy_axis])
113geom_image = geom.to_image()
114exclusion_mask = ~geom_image.region_mask([exclude_region])
115# exclusion_mask.sum_over_axes().plot();
117stacked = MapDataset.create(geom=geom, name=source_name + "_stacked")
118unstacked = Datasets()
119maker = MapDatasetMaker(selection=["counts"])
120maker_safe_mask = SafeMaskMaker(methods=["offset-max"], offset_max="3 deg")
122for obs in obs_collection:
123 cutout = stacked.cutout(obs.pointing_radec, width="5 deg")
124 # cutout = stacked.cutout(obs.get_pointing_icrs(), width="5 deg")
125 dataset = maker.run(cutout, obs)
126 dataset = maker_safe_mask.run(dataset, obs)
127 stacked.stack(dataset)
128 unstacked.append(dataset)
130stacked.to_image().counts.smooth("0.12 deg", kernel="disk").plot(add_cbar=True, cmap="coolwarm")
131cbar = plt.gca().images[-1].colorbar
132cbar.ax.get_yaxis().labelpad = 5
133cbar.ax.set_ylabel(r"events/0.0004 deg$^2$", rotation=90)
134plt.gca().scatter(
135 source_position.ra,
136 source_position.dec,
137 transform=plt.gca().get_transform("world"),
138 marker="+",
139 c="k",
140 label=source_name,
141 s=100,
142)
143plt.legend()
145plt.savefig(os.path.join(path_plot, "{}__countmap.png".format(filename_output)), dpi=300)
146plt.close()
148dataset_maker_spectrum_significance = SpectrumDatasetMaker(selection=["counts", "exposure"], use_region_center=True)
149spectrum_dataset_empty_significance = SpectrumDataset.create(geom=geom_on_region)
150bkg_maker_spectrum_significance = ReflectedRegionsBackgroundMaker(exclusion_mask=exclusion_mask)
152datasets_significance = Datasets()
153for obs in obs_collection:
154 dataset_spectrum_significance = dataset_maker_spectrum_significance.run(
155 spectrum_dataset_empty_significance.copy(name=f"obs-{obs.obs_id}"), obs
156 )
157 dataset_on_off_spectrum_significance = bkg_maker_spectrum_significance.run(
158 observation=obs, dataset=dataset_spectrum_significance
159 )
160 datasets_significance.append(dataset_on_off_spectrum_significance)
162info_table = datasets_significance.info_table(cumulative=True)
165print("Livetime :", info_table["livetime"].to("s")[-1])
166print("Significance :", info_table["sqrt_ts"][-1])
167print("Excess :", info_table["excess"][-1])
168print("S/B :", info_table["excess"][-1] / info_table["background"][-1])
169print("Rate :", info_table["excess_rate"].to(1.0 / u.min)[-1].value, "per min")
171theta2_axis = MapAxis.from_bounds(
172 0, max_theta_squared, nbin=int(n_bin_per_square_deg * max_theta_squared), interp="lin", unit="deg2"
173)
175theta2_table = make_theta_squared_table(
176 observations=obs_collection,
177 position=source_position,
178 theta_squared_axis=theta2_axis,
179)
181# Values to be used for the plots
182counts_on = theta2_table["counts"]
183counts_on_err = np.sqrt(counts_on)
185counts_off = theta2_table["counts_off"]
186counts_off_err = np.sqrt(counts_off)
188excess = theta2_table["excess"]
189excess_err = np.append(
190 np.abs(theta2_table["excess_errn"].data[np.newaxis, :]), theta2_table["excess_errp"].data[np.newaxis, :], axis=0
191)
193bin_center = theta2_axis.center
194theta2_cut = on_radius.to_value(u.deg) ** 2
196N_off_err = np.sqrt(info_table["counts_off"][-1])
197significance_lima = info_table["sqrt_ts"][-1]
198text = (
199 r"Livetime = {:.1f}"
200 f"\n"
201 r"N$_{{\rm on}}$ = {:.0f}"
202 f"\n"
203 r"1/$\alpha$ = {:.1f}"
204 f"\n"
205 r"N$_{{\rm off}}$ = {:.1f} $\pm$ {:.1f}"
206 f"\n"
207 r"N_excess = {:.1f}"
208 f"\n"
209 r"Significance (Li&Ma) = {:.1f} $\sigma$ ".format(
210 info_table["livetime"].to("s")[-1],
211 info_table["counts"][-1],
212 1 / info_table["alpha"][-1],
213 info_table["counts_off"][-1],
214 N_off_err,
215 info_table["excess"][-1],
216 significance_lima,
217 )
218)
220fig, ax = plt.subplots(figsize=(8, 6))
222ax.errorbar(
223 bin_center,
224 counts_on,
225 xerr=theta2_axis.bin_width / 2,
226 yerr=counts_on_err,
227 fmt="o",
228 label="ON data",
229 ms=10,
230 color="tab:orange",
231)
232ax.errorbar(
233 bin_center,
234 counts_off,
235 xerr=theta2_axis.bin_width / 2,
236 yerr=counts_off_err,
237 fmt="s",
238 label="OFF data",
239 ms=10,
240 color="tab:blue",
241)
242ax.set_xlim(left=0)
243ax.grid(ls="-")
244ax.axvline(theta2_cut, color="black", ls="--", alpha=0.75)
245ax.set_xlabel("$\\theta^2$ [deg$^2$]")
246ax.set_ylabel("Counts")
247ax.legend(loc="center right")
248ax.set_title(source_name)
250box_prop = dict(boxstyle="Round", facecolor="white", alpha=0.5)
251text_prop = dict(fontsize="x-large", bbox=box_prop)
252txt = AnchoredText(text, loc=1, transform=ax.transAxes, prop=text_prop, frameon=False)
253ax.add_artist(txt)
255plt.savefig(os.path.join(path_plot, "{}_theta2.png".format(filename_output)), dpi=300)
256plt.close()
258# Excess plot
259fig, ax = plt.subplots(figsize=(12, 10))
261ax.errorbar(bin_center, excess, yerr=excess_err, fmt="o", color="tab:blue", label="Excess")
262ax.bar(bin_center, excess, theta2_axis.bin_width.to_value(u.deg**2), color="lightblue", ec="lightblue", alpha=0.5)
263# ax.axhline(0)
264ax.set_xlim(left=0)
265ax.axvline(theta2_cut, color="black", ls="--", alpha=0.75)
266ax.grid(ls="-")
267ax.set_xlabel("$\\theta^2$ [deg$^2$]")
268ax.set_ylabel("Excess")
269ax.legend(loc="upper right", title=f"Significance = {significance_lima:.1f} $\sigma$")
271plt.savefig(os.path.join(path_plot, "{}_theta2_excess.png".format(filename_output)), dpi=300)
272plt.close()
274excess_err_livetime = np.sqrt(
275 np.sqrt(info_table["counts"]) ** 2 + (np.sqrt(info_table["counts_off"]) * info_table["alpha"]) ** 2
276)
278excess_livetime_function = lambda x, a: a * x
279significance_livetime_function = lambda x, a: a * np.sqrt(x)
281popt_excess_livetime, pcov_excess_livetime = curve_fit(
282 excess_livetime_function, info_table["livetime"].to(u.h), info_table["excess"], sigma=excess_err_livetime
283)
284popt_significance_livetime, pcov_significance_livetime = curve_fit(
285 significance_livetime_function, info_table["livetime"].to(u.h), info_table["sqrt_ts"]
286)
288popt_excess_livetime = popt_excess_livetime[0] * u.Unit("h^-1")
289pcov_excess_livetime = pcov_excess_livetime[0, 0] * u.Unit("h^-2")
290popt_significance_livetime = popt_significance_livetime[0] * u.Unit("h^(-1/2)")
291pcov_significance_livetime = pcov_significance_livetime[0, 0] * u.Unit("h^(-1)")
293x_plot_livetime = np.linspace(0, np.max(info_table["livetime"].to(u.h)) * 1.2, num=len(info_table["livetime"]) * 10)
295y_plot_fit_excess_livetime = excess_livetime_function(x_plot_livetime, popt_excess_livetime)
296y_err_plot_fit_excess_livetime = x_plot_livetime * np.sqrt(pcov_excess_livetime)
298y_plot_fit_significance_livetime = significance_livetime_function(x_plot_livetime, popt_significance_livetime)
299y_err_plot_fit_significance_livetime = np.sqrt(x_plot_livetime) * np.sqrt(pcov_significance_livetime)
301fig, ax = plt.subplots(figsize=(15, 7), ncols=2)
302ax1, ax2 = ax[0], ax[1]
303ax1.fill_between(
304 x_plot_livetime.to_value(u.h),
305 y_plot_fit_excess_livetime - y_err_plot_fit_excess_livetime,
306 y_plot_fit_excess_livetime + y_err_plot_fit_excess_livetime,
307 color="silver",
308 alpha=0.5,
309)
310ax1.plot(x_plot_livetime, y_plot_fit_excess_livetime, c="silver")
311ax1.errorbar(
312 info_table["livetime"].to(u.h),
313 info_table["excess"],
314 yerr=excess_err_livetime,
315 fmt="o",
316 color="tab:blue",
317 label="Excess",
318)
319ax1.set_xlabel("Livetime [h]")
320ax1.set_ylabel("Excess")
321txt_ax1 = AnchoredText(
322 "Excess rate : {:.1f} +/- {:.1f} ".format(popt_excess_livetime.value, np.sqrt(pcov_excess_livetime).value)
323 + str(popt_excess_livetime.unit),
324 loc=2,
325 transform=ax1.transAxes,
326 prop=text_prop,
327 frameon=False,
328)
329ax1.add_artist(txt_ax1)
331ax2.fill_between(
332 x_plot_livetime.to_value(u.h),
333 y_plot_fit_significance_livetime - y_err_plot_fit_significance_livetime,
334 y_plot_fit_significance_livetime + y_err_plot_fit_significance_livetime,
335 color="silver",
336 alpha=0.5,
337)
338ax2.plot(x_plot_livetime, y_plot_fit_significance_livetime, c="silver")
339ax2.errorbar(info_table["livetime"].to(u.h), info_table["sqrt_ts"], fmt="o", color="tab:blue", label="Excess")
340txt_ax2 = AnchoredText(
341 "Significance rate : {:.1f} +/- {:.1f} ".format(
342 popt_significance_livetime.value, np.sqrt(pcov_significance_livetime).value
343 )
344 + str(popt_significance_livetime.unit),
345 loc=2,
346 transform=ax2.transAxes,
347 prop=text_prop,
348 frameon=False,
349)
350ax2.add_artist(txt_ax2)
351ax2.set_xlabel("Livetime [h]")
352ax2.set_ylabel("Significance")
354plt.savefig(os.path.join(path_plot, "{}__excess_significance_over_time.png".format(filename_output)), dpi=300)
355plt.close()
357ax = stacked.to_image().counts.smooth("0.12 deg", kernel="disk").plot(add_cbar=True, cmap="coolwarm")
358on_region.to_pixel(ax.wcs).plot(ax=ax, color="black")
359plot_spectrum_datasets_off_regions(datasets=datasets_significance, ax=ax)
360plt.gca().scatter(
361 source_position.ra,
362 source_position.dec,
363 transform=plt.gca().get_transform("world"),
364 marker="+",
365 c="k",
366 label=source_name,
367 s=[100],
368)
369plt.savefig(os.path.join(path_plot, "{}__countmap_with_OFF.png".format(filename_output)), dpi=300)
370plt.close()