class Flux:
def __init__(self, source_position: SkyCoord, e_min_reco: float = 0.2, e_max_reco: float = 10) -> None:
self.source_position = source_position
self.e_min_reco = e_min_reco * u.TeV
self.e_max_reco = e_max_reco * u.TeV
self.datasets_spectrum_joint_with_model = None
self.model = None
self.selected_model = None
self.energy_edges_flux_point = None
self.required_irf = ["aeff", "edisp"]
def initialize(self):
n_bin_per_decade_e_reco = 10
e_min_true = 0.005 * u.TeV
e_max_true = 40.0 * u.TeV
n_bin_per_decade_e_true = 20
e_min_flux_point = self.e_min_reco
e_max_flux_point = self.e_max_reco
n_bin_per_decade_flux_point = 10
n_bin_energy_reco = int(
(np.log10(self.e_max_reco.to_value(u.TeV)) - np.log10(self.e_min_reco.to_value(u.TeV)))
* n_bin_per_decade_e_reco
)
energy_axis_reco = MapAxis.from_edges(
np.logspace(
np.log10(self.e_min_reco.to_value(u.TeV)),
np.log10(self.e_max_reco.to_value(u.TeV)),
n_bin_energy_reco,
),
unit="TeV",
name="energy",
interp="log",
)
n_bin_energy_true = int(
(np.log10(e_max_true.to_value(u.TeV)) - np.log10(e_min_true.to_value(u.TeV))) * n_bin_per_decade_e_true
)
energy_axis_true = MapAxis.from_edges(
np.logspace(
np.log10(e_min_true.to_value(u.TeV)), np.log10(e_max_true.to_value(u.TeV)), n_bin_energy_true
),
unit="TeV",
name="energy_true",
interp="log",
)
on_radius = 0.2 * u.deg
on_region = CircleSkyRegion(center=self.source_position, radius=on_radius)
geom_on_region = RegionGeom.create(region=on_region, axes=[energy_axis_reco])
geom = WcsGeom.create(
skydir=self.source_position, npix=(200, 200), binsz=0.02, frame="icrs", axes=[energy_axis_reco]
)
geom_image = geom.to_image()
exclusion_radius = 0.35 * u.deg
exclude_region = CircleSkyRegion(center=self.source_position, radius=exclusion_radius)
exclusion_mask = ~geom_image.region_mask([exclude_region])
n_bin_energy_flux_point = int(
(np.log10(e_max_flux_point.to_value(u.TeV)) - np.log10(e_min_flux_point.to_value(u.TeV)))
* n_bin_per_decade_flux_point
)
energy_edges_flux_point = (
np.logspace(
np.log10(e_min_flux_point.to_value(u.TeV)),
np.log10(e_max_flux_point.to_value(u.TeV)),
n_bin_energy_flux_point,
)
* u.TeV
)
return geom_on_region, energy_axis_true, exclusion_mask, energy_axis_reco, energy_edges_flux_point
def get_datasets(self, obs_ids: list, obs_collection: Observations):
geom_on_region, energy_axis_true, exclusion_mask, energy_axis_reco, energy_edges_flux_point = (
self.initialize()
)
datasets_spectrum_joint = Datasets()
dataset_maker_spectrum = SpectrumDatasetMaker(
selection=["counts", "exposure", "edisp"], use_region_center=True
)
spectrum_dataset_empty = SpectrumDataset.create(geom=geom_on_region, energy_axis_true=energy_axis_true)
bkg_maker_spectrum = ReflectedRegionsBackgroundMaker(exclusion_mask=exclusion_mask)
for i, obs in enumerate(obs_collection):
dataset_spectrum = dataset_maker_spectrum.run(spectrum_dataset_empty.copy(name=f"obs-{i}".format(i)), obs)
dataset_on_off_spectrum = bkg_maker_spectrum.run(observation=obs, dataset=dataset_spectrum)
datasets_spectrum_joint.append(dataset_on_off_spectrum)
return datasets_spectrum_joint, energy_axis_reco, energy_edges_flux_point
def compute_models(self, energy_axis_reco):
reference_energy = energy_axis_reco.center[int(energy_axis_reco.nbin / 2)]
spectral_model = {
"power_law": PowerLawSpectralModel(
amplitude=1e-12 * u.Unit("cm-2 s-1 TeV-1"), index=2, reference=reference_energy
),
"log_parabola": LogParabolaSpectralModel(
amplitude=1e-12 * u.Unit("cm-2 s-1 TeV-1"), alpha=2, beta=0, reference=reference_energy
),
}
model = {}
for k in spectral_model:
model[k] = SkyModel(spectral_model=spectral_model[k], name=k)
return model, spectral_model
def select_model(self, khi2_flux) -> None:
"""Select a model fit for the flux if pvalue >5 sigma with a power law as H0 hypothesis."""
challenger_model = min(khi2_flux, key=khi2_flux.get)
LR_statistic = khi2_flux["power_law"] - khi2_flux[challenger_model]
dof = 1
p_val = chi2.sf(LR_statistic, dof)
significance_challenger = np.sqrt(2) * scipy.special.erfinv(1 - p_val)
if significance_challenger > 5.0:
self.selected_model = challenger_model
else:
self.selected_model = "power_law"
def get_time_unix_from_lst_epoch(self, time):
if self.version == "lstchain":
return time
LST_EPOCH = Time("2018-10-01T00:00:00", scale="utc")
time_lst_unix = Time(LST_EPOCH, format="unix", scale="utc").value
if type(time) is list:
return [t + time_lst_unix for t in time]
return time + time_lst_unix
def get_dataset_joint(self, obs_collection: Observations, obs_ids: list) -> None:
datasets_spectrum_joint, energy_axis_reco, energy_edges_flux_point = self.get_datasets(
obs_ids, obs_collection
)
model, spectral_model = self.compute_models(energy_axis_reco)
datasets_spectrum_joint_with_model = {}
fit_joint = {}
result_fit_joint = {}
khi2_flux = {}
for k in spectral_model.keys():
datasets_spectrum_joint_with_model[k] = datasets_spectrum_joint.copy()
datasets_spectrum_joint_with_model[k].models = [model[k]]
fit_joint[k] = Fit()
result_fit_joint[k] = fit_joint[k].run(datasets=datasets_spectrum_joint_with_model[k])
khi2_flux[k] = result_fit_joint[k].total_stat
for k in spectral_model.keys():
if not result_fit_joint[k].success:
logging.warning("Fit of " + k + " did not converge.")
else:
self.select_model(khi2_flux)
self.datasets_spectrum_joint_with_model = datasets_spectrum_joint_with_model
self.model = model
self.energy_edges_flux_point = energy_edges_flux_point
return datasets_spectrum_joint_with_model
def get_plot(self, data: SourceData, output_file: Path):
self.required_irf = data.required_irf
obs_collection = data.get_obs_collection()
obs_ids = data.get_obs_ids()
self.get_dataset_joint(obs_collection, obs_ids)
fpe = FluxPointsEstimator(
energy_edges=self.energy_edges_flux_point,
source=self.selected_model,
selection_optional="all",
n_sigma_ul=3,
reoptimize=False,
)
flux_points = fpe.run(datasets=self.datasets_spectrum_joint_with_model[self.selected_model])
flux_point_dataset_loop = FluxPointsDataset(data=flux_points, models=self.model[self.selected_model])
flux_point_dataset_loop.plot_spectrum()
plt.legend([f"{self.selected_model} fit", "Stat. err.", "Flux"], loc=0, frameon=True, prop={"size": 12})
plt.tight_layout()
if output_file:
plt.savefig(output_file)
# plt.show()
plt.close()