Coverage for src / lstautorta / Ring_Background_Maps_ACADA.py: 0%
144 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-08-10 11:56 +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 numpy as np
20from astropy.coordinates import SkyCoord
21from matplotlib import style
23style.use("tableau-colorblind10")
24import argparse
25from pathlib import Path
27from acceptance_modelisation import RadialAcceptanceMapCreator
28from gammapy.data import DataStore
29from gammapy.datasets import (
30 Datasets,
31 MapDataset,
32)
33from gammapy.estimators import ExcessMapEstimator
34from gammapy.makers import (
35 MapDatasetMaker,
36 RingBackgroundMaker,
37 SafeMaskMaker,
38)
39from gammapy.maps import MapAxis, WcsGeom
40from matplotlib.offsetbox import AnchoredText
41from regions import CircleSkyRegion
42from scipy.stats import norm
44parser = argparse.ArgumentParser(
45 description="Automatic Script for the DL1 check", formatter_class=argparse.ArgumentDefaultsHelpFormatter
46)
47parser.add_argument("-d", "--directory", default="/fefs/onsite/pipeline/rta/data/", help="Directory for data")
48parser.add_argument("-da", "--date", default="20230705", help="Date of the run to check")
49parser.add_argument("-r", "--run-id", default="13600", help="run id to check")
50parser.add_argument("-add", "--add-string", default="reco/", help="add a string to the path")
51parser.add_argument("-RA", "--right-ascension", default="270.19042", help="right-ascension in deg")
52parser.add_argument("-DEC", "--declination", default="78.46806", help="declination in deg")
54args = parser.parse_args()
55config = vars(args)
57location_data = (
58 config["directory"] + config["date"] + "/" + config["run_id"] + "/" + config["add_string"] + "/plots"
59) # path to DL3 folder
60source_name = config["run_id"] # e.g., Crab, GRB210807A
61cut_type = "standard" # e.g., loose, hard, ...
62filename_output = f"{source_name}_{cut_type}"
64source_position = SkyCoord(ra=config["right_ascension"], dec=config["declination"], unit="deg", frame="icrs")
65max_offset_run = 5 * u.deg
66work_directory = location_data
67path_plot = Path(work_directory + "/")
68print(work_directory + "/")
69path_plot.mkdir(exist_ok=True)
70path_background = Path(work_directory + "/")
71path_background.mkdir(exist_ok=True)
73e_min = 0.05 * u.TeV
74e_max = 10.0 * u.TeV
75n_bin_per_decade = 10
76on_radius = 0.2 * u.deg
77exclusion_radius = 0.35 * u.deg
78fov_observation = 4.5 * u.deg
80r_in = 0.5 * u.deg
81width = 0.4 * u.deg
82correlation_radius = 0.2 * u.deg
84n_bin_per_decade_acceptance = 2.5
85offset_bin_size_acceptance = 0.4 * u.deg
87data_store = DataStore.from_dir(location_data)
88data_store.info()
90obs_ids = data_store.obs_table[source_position.separation(data_store.obs_table.pointing_radec) < max_offset_run][
91 "OBS_ID"
92]
93obs_collection = data_store.get_observations(obs_ids, required_irf=None)
95exclude_region = CircleSkyRegion(center=source_position, radius=exclusion_radius)
97n_bin_energy = int((np.log10(e_max.to_value(u.TeV)) - np.log10(e_min.to_value(u.TeV))) * n_bin_per_decade)
98energy_axis = MapAxis.from_edges(
99 np.logspace(np.log10(e_min.to_value(u.TeV)), np.log10(e_max.to_value(u.TeV)), n_bin_energy + 1),
100 unit="TeV",
101 name="energy",
102 interp="log",
103)
104maximal_run_separation = np.max(
105 source_position.separation(data_store.obs_table[np.isin(data_store.obs_table["OBS_ID"], obs_ids)].pointing_radec)
106)
107geom = WcsGeom.create(
108 skydir=source_position,
109 width=((maximal_run_separation + fov_observation) * 1.5, (maximal_run_separation + fov_observation) * 1.5),
110 binsz=0.02,
111 frame="icrs",
112 axes=[energy_axis],
113)
115geom_image = geom.to_image()
116exclusion_mask = ~geom_image.region_mask([exclude_region])
118stacked = MapDataset.create(geom=geom, name=source_name + "_stacked")
119unstacked = Datasets()
120maker = MapDatasetMaker(selection=["counts"])
121maker_safe_mask = SafeMaskMaker(methods=["offset-max"], offset_max=fov_observation)
123for obs in obs_collection:
124 cutout = stacked.cutout(obs.pointing_radec, width="6.5 deg")
125 dataset = maker.run(cutout, obs)
126 dataset = maker_safe_mask.run(dataset, obs)
127 stacked.stack(dataset)
128 unstacked.append(dataset)
130n_bin_energy_acceptance = int(
131 (np.log10(e_max.to_value(u.TeV)) - np.log10(e_min.to_value(u.TeV))) * n_bin_per_decade_acceptance
132)
133energyAxisAcceptance = MapAxis.from_edges(
134 np.logspace(np.log10(e_min.to_value(u.TeV)), np.log10(e_max.to_value(u.TeV)), 1 + n_bin_energy_acceptance),
135 unit="TeV",
136 name="energy",
137 interp="log",
138)
139n_bin_offset_acceptance = int(fov_observation.to_value(u.deg) / offset_bin_size_acceptance.to_value(u.deg))
140offsetAxisAcceptance = MapAxis.from_edges(
141 np.linspace(0.0, fov_observation.to_value(u.deg), 1 + n_bin_offset_acceptance),
142 unit="deg",
143 name="offset",
144 interp="lin",
145)
147background_creator = RadialAcceptanceMapCreator(
148 energyAxisAcceptance,
149 offsetAxisAcceptance,
150 exclude_regions=[
151 exclude_region,
152 ],
153 oversample_map=10,
154)
155background = background_creator.create_radial_acceptance_map_per_observation(obs_collection)
156# background[list(background.keys())[0]].peek()
157for obs_id in background.keys():
158 hdu_background = background[obs_id].to_table_hdu()
159 hdu_background.writeto(
160 os.path.join(path_background, filename_output + "_" + str(obs_id).zfill(5) + "_background.fits"),
161 overwrite=True,
162 )
164data_store.hdu_table.remove_rows(data_store.hdu_table["HDU_TYPE"] == "bkg")
166for obs_id in np.unique(data_store.hdu_table["OBS_ID"]):
167 data_store.hdu_table.add_row(
168 {
169 "OBS_ID": obs_id,
170 "HDU_TYPE": "bkg",
171 "HDU_CLASS": "bkg_2d",
172 "FILE_DIR": "",
173 "FILE_NAME": os.path.join(path_background, filename_output + "_" + str(obs_id) + "_background.fits"),
174 "HDU_NAME": "BACKGROUND",
175 "SIZE": hdu_background.size,
176 }
177 )
179data_store.hdu_table = data_store.hdu_table.copy()
180obs_collection = data_store.get_observations(obs_ids, required_irf=None)
182stacked = MapDataset.create(geom=geom)
183unstacked = Datasets()
184maker = MapDatasetMaker(selection=["counts", "background"])
185maker_safe_mask = SafeMaskMaker(methods=["offset-max"], offset_max=fov_observation)
187for obs in obs_collection:
188 cutout = stacked.cutout(obs.pointing_radec, width="6.5 deg")
189 dataset = maker.run(cutout, obs)
190 dataset = maker_safe_mask.run(dataset, obs)
191 stacked.stack(dataset)
192 unstacked.append(dataset)
194ring_bkg_maker = RingBackgroundMaker(r_in=r_in, width=width) # , exclusion_mask=exclusion_mask)
195stacked_ring = ring_bkg_maker.run(stacked.to_image())
196estimator = ExcessMapEstimator(correlation_radius, correlate_off=False)
197lima_maps = estimator.run(stacked_ring)
199significance_all = lima_maps["sqrt_ts"].data[np.isfinite(lima_maps["sqrt_ts"].data)]
200significance_background = lima_maps["sqrt_ts"].data[
201 np.logical_and(np.isfinite(lima_maps["sqrt_ts"].data), exclusion_mask.data)
202]
204bins = np.linspace(
205 np.min(significance_all),
206 np.max(significance_all),
207 num=int((np.max(significance_all) - np.min(significance_all)) * 3),
208)
210# Now, fit the off distribution with a Gaussian
211mu, std = norm.fit(significance_background)
212x = np.linspace(-8, 8, 50)
213p = norm.pdf(x, mu, std)
215plt.figure(figsize=(8, 21))
216ax1 = plt.subplot(3, 1, 1, projection=lima_maps["sqrt_ts"].geom.wcs)
217ax2 = plt.subplot(3, 1, 2, projection=lima_maps["sqrt_ts"].geom.wcs)
218ax3 = plt.subplot(3, 1, 3)
220ax2.set_title("Significance map")
221lima_maps["sqrt_ts"].plot(ax=ax2, add_cbar=True)
222ax2.scatter(
223 source_position.ra,
224 source_position.dec,
225 transform=ax2.get_transform("world"),
226 marker="+",
227 c="red",
228 label=filename_output,
229 s=[300],
230 linewidths=3,
231)
232ax2.legend()
234ax1.set_title("Excess map")
235lima_maps["npred_excess"].plot(ax=ax1, add_cbar=True)
236ax1.scatter(
237 source_position.ra,
238 source_position.dec,
239 transform=ax1.get_transform("world"),
240 marker="+",
241 c="red",
242 label=filename_output,
243 s=[300],
244 linewidths=3,
245)
246ax1.legend()
248ax3.set_title("Significance distribution")
249ax3.hist(significance_all, density=True, alpha=0.5, color="red", label="All bins", bins=bins)
250ax3.hist(significance_background, density=True, alpha=0.5, color="blue", label="Background bins", bins=bins)
252ax3.plot(x, p, lw=2, color="black")
253ax3.legend()
254ax3.set_xlabel("Significance")
255ax3.set_yscale("log")
256ax3.set_ylim(1e-5, 1)
257xmin, xmax = np.min(significance_all), np.max(significance_all)
258ax3.set_xlim(xmin, xmax)
260text = text = rf"$\mu$ = {mu:.2f}" "\n" rf"$\sigma$ = {std:.2f}"
261box_prop = dict(boxstyle="Round", facecolor="white", alpha=0.5)
262text_prop = dict(fontsize="x-large", bbox=box_prop)
263# txt = AnchoredText(text, loc=2, transform=ax3.transAxes, prop=text_prop, frameon=False)
264txt = AnchoredText(text, loc=2, prop=text_prop, frameon=False)
265ax3.add_artist(txt)
267plt.savefig(os.path.join(path_plot, f"{filename_output}__sky_map.png"), dpi=300)
268print(f"Fit results: mu = {mu:.2f}, std = {std:.2f}")