Coverage for lst_auto_rta/Theta_square.py: 0%
181 statements
« prev ^ index » next coverage.py v7.6.10, created at 2025-01-26 14:48 +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
78#e_min = 0.5 * u.TeV
79e_max = 20.0 * u.TeV
80on_radius = 0.2 * u.deg
81exclusion_radius = 0.2 * u.deg
82max_theta_squared = (3 * exclusion_radius.to(u.deg).value) ** 2
83n_bin_per_square_deg = 120
85data_store = DataStore.from_dir(location_data)
87print(data_store.obs_table)
89obs_ids = data_store.obs_table[source_position.separation(data_store.obs_table.pointing_radec) < max_offset_run][
90 "OBS_ID"
91]
92print(len(obs_ids))
93obs_collection = data_store.get_observations(obs_ids, required_irf="point-like")
95print("OBS_ID ", obs_ids)
97for i in range(len(obs_collection)):
98 obs_collection[i].events.peek()
99 plt.savefig(os.path.join(path_plot, "{}__obs.png".format(filename_output)), dpi=300)
100 plt.close()
102on_region = CircleSkyRegion(center=source_position, radius=on_radius)
103exclude_region = CircleSkyRegion(center=source_position, radius=exclusion_radius)
105energy_axis = MapAxis.from_edges(
106 np.logspace(np.log10(e_min.to_value(u.TeV)), np.log10(e_max.to_value(u.TeV)), 2),
107 unit="TeV",
108 name="energy",
109 interp="log",
110)
111geom_on_region = RegionGeom.create(region=on_region, axes=[energy_axis])
112geom = WcsGeom.create(skydir=source_position, npix=(200, 200), binsz=0.05, frame="icrs", axes=[energy_axis])
114geom_image = geom.to_image()
115exclusion_mask = ~geom_image.region_mask([exclude_region])
116# exclusion_mask.sum_over_axes().plot();
118stacked = MapDataset.create(geom=geom, name=source_name + "_stacked")
119unstacked = Datasets()
120maker = MapDatasetMaker(selection=["counts"])
121maker_safe_mask = SafeMaskMaker(methods=["offset-max"], offset_max="3 deg")
123for obs in obs_collection:
124 cutout = stacked.cutout(obs.pointing_radec, width="5 deg")
125 # cutout = stacked.cutout(obs.get_pointing_icrs(), width="5 deg")
126 dataset = maker.run(cutout, obs)
127 dataset = maker_safe_mask.run(dataset, obs)
128 stacked.stack(dataset)
129 unstacked.append(dataset)
131stacked.to_image().counts.smooth("0.12 deg", kernel="disk").plot(add_cbar=True, cmap="coolwarm")
132cbar = plt.gca().images[-1].colorbar
133cbar.ax.get_yaxis().labelpad = 5
134cbar.ax.set_ylabel(r"events/0.0004 deg$^2$", rotation=90)
135plt.gca().scatter(
136 source_position.ra,
137 source_position.dec,
138 transform=plt.gca().get_transform("world"),
139 marker="+",
140 c="k",
141 label=source_name,
142 s=100,
143)
144plt.legend()
146plt.savefig(os.path.join(path_plot, "{}__countmap.png".format(filename_output)), dpi=300)
147plt.close()
149dataset_maker_spectrum_significance = SpectrumDatasetMaker(selection=["counts", "exposure"], use_region_center=True)
150spectrum_dataset_empty_significance = SpectrumDataset.create(geom=geom_on_region)
151bkg_maker_spectrum_significance = ReflectedRegionsBackgroundMaker(exclusion_mask=exclusion_mask)
153datasets_significance = Datasets()
154for obs in obs_collection:
155 dataset_spectrum_significance = dataset_maker_spectrum_significance.run(
156 spectrum_dataset_empty_significance.copy(name=f"obs-{obs.obs_id}"), obs
157 )
158 dataset_on_off_spectrum_significance = bkg_maker_spectrum_significance.run(
159 observation=obs, dataset=dataset_spectrum_significance
160 )
161 datasets_significance.append(dataset_on_off_spectrum_significance)
163info_table = datasets_significance.info_table(cumulative=True)
166print("Livetime :", info_table["livetime"].to("s")[-1])
167print("Significance :", info_table["sqrt_ts"][-1])
168print("Excess :", info_table["excess"][-1])
169print("S/B :", info_table["excess"][-1] / info_table["background"][-1])
170print("Rate :", info_table["excess_rate"].to(1.0 / u.min)[-1].value, "per min")
172theta2_axis = MapAxis.from_bounds(
173 0, max_theta_squared, nbin=int(n_bin_per_square_deg * max_theta_squared), interp="lin", unit="deg2"
174)
176theta2_table = make_theta_squared_table(
177 observations=obs_collection,
178 position=source_position,
179 theta_squared_axis=theta2_axis,
180)
182# Values to be used for the plots
183counts_on = theta2_table["counts"]
184counts_on_err = np.sqrt(counts_on)
186counts_off = theta2_table["counts_off"]
187counts_off_err = np.sqrt(counts_off)
189excess = theta2_table["excess"]
190excess_err = np.append(
191 np.abs(theta2_table["excess_errn"].data[np.newaxis, :]), theta2_table["excess_errp"].data[np.newaxis, :], axis=0
192)
194bin_center = theta2_axis.center
195theta2_cut = on_radius.to_value(u.deg) ** 2
197N_off_err = np.sqrt(info_table["counts_off"][-1])
198significance_lima = info_table["sqrt_ts"][-1]
199text = (
200 r"Livetime = {:.1f}"
201 f"\n"
202 r"N$_{{\rm on}}$ = {:.0f}"
203 f"\n"
204 r"1/$\alpha$ = {:.1f}"
205 f"\n"
206 r"N$_{{\rm off}}$ = {:.1f} $\pm$ {:.1f}"
207 f"\n"
208 r"N_excess = {:.1f}"
209 f"\n"
210 r"Significance (Li&Ma) = {:.1f} $\sigma$ ".format(
211 info_table["livetime"].to("s")[-1],
212 info_table["counts"][-1],
213 1 / info_table["alpha"][-1],
214 info_table["counts_off"][-1],
215 N_off_err,
216 info_table["excess"][-1],
217 significance_lima,
218 )
219)
221fig, ax = plt.subplots(figsize=(8, 6))
223ax.errorbar(
224 bin_center,
225 counts_on,
226 xerr=theta2_axis.bin_width / 2,
227 yerr=counts_on_err,
228 fmt="o",
229 label="ON data",
230 ms=10,
231 color="tab:orange",
232)
233ax.errorbar(
234 bin_center,
235 counts_off,
236 xerr=theta2_axis.bin_width / 2,
237 yerr=counts_off_err,
238 fmt="s",
239 label="OFF data",
240 ms=10,
241 color="tab:blue",
242)
243ax.set_xlim(left=0)
244ax.grid(ls="-")
245ax.axvline(theta2_cut, color="black", ls="--", alpha=0.75)
246ax.set_xlabel("$\\theta^2$ [deg$^2$]")
247ax.set_ylabel("Counts")
248ax.legend(loc="center right")
249ax.set_title(source_name)
251box_prop = dict(boxstyle="Round", facecolor="white", alpha=0.5)
252text_prop = dict(fontsize="x-large", bbox=box_prop)
253txt = AnchoredText(text, loc=1, transform=ax.transAxes, prop=text_prop, frameon=False)
254ax.add_artist(txt)
256plt.savefig(os.path.join(path_plot, "{}_theta2.png".format(filename_output)), dpi=300)
257plt.close()
259# Excess plot
260fig, ax = plt.subplots(figsize=(12, 10))
262ax.errorbar(bin_center, excess, yerr=excess_err, fmt="o", color="tab:blue", label="Excess")
263ax.bar(bin_center, excess, theta2_axis.bin_width.to_value(u.deg**2), color="lightblue", ec="lightblue", alpha=0.5)
264# ax.axhline(0)
265ax.set_xlim(left=0)
266ax.axvline(theta2_cut, color="black", ls="--", alpha=0.75)
267ax.grid(ls="-")
268ax.set_xlabel("$\\theta^2$ [deg$^2$]")
269ax.set_ylabel("Excess")
270ax.legend(loc="upper right", title=f"Significance = {significance_lima:.1f} $\sigma$")
272plt.savefig(os.path.join(path_plot, "{}_theta2_excess.png".format(filename_output)), dpi=300)
273plt.close()
275excess_err_livetime = np.sqrt(
276 np.sqrt(info_table["counts"]) ** 2 + (np.sqrt(info_table["counts_off"]) * info_table["alpha"]) ** 2
277)
279excess_livetime_function = lambda x, a: a * x
280significance_livetime_function = lambda x, a: a * np.sqrt(x)
282popt_excess_livetime, pcov_excess_livetime = curve_fit(
283 excess_livetime_function, info_table["livetime"].to(u.h), info_table["excess"], sigma=excess_err_livetime
284)
285popt_significance_livetime, pcov_significance_livetime = curve_fit(
286 significance_livetime_function, info_table["livetime"].to(u.h), info_table["sqrt_ts"]
287)
289popt_excess_livetime = popt_excess_livetime[0] * u.Unit("h^-1")
290pcov_excess_livetime = pcov_excess_livetime[0, 0] * u.Unit("h^-2")
291popt_significance_livetime = popt_significance_livetime[0] * u.Unit("h^(-1/2)")
292pcov_significance_livetime = pcov_significance_livetime[0, 0] * u.Unit("h^(-1)")
294x_plot_livetime = np.linspace(0, np.max(info_table["livetime"].to(u.h)) * 1.2, num=len(info_table["livetime"]) * 10)
296y_plot_fit_excess_livetime = excess_livetime_function(x_plot_livetime, popt_excess_livetime)
297y_err_plot_fit_excess_livetime = x_plot_livetime * np.sqrt(pcov_excess_livetime)
299y_plot_fit_significance_livetime = significance_livetime_function(x_plot_livetime, popt_significance_livetime)
300y_err_plot_fit_significance_livetime = np.sqrt(x_plot_livetime) * np.sqrt(pcov_significance_livetime)
302fig, ax = plt.subplots(figsize=(15, 7), ncols=2)
303ax1, ax2 = ax[0], ax[1]
304ax1.fill_between(
305 x_plot_livetime.to_value(u.h),
306 y_plot_fit_excess_livetime - y_err_plot_fit_excess_livetime,
307 y_plot_fit_excess_livetime + y_err_plot_fit_excess_livetime,
308 color="silver",
309 alpha=0.5,
310)
311ax1.plot(x_plot_livetime, y_plot_fit_excess_livetime, c="silver")
312ax1.errorbar(
313 info_table["livetime"].to(u.h),
314 info_table["excess"],
315 yerr=excess_err_livetime,
316 fmt="o",
317 color="tab:blue",
318 label="Excess",
319)
320ax1.set_xlabel("Livetime [h]")
321ax1.set_ylabel("Excess")
322txt_ax1 = AnchoredText(
323 "Excess rate : {:.1f} +/- {:.1f} ".format(popt_excess_livetime.value, np.sqrt(pcov_excess_livetime).value)
324 + str(popt_excess_livetime.unit),
325 loc=2,
326 transform=ax1.transAxes,
327 prop=text_prop,
328 frameon=False,
329)
330ax1.add_artist(txt_ax1)
332ax2.fill_between(
333 x_plot_livetime.to_value(u.h),
334 y_plot_fit_significance_livetime - y_err_plot_fit_significance_livetime,
335 y_plot_fit_significance_livetime + y_err_plot_fit_significance_livetime,
336 color="silver",
337 alpha=0.5,
338)
339ax2.plot(x_plot_livetime, y_plot_fit_significance_livetime, c="silver")
340ax2.errorbar(info_table["livetime"].to(u.h), info_table["sqrt_ts"], fmt="o", color="tab:blue", label="Excess")
341txt_ax2 = AnchoredText(
342 "Significance rate : {:.1f} +/- {:.1f} ".format(
343 popt_significance_livetime.value, np.sqrt(pcov_significance_livetime).value
344 )
345 + str(popt_significance_livetime.unit),
346 loc=2,
347 transform=ax2.transAxes,
348 prop=text_prop,
349 frameon=False,
350)
351ax2.add_artist(txt_ax2)
352ax2.set_xlabel("Livetime [h]")
353ax2.set_ylabel("Significance")
355plt.savefig(os.path.join(path_plot, "{}__excess_significance_over_time.png".format(filename_output)), dpi=300)
356plt.close()
358ax = stacked.to_image().counts.smooth("0.12 deg", kernel="disk").plot(add_cbar=True, cmap="coolwarm")
359on_region.to_pixel(ax.wcs).plot(ax=ax, color="black")
360plot_spectrum_datasets_off_regions(datasets=datasets_significance, ax=ax)
361plt.gca().scatter(
362 source_position.ra,
363 source_position.dec,
364 transform=plt.gca().get_transform("world"),
365 marker="+",
366 c="k",
367 label=source_name,
368 s=[100],
369)
370plt.savefig(os.path.join(path_plot, "{}__countmap_with_OFF.png".format(filename_output)), dpi=300)
371plt.close()