Coverage for src / lstautorta / Spectra.py: 0%
189 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
20import scipy
21from astropy.coordinates import SkyCoord
22from matplotlib import style
24style.use("tableau-colorblind10")
25import argparse
26from pathlib import Path
28from gammapy.data import DataStore
29from gammapy.datasets import Datasets, FluxPointsDataset, MapDataset, SpectrumDataset
30from gammapy.estimators import FluxPointsEstimator
31from gammapy.makers import (
32 MapDatasetMaker,
33 ReflectedRegionsBackgroundMaker,
34 SafeMaskMaker,
35 SpectrumDatasetMaker,
36)
37from gammapy.maps import MapAxis, RegionGeom, WcsGeom
38from gammapy.modeling import Fit
39from gammapy.modeling.models import (
40 ExpCutoffPowerLawSpectralModel,
41 LogParabolaSpectralModel,
42 PowerLawSpectralModel,
43 SkyModel,
44)
45from regions import CircleSkyRegion
47parser = argparse.ArgumentParser(
48 description="Automatic Script for the DL1 check", formatter_class=argparse.ArgumentDefaultsHelpFormatter
49)
50parser.add_argument("-d", "--directory", default="/fefs/onsite/pipeline/rta/data/", help="Directory for data")
51parser.add_argument("-da", "--date", default="20230705", help="Date of the run to check")
52parser.add_argument("-r", "--run-id", default="13600", help="run id to check")
53parser.add_argument("-add", "--add-string", default="", 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 = f"{source_name}_{cut_type}"
70source_position = SkyCoord(ra=config["right_ascension"], dec=config["declination"], unit="deg", frame="icrs")
71max_offset_run = 5 * u.deg
72work_directory = location_data
73path_plot = Path(work_directory + "/../plots")
74print(work_directory + "/../plots")
75path_plot.mkdir(exist_ok=True)
76path_background = Path(work_directory + "/../plots")
77path_background.mkdir(exist_ok=True)
79on_radius = 0.2 * u.deg
80exclusion_radius = 0.35 * u.deg
82e_min_reco = 0.02 * u.TeV
83e_max_reco = 10.0 * u.TeV
84n_bin_per_decade_e_reco = 10
86e_min_true = 0.005 * u.TeV
87e_max_true = 40.0 * u.TeV
88n_bin_per_decade_e_true = 20
90e_min_flux_point = 0.02 * u.TeV
91e_max_flux_point = 10.0 * u.TeV
92n_bin_per_decade_flux_point = 5
95def plot_excess(dataset_to_plot, ax=None, kwargs_excess=None, kwargs_npred_signal=None, **kwargs):
96 from gammapy.stats import CashCountsStatistic, WStatCountsStatistic
98 region = None
100 kwargs_excess = kwargs_excess or {}
101 kwargs_npred_signal = kwargs_npred_signal or {}
103 counts, npred = dataset_to_plot.counts.copy(), dataset_to_plot.npred()
105 if dataset_to_plot.mask is None:
106 mask = dataset_to_plot.copy()
107 mask.data = 1
108 else:
109 mask = dataset_to_plot.mask
110 counts *= mask
111 npred *= mask
113 if dataset_to_plot.stat_type == "wstat":
114 counts_off = dataset_to_plot.counts_off * mask
115 with np.errstate(invalid="ignore"):
116 alpha = dataset_to_plot.alpha * mask
118 stat = WStatCountsStatistic(
119 n_on=counts,
120 n_off=counts_off,
121 alpha=alpha,
122 # mu_sig=npred,
123 )
124 elif dataset_to_plot.stat_type == "cash":
125 stat = CashCountsStatistic(counts.data, npred.data)
127 yerr = np.zeros((2,) + counts.data.shape)
128 yerr[0], yerr[1] = -stat.compute_errn(), stat.compute_errp()
129 yerr = stat.error
131 plot_kwargs = kwargs.copy()
132 plot_kwargs.update(kwargs_excess)
133 plot_kwargs.setdefault("label", "Excess counts")
134 ax = dataset_to_plot.excess.plot(ax, yerr=yerr, **plot_kwargs)
136 plot_kwargs = kwargs.copy()
137 plot_kwargs.update(kwargs_npred_signal)
138 plot_kwargs.setdefault("label", "Predicted signal counts")
139 dataset_to_plot.npred_signal().plot_hist(ax, **plot_kwargs)
141 ax.legend(numpoints=1)
142 return ax
145def plot_fit(
146 dataset_to_plot,
147 ax_spectrum=None,
148 ax_residuals=None,
149 kwargs_spectrum=None,
150 kwargs_residuals=None,
151):
152 from gammapy.datasets.utils import get_axes
153 from matplotlib.gridspec import GridSpec
155 gs = GridSpec(7, 1)
156 ax_spectrum, ax_residuals = get_axes(
157 ax_spectrum,
158 ax_residuals,
159 8,
160 7,
161 [gs[:5, :]],
162 [gs[5:, :]],
163 kwargs2={"sharex": ax_spectrum},
164 )
165 kwargs_spectrum = kwargs_spectrum or {}
166 kwargs_residuals = kwargs_residuals or {}
168 plot_excess(dataset_to_plot, ax_spectrum, **kwargs_spectrum)
170 dataset_to_plot.plot_residuals_spectral(ax_residuals, **kwargs_residuals)
172 method = kwargs_residuals.get("method", "diff")
173 label = dataset_to_plot._residuals_labels[method]
174 ax_residuals.set_ylabel(f"Residuals\n{label}")
176 return ax_spectrum, ax_residuals
179data_store = DataStore.from_dir(location_data)
180obs_ids = data_store.obs_table[source_position.separation(data_store.obs_table.pointing_radec) < max_offset_run][
181 "OBS_ID"
182]
183obs_collection = data_store.get_observations(obs_ids, required_irf="point-like")
185on_region = CircleSkyRegion(center=source_position, radius=on_radius)
186exclude_region = CircleSkyRegion(center=source_position, radius=exclusion_radius)
188n_bin_energy_reco = int(
189 (np.log10(e_max_reco.to_value(u.TeV)) - np.log10(e_min_reco.to_value(u.TeV))) * n_bin_per_decade_e_reco
190)
191energy_axis_reco = MapAxis.from_edges(
192 np.logspace(np.log10(e_min_reco.to_value(u.TeV)), np.log10(e_max_reco.to_value(u.TeV)), n_bin_energy_reco),
193 unit="TeV",
194 name="energy",
195 interp="log",
196)
197n_bin_energy_true = int(
198 (np.log10(e_max_true.to_value(u.TeV)) - np.log10(e_min_true.to_value(u.TeV))) * n_bin_per_decade_e_true
199)
200energy_axis_true = MapAxis.from_edges(
201 np.logspace(np.log10(e_min_true.to_value(u.TeV)), np.log10(e_max_true.to_value(u.TeV)), n_bin_energy_true),
202 unit="TeV",
203 name="energy_true",
204 interp="log",
205)
207geom_on_region = RegionGeom.create(region=on_region, axes=[energy_axis_reco])
208geom = WcsGeom.create(
209 skydir=source_position,
210 npix=(200, 200),
211 binsz=0.05,
212 frame="icrs",
213 axes=[energy_axis_reco],
214)
216geom_image = geom.to_image()
217exclusion_mask = ~geom_image.region_mask([exclude_region])
220dataset_maker_spectrum = SpectrumDatasetMaker(selection=["counts", "exposure", "edisp"], use_region_center=True)
221spectrum_dataset_empty = SpectrumDataset.create(geom=geom_on_region, energy_axis_true=energy_axis_true)
222bkg_maker_spectrum = ReflectedRegionsBackgroundMaker(exclusion_mask=exclusion_mask)
224dataset_stacked_map = MapDataset.create(geom=geom, name=source_name + "_stacked")
225map_dataset_maker = MapDatasetMaker(selection=["counts"])
226map_maker_safe_mask = SafeMaskMaker(methods=["offset-max"], offset_max="4 deg")
228datasets_spectrum_joint = Datasets()
229for obs in obs_collection:
230 dataset_spectrum = dataset_maker_spectrum.run(spectrum_dataset_empty.copy(name=f"obs-{obs.obs_id}"), obs)
231 dataset_on_off_spectrum = bkg_maker_spectrum.run(observation=obs, dataset=dataset_spectrum)
232 datasets_spectrum_joint.append(dataset_on_off_spectrum)
234 map_cutout = dataset_stacked_map.cutout(obs.pointing_radec, width="6.5 deg")
235 map_dataset = map_dataset_maker.run(map_cutout, obs)
236 map_dataset = map_maker_safe_mask.run(map_dataset, obs)
237 dataset_stacked_map.stack(map_dataset)
239info_table = datasets_spectrum_joint.info_table(cumulative=True)
241reference_energy = energy_axis_reco.center[int(energy_axis_reco.nbin / 2)]
243spectral_model = {
244 "power_law": PowerLawSpectralModel(
245 amplitude=1e-12 * u.Unit("cm-2 s-1 TeV-1"), index=2, reference=reference_energy
246 ),
247 "log_parabola": LogParabolaSpectralModel(
248 amplitude=1e-12 * u.Unit("cm-2 s-1 TeV-1"), alpha=2, beta=0, reference=reference_energy
249 ),
250 "power_law_exponantial_cutoff": ExpCutoffPowerLawSpectralModel(
251 amplitude=1e-12 * u.Unit("cm-2 s-1 TeV-1"),
252 index=2,
253 lambda_=1.0 / reference_energy,
254 reference=reference_energy,
255 ),
256}
257model = {}
258for k in spectral_model:
259 model[k] = SkyModel(spectral_model=spectral_model[k], name=source_name + "_" + k)
261datasets_spectrum_joint_with_model = {}
262datasets_spectrum_joint_with_model_stacked = {}
263fit_joint = {}
264result_fit_joint = {}
265for k in spectral_model:
266 datasets_spectrum_joint_with_model[k] = datasets_spectrum_joint.copy()
267 datasets_spectrum_joint_with_model[k].models = [model[k]]
268 fit_joint[k] = Fit()
269 result_fit_joint[k] = fit_joint[k].run(datasets=datasets_spectrum_joint_with_model[k])
271 datasets_spectrum_joint_with_model_stacked[k] = datasets_spectrum_joint_with_model[k].stack_reduce()
272 datasets_spectrum_joint_with_model_stacked[k].models = [model[k]]
274if result_fit_joint["power_law_exponantial_cutoff"].total_stat > result_fit_joint["log_parabola"].total_stat:
275 chalenger_model = "log_parabola"
276 print("Log parabola perform better than power law than exponantial cut off")
277else:
278 chalenger_model = "power_law_exponantial_cutoff"
279 print("Power law than exponantial cut off perform better than log parabola")
282# significance_chalenger = np.sqrt((result_fit_joint['power_law'].total_stat-result_fit_joint[chalenger_model].total_stat))
283print("\n\nLikelihood Ratio test\n")
284LR_statistic = result_fit_joint["power_law"].total_stat - result_fit_joint[chalenger_model].total_stat
285dof = 1
286p_val = scipy.stats.chi2.sf(LR_statistic, dof)
288significance_chalenger = np.sqrt(2) * scipy.special.erfinv(1 - p_val)
289print(f"The significance of the {chalenger_model} model over power law is {significance_chalenger:.5f} sigma")
291if significance_chalenger > 5.0:
292 selected_model = chalenger_model
293 print("The spectral model " + chalenger_model + " is selected as significant")
294else:
295 selected_model = "power_law"
296 print("The spectral model " + chalenger_model + " is rejected as not significant over power law")
297 print("The selected spectral model is the power law")
299ax1, ax2 = plot_fit(datasets_spectrum_joint_with_model_stacked[k])
300excess = datasets_spectrum_joint_with_model_stacked[k].excess.data
301n_on = datasets_spectrum_joint_with_model_stacked[k].counts.data
302n_off = datasets_spectrum_joint_with_model_stacked[k].counts_off.data
303alpha = datasets_spectrum_joint_with_model_stacked[k].alpha.data
304low_excess = excess - np.sqrt(n_on + alpha * alpha * n_off)
305low_excess = low_excess[low_excess > 0.0]
306ax1.set_ylim(bottom=np.min(0.8 * low_excess))
308print(datasets_spectrum_joint_with_model[selected_model].models)
310n_bin_energy_flux_point = int(
311 (np.log10(e_max_flux_point.to_value(u.TeV)) - np.log10(e_min_flux_point.to_value(u.TeV)))
312 * n_bin_per_decade_flux_point
313)
314energy_edges_flux_point = (
315 np.logspace(
316 np.log10(e_min_flux_point.to_value(u.TeV)),
317 np.log10(e_max_flux_point.to_value(u.TeV)),
318 n_bin_energy_flux_point,
319 )
320 * u.TeV
321)
323fpe = FluxPointsEstimator(
324 energy_edges=energy_edges_flux_point,
325 source=source_name + "_" + selected_model,
326 selection_optional="all",
327 n_sigma_ul=3,
328 reoptimize=False,
329)
330flux_points = fpe.run(datasets=datasets_spectrum_joint_with_model[selected_model])
332flux_points.to_table(sed_type="e2dnde", formatted=True)
334plt.figure(figsize=(8, 5))
335ax = flux_points.plot(sed_type="e2dnde", color="darkorange")
336flux_points.plot_ts_profiles(ax=ax, sed_type="e2dnde")
337plt.savefig(os.path.join(path_plot, f"{filename_output}__flux_point.png"), dpi=300)
339flux_points_dataset = FluxPointsDataset(data=flux_points, models=model[selected_model])
340flux_points_dataset.plot_fit()
341plt.savefig(os.path.join(path_plot, f"{filename_output}__spectra.png"), dpi=300)