Coverage for src / lstautorta / Ring_Background_Maps.py: 0%
170 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.estimators.utils import find_peaks
35from gammapy.makers import (
36 MapDatasetMaker,
37 RingBackgroundMaker,
38 SafeMaskMaker,
39)
40from gammapy.maps import MapAxis, WcsGeom
41from matplotlib.offsetbox import AnchoredText
42from regions import CircleSkyRegion
43from scipy.stats import norm
45parser = argparse.ArgumentParser(
46 description="Automatic Script for the DL1 check", formatter_class=argparse.ArgumentDefaultsHelpFormatter
47)
48parser.add_argument("-d", "--directory", default="/fefs/onsite/pipeline/rta/data/", help="Directory for data")
49parser.add_argument("-da", "--date", default="20230705", help="Date of the run to check")
50parser.add_argument("-r", "--run-id", default="13600", help="run id to check")
51parser.add_argument("-add", "--add-string", default="", help="add a string to the path")
52parser.add_argument("-RA", "--right-ascension", default="270.19042", help="right-ascension in deg")
53parser.add_argument("-DEC", "--declination", default="78.46806", help="declination in deg")
55args = parser.parse_args()
56config = vars(args)
58location_data = (
59 config["directory"] + config["date"] + "/" + config["run_id"] + "/" + config["add_string"] + "/DL3"
60) # path to DL3 folder
61source_name = config["run_id"] # e.g., Crab, GRB210807A
62cut_type = "standard" # e.g., loose, hard, ...
63filename_output = f"{source_name}_{cut_type}"
65source_position = SkyCoord(ra=config["right_ascension"], dec=config["declination"], unit="deg", frame="icrs")
66max_offset_run = 10 * u.deg
67work_directory = location_data
68path_plot = Path(work_directory + "/../plots")
69print(work_directory + "/../plots")
70path_plot.mkdir(exist_ok=True)
71path_background = Path(work_directory + "/../plots")
72path_background.mkdir(exist_ok=True)
74e_min = 0.05 * u.TeV
75e_max = 10.0 * u.TeV
76n_bin_per_decade = 10
77on_radius = 0.2 * u.deg
78exclusion_radius = 0.0 * u.deg
79fov_observation = 7 * u.deg
81r_in = 0.5 * u.deg
82width = 0.4 * u.deg
83correlation_radius = 0.2 * u.deg
85n_bin_per_decade_acceptance = 1.0
86offset_bin_size_acceptance = 0.4 * u.deg
88data_store = DataStore.from_dir(location_data)
89data_store.info()
91obs_ids = data_store.obs_table[source_position.separation(data_store.obs_table.pointing_radec) < max_offset_run][
92 "OBS_ID"
93]
94obs_collection = data_store.get_observations(obs_ids, required_irf=None)
96exclude_region = CircleSkyRegion(center=source_position, radius=exclusion_radius)
98n_bin_energy = int((np.log10(e_max.to_value(u.TeV)) - np.log10(e_min.to_value(u.TeV))) * n_bin_per_decade)
99energy_axis = MapAxis.from_edges(
100 np.logspace(np.log10(e_min.to_value(u.TeV)), np.log10(e_max.to_value(u.TeV)), n_bin_energy + 1),
101 unit="TeV",
102 name="energy",
103 interp="log",
104)
105# maximal_run_separation = np.max(
106# source_position.separation(data_store.obs_table[np.isin(data_store.obs_table["OBS_ID"], obs_ids)].pointing_radec)
107# )
108# geom = WcsGeom.create(
109# skydir=source_position,
110# width=((maximal_run_separation + fov_observation)*1.2, (maximal_run_separation + fov_observation)*1.2),
111# binsz=0.05,
112# frame="icrs",
113# axes=[energy_axis],
114# )
116pointing_coords = SkyCoord([obs.pointing_radec for obs in obs_collection])
118x = np.cos(pointing_coords.ra.radian) * np.cos(pointing_coords.dec.radian)
119y = np.sin(pointing_coords.ra.radian) * np.cos(pointing_coords.dec.radian)
120z = np.sin(pointing_coords.dec.radian)
122x_mean = np.mean(x)
123y_mean = np.mean(y)
124z_mean = np.mean(z)
126r = np.sqrt(x_mean**2 + y_mean**2 + z_mean**2)
127x_mean /= r
128y_mean /= r
129z_mean /= r
131ra_mean = np.arctan2(y_mean, x_mean) * u.rad
132dec_mean = np.arcsin(z_mean) * u.rad
134center_coord = SkyCoord(ra=ra_mean, dec=dec_mean, frame="icrs")
135max_separation = center_coord.separation(pointing_coords).max()
136map_width = 2 * max_separation + fov_observation
138geom = WcsGeom.create(
139 skydir=center_coord,
140 width=(map_width, map_width),
141 binsz=0.05,
142 frame="icrs",
143 axes=[energy_axis],
144)
146geom_image = geom.to_image()
147exclusion_mask = ~geom_image.region_mask([exclude_region])
149stacked = MapDataset.create(geom=geom, name=source_name + "_stacked")
150unstacked = Datasets()
151maker = MapDatasetMaker(selection=["counts"])
152maker_safe_mask = SafeMaskMaker(methods=["offset-max"], offset_max=fov_observation)
154for obs in obs_collection:
155 cutout = stacked.cutout(obs.pointing_radec, width="6.5 deg")
156 dataset = maker.run(cutout, obs)
157 dataset = maker_safe_mask.run(dataset, obs)
158 stacked.stack(dataset)
159 unstacked.append(dataset)
161n_bin_energy_acceptance = int(
162 (np.log10(e_max.to_value(u.TeV)) - np.log10(e_min.to_value(u.TeV))) * n_bin_per_decade_acceptance
163)
164energyAxisAcceptance = MapAxis.from_edges(
165 np.logspace(np.log10(e_min.to_value(u.TeV)), np.log10(e_max.to_value(u.TeV)), 1 + n_bin_energy_acceptance),
166 unit="TeV",
167 name="energy",
168 interp="log",
169)
170n_bin_offset_acceptance = int(fov_observation.to_value(u.deg) / offset_bin_size_acceptance.to_value(u.deg))
171offsetAxisAcceptance = MapAxis.from_edges(
172 np.linspace(0.0, fov_observation.to_value(u.deg), 1 + n_bin_offset_acceptance),
173 unit="deg",
174 name="offset",
175 interp="lin",
176)
178background_creator = RadialAcceptanceMapCreator(
179 energyAxisAcceptance,
180 offsetAxisAcceptance,
181 exclude_regions=[
182 exclude_region,
183 ],
184 oversample_map=10,
185)
186background = background_creator.create_radial_acceptance_map_per_observation(obs_collection)
187# background[list(background.keys())[0]].peek()
188for obs_id in background.keys():
189 hdu_background = background[obs_id].to_table_hdu()
190 hdu_background.writeto(
191 os.path.join(path_background, filename_output + "_" + str(obs_id) + "_background.fits"), overwrite=True
192 )
194data_store.hdu_table.remove_rows(data_store.hdu_table["HDU_TYPE"] == "bkg")
196for obs_id in np.unique(data_store.hdu_table["OBS_ID"]):
197 data_store.hdu_table.add_row(
198 {
199 "OBS_ID": obs_id,
200 "HDU_TYPE": "bkg",
201 "HDU_CLASS": "bkg_2d",
202 "FILE_DIR": "",
203 "FILE_NAME": os.path.join(path_background, filename_output + "_" + str(obs_id) + "_background.fits"),
204 "HDU_NAME": "BACKGROUND",
205 "SIZE": hdu_background.size,
206 }
207 )
209data_store.hdu_table = data_store.hdu_table.copy()
210obs_collection = data_store.get_observations(obs_ids, required_irf=None)
212stacked = MapDataset.create(geom=geom)
213unstacked = Datasets()
214maker = MapDatasetMaker(selection=["counts", "background"])
215maker_safe_mask = SafeMaskMaker(methods=["offset-max"], offset_max=fov_observation)
217for obs in obs_collection:
218 cutout = stacked.cutout(obs.pointing_radec, width="6.5 deg")
219 dataset = maker.run(cutout, obs)
220 dataset = maker_safe_mask.run(dataset, obs)
221 stacked.stack(dataset)
222 unstacked.append(dataset)
224ring_bkg_maker = RingBackgroundMaker(r_in=r_in, width=width) # , exclusion_mask=exclusion_mask)
225stacked_ring = ring_bkg_maker.run(stacked.to_image())
226estimator = ExcessMapEstimator(correlation_radius, correlate_off=False)
227lima_maps = estimator.run(stacked_ring)
229significance_all = lima_maps["sqrt_ts"].data[np.isfinite(lima_maps["sqrt_ts"].data)]
230significance_background = lima_maps["sqrt_ts"].data[
231 np.logical_and(np.isfinite(lima_maps["sqrt_ts"].data), exclusion_mask.data)
232]
234bins = np.linspace(
235 np.min(significance_all),
236 np.max(significance_all),
237 num=int((np.max(significance_all) - np.min(significance_all)) * 3),
238)
240# Now, fit the off distribution with a Gaussian
241mu, std = norm.fit(significance_background)
242x = np.linspace(-8, 8, 50)
243p = norm.pdf(x, mu, std)
245plt.figure(figsize=(8, 21))
246ax1 = plt.subplot(3, 1, 1, projection=lima_maps["sqrt_ts"].geom.wcs)
247ax2 = plt.subplot(3, 1, 2, projection=lima_maps["sqrt_ts"].geom.wcs)
248ax3 = plt.subplot(3, 1, 3)
250ax2.set_title("Significance map")
251lima_maps["sqrt_ts"].plot(ax=ax2, add_cbar=True)
252ax2.scatter(
253 source_position.ra,
254 source_position.dec,
255 transform=ax2.get_transform("world"),
256 marker="+",
257 c="red",
258 label=filename_output,
259 s=[300],
260 linewidths=3,
261)
262ax2.legend()
264sources = find_peaks(
265 lima_maps["sqrt_ts"].get_image_by_idx((0,)),
266 threshold=5,
267 min_distance="0.2 deg",
268)
269# sources = find_peaks(
270# lima_maps["npred_excess"].get_image_by_idx((0,)),
271# threshold=150,
272# min_distance="0.2 deg",
273# )
274print(sources)
275# now = dt.datetime.now()
276# timestamp_str = now.strftime("%Y-%m-%d %H:%M:%S")
277# ax1.text(0.02, 0.98, timestamp_str, transform=ax1.transAxes,
278# fontsize=11, fontweight='bold', va='top', ha='left')
279# ax2.text(0.02, 0.98, timestamp_str, transform=ax2.transAxes,
280# fontsize=11, fontweight='bold', va='top', ha='left')
281if len(sources) > 0:
282 ax2.scatter(
283 sources["ra"],
284 sources["dec"],
285 # transform=plt.gca().get_transform("icrs"),
286 transform=ax2.get_transform("icrs"),
287 color="none",
288 edgecolor="red",
289 marker="o",
290 s=300,
291 lw=1.5,
292 )
294for obs in obs_collection:
295 ax2.scatter(
296 obs.pointing_radec.ra,
297 obs.pointing_radec.dec,
298 transform=ax2.get_transform("icrs"),
299 marker="x",
300 c="black",
301 s=100,
302 linewidths=1.5,
303 )
305# Dummy entry for legend
306ax2.scatter([], [], marker="x", c="black", label="Pointings")
309ax1.set_title("Excess map")
310lima_maps["npred_excess"].plot(ax=ax1, add_cbar=True)
311ax1.scatter(
312 source_position.ra,
313 source_position.dec,
314 transform=ax1.get_transform("icrs"),
315 marker="+",
316 c="red",
317 label=filename_output,
318 s=[300],
319 linewidths=3,
320)
322for obs in obs_collection:
323 ax1.scatter(
324 obs.pointing_radec.ra,
325 obs.pointing_radec.dec,
326 transform=ax1.get_transform("icrs"),
327 marker="x",
328 c="black",
329 s=100,
330 linewidths=1.5,
331 )
333ax1.scatter([], [], marker="x", c="black", label="Pointings")
336ax1.legend()
338ax3.set_title("Significance distribution")
339ax3.hist(significance_all, density=True, alpha=0.5, color="red", label="All bins", bins=bins)
340ax3.hist(significance_background, density=True, alpha=0.5, color="blue", label="Background bins", bins=bins)
342ax3.plot(x, p, lw=2, color="black")
343ax3.legend()
344ax3.set_xlabel("Significance")
345ax3.set_yscale("log")
346ax3.set_ylim(1e-5, 1)
347xmin, xmax = np.min(significance_all), np.max(significance_all)
348ax3.set_xlim(xmin, xmax)
350text = text = rf"$\mu$ = {mu:.2f}" "\n" rf"$\sigma$ = {std:.2f}"
351box_prop = dict(boxstyle="Round", facecolor="white", alpha=0.5)
352text_prop = dict(fontsize="x-large", bbox=box_prop)
353# txt = AnchoredText(text, loc=2, transform=ax3.transAxes, prop=text_prop, frameon=False)
354txt = AnchoredText(text, loc=2, prop=text_prop, frameon=False)
355ax3.add_artist(txt)
357plt.savefig(os.path.join(path_plot, f"{filename_output}__sky_map.png"), dpi=300)
358print(f"Fit results: mu = {mu:.2f}, std = {std:.2f}")