--- a/Spectrum_gammapy.py	Mon Feb 19 10:56:44 2024 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,266 +0,0 @@
-#!/usr/bin/env python
-# coding: utf-8
-
-# flake8: noqa
-
-import json
-import os
-import shutil
-
-import astropy.units as u
-import matplotlib.pyplot as plt
-import numpy as np
-from astropy.coordinates import Angle, SkyCoord
-from astropy.time import Time
-from gammapy.data import DataStore
-from gammapy.datasets import (
-    Datasets,
-    FluxPointsDataset,
-    SpectrumDataset,
-    SpectrumDatasetOnOff,
-)
-from gammapy.makers import (
-    ReflectedRegionsBackgroundMaker,
-    SafeMaskMaker,
-    SpectrumDatasetMaker,
-)
-
-# from gammapy.makers.utils import make_theta_squared_table
-from gammapy.maps import MapAxis, RegionGeom, WcsGeom
-from oda_api.data_products import ODAAstropyTable, PictureProduct
-from oda_api.json import CustomJSONEncoder
-from regions import CircleSkyRegion
-
-hess_data = "gammapy-datasets/1.1/hess-dl3-dr1/"
-if not (os.path.exists(hess_data)):
-    get_ipython().system("gammapy download datasets")   # noqa: F821
-
-data_store = DataStore.from_dir(hess_data)
-
-# src_name='Crab' #http://odahub.io/ontology#AstrophysicalObject
-# RA = 83.628700  # http://odahub.io/ontology#PointOfInterestRA
-# DEC = 22.014700 # http://odahub.io/ontology#PointOfInterestDEC
-src_name = "PKS 2155-304"
-RA = 329.716938  # http://odahub.io/ontology#PointOfInterestRA
-DEC = -30.225588  # http://odahub.io/ontology#PointOfInterestDEC
-T1 = "2000-10-09T13:16:00.0"  # http://odahub.io/ontology#StartTime
-T2 = "2022-10-10T13:16:00.0"  # http://odahub.io/ontology#EndTime
-Radius = 2.5  # http://odahub.io/ontology#AngleDegrees
-R_s = 0.5  # http://odahub.io/ontology#AngleDegrees
-
-Emin = 100.0  # http://odahub.io/ontology#Energy_GeV
-Emax = 10000.0  # http://odahub.io/ontology#Energy_GeV
-NEbins = 20  # http://odahub.io/ontology#Integer
-
-_galaxy_wd = os.getcwd()
-
-with open("inputs.json", "r") as fd:
-    inp_dic = json.load(fd)
-if "_data_product" in inp_dic.keys():
-    inp_pdic = inp_dic["_data_product"]
-else:
-    inp_pdic = inp_dic
-
-for vn, vv in inp_pdic.items():
-    if vn != "_selector":
-        globals()[vn] = type(globals()[vn])(vv)
-
-T1 = Time(T1, format="isot", scale="utc").mjd
-T2 = Time(T2, format="isot", scale="utc").mjd
-
-dates1 = data_store.obs_table["DATE-OBS"]
-dates2 = data_store.obs_table["DATE-END"]
-times1 = data_store.obs_table["TIME-OBS"]
-times2 = data_store.obs_table["TIME-END"]
-OBSIDs = data_store.obs_table["OBS_ID"]
-Tstart = []
-Tstop = []
-for i in range(len(dates1)):
-    Tstart.append(
-        Time(dates1[i] + "T" + times1[i], format="isot", scale="utc").mjd
-    )
-    Tstop.append(
-        Time(dates2[i] + "T" + times2[i], format="isot", scale="utc").mjd
-    )
-
-RA_pnts = np.array(data_store.obs_table["RA_PNT"])
-DEC_pnts = np.array(data_store.obs_table["DEC_PNT"])
-
-Coords_s = SkyCoord(RA, DEC, unit="degree")
-COORDS_pnts = SkyCoord(RA_pnts, DEC_pnts, unit="degree")
-seps = COORDS_pnts.separation(Coords_s).deg
-
-mask = np.where((seps < Radius) & (Tstart > T1) & (Tstop < T2))[0]
-OBSlist = []
-obs_ids = OBSIDs[mask]
-if len(obs_ids) == 0:
-    message = "No data found"
-    raise RuntimeError("No data found")
-obs_ids
-
-observations = data_store.get_observations(obs_ids)
-
-target_position = Coords_s
-on_region_radius = Angle(str(R_s) + " deg")
-on_region = CircleSkyRegion(center=target_position, radius=on_region_radius)
-skydir = target_position.galactic
-geom = WcsGeom.create(
-    npix=(150, 150), binsz=0.05, skydir=skydir, proj="TAN", frame="icrs"
-)
-
-Emin = 100.0  # http://odahub.io/ontology#Energy_GeV
-Emax = 10000.0  # http://odahub.io/ontology#Energy_GeV
-NEbins = 20  # http://odahub.io/ontology#Integer
-
-energy_axis = MapAxis.from_energy_bounds(
-    Emin * 1e-3,
-    Emax * 1e-3,
-    nbin=NEbins,
-    per_decade=True,
-    unit="TeV",
-    name="energy",
-)
-energy_axis_true = MapAxis.from_energy_bounds(
-    0.05, 100, nbin=20, per_decade=True, unit="TeV", name="energy_true"
-)
-
-geom = RegionGeom.create(region=on_region, axes=[energy_axis])
-dataset_empty = SpectrumDataset.create(
-    geom=geom, energy_axis_true=energy_axis_true
-)
-
-dataset_maker = SpectrumDatasetMaker(
-    containment_correction=True, selection=["counts", "exposure", "edisp"]
-)
-bkg_maker = ReflectedRegionsBackgroundMaker()
-safe_mask_masker = SafeMaskMaker(methods=["aeff-max"], aeff_percent=10)
-
-datasets = Datasets()
-
-for obs_id, observation in zip(obs_ids, observations):
-    dataset = dataset_maker.run(
-        dataset_empty.copy(name=str(obs_id)), observation
-    )
-    dataset_on_off = bkg_maker.run(dataset, observation)
-    # dataset_on_off = safe_mask_masker.run(dataset_on_off, observation)
-    datasets.append(dataset_on_off)
-
-print(datasets)
-
-from pathlib import Path
-
-path = Path("spectrum_analysis")
-path.mkdir(exist_ok=True)
-
-for dataset in datasets:
-    dataset.write(
-        filename=path / f"obs_{dataset.name}.fits.gz", overwrite=True
-    )
-
-datasets = Datasets()
-
-for obs_id in obs_ids:
-    filename = path / f"obs_{obs_id}.fits.gz"
-    datasets.append(SpectrumDatasetOnOff.read(filename))
-
-from gammapy.modeling import Fit
-from gammapy.modeling.models import ExpCutoffPowerLawSpectralModel, SkyModel
-
-spectral_model = ExpCutoffPowerLawSpectralModel(
-    amplitude=1e-12 * u.Unit("cm-2 s-1 TeV-1"),
-    index=2,
-    lambda_=0.1 * u.Unit("TeV-1"),
-    reference=1 * u.TeV,
-)
-model = SkyModel(spectral_model=spectral_model, name="crab")
-
-datasets.models = [model]
-
-fit_joint = Fit()
-result_joint = fit_joint.run(datasets=datasets)
-
-# we make a copy here to compare it later
-model_best_joint = model.copy()
-
-print(result_joint)
-
-display(result_joint.models.to_parameters_table())   # noqa: F821
-
-e_min, e_max = Emin * 1e-3, Emax * 1e-3
-energy_edges = np.geomspace(e_min, e_max, NEbins) * u.TeV
-
-from gammapy.estimators import FluxPointsEstimator
-
-fpe = FluxPointsEstimator(
-    energy_edges=energy_edges, source="crab", selection_optional="all"
-)
-flux_points = fpe.run(datasets=datasets)
-
-flux_points_dataset = FluxPointsDataset(
-    data=flux_points, models=model_best_joint
-)
-flux_points_dataset.plot_fit()
-# plt.show()
-plt.savefig("Spectrum.png", format="png", bbox_inches="tight")
-
-res = flux_points.to_table(sed_type="dnde", formatted=True)
-np.array(res["dnde"])
-
-bin_image = PictureProduct.from_file("Spectrum.png")
-from astropy.table import Table
-
-Emean = np.array(res["e_ref"])
-Emin = np.array(res["e_min"])
-Emax = np.array(res["e_max"])
-flux = Emean**2 * np.array(res["dnde"])
-flux_err = Emean**2 * np.array(res["dnde_err"])
-data = [Emean, Emin, Emax, flux, flux_err]
-names = (
-    "Emean[TeV]",
-    "Emin[TeV]",
-    "Emax[TeV]",
-    "Flux[TeV/cm2s]",
-    "Flux_error[TeV/cm2s]",
-)
-spec = ODAAstropyTable(Table(data, names=names))
-
-picture_png = bin_image  # http://odahub.io/ontology#ODAPictureProduct
-spectrum_astropy_table = spec  # http://odahub.io/ontology#ODAAstropyTable
-
-# output gathering
-_galaxy_meta_data = {}
-_oda_outs = []
-_oda_outs.append(
-    (
-        "out_Spectrum_gammapy_picture_png",
-        "picture_png_galaxy.output",
-        picture_png,
-    )
-)
-_oda_outs.append(
-    (
-        "out_Spectrum_gammapy_spectrum_astropy_table",
-        "spectrum_astropy_table_galaxy.output",
-        spectrum_astropy_table,
-    )
-)
-
-for _outn, _outfn, _outv in _oda_outs:
-    _galaxy_outfile_name = os.path.join(_galaxy_wd, _outfn)
-    if isinstance(_outv, str) and os.path.isfile(_outv):
-        shutil.move(_outv, _galaxy_outfile_name)
-        _galaxy_meta_data[_outn] = {"ext": "_sniff_"}
-    elif getattr(_outv, "write_fits_file", None):
-        _outv.write_fits_file(_galaxy_outfile_name)
-        _galaxy_meta_data[_outn] = {"ext": "fits"}
-    elif getattr(_outv, "write_file", None):
-        _outv.write_file(_galaxy_outfile_name)
-        _galaxy_meta_data[_outn] = {"ext": "_sniff_"}
-    else:
-        with open(_galaxy_outfile_name, "w") as fd:
-            json.dump(_outv, fd, cls=CustomJSONEncoder)
-        _galaxy_meta_data[_outn] = {"ext": "json"}
-
-with open(os.path.join(_galaxy_wd, "galaxy.json"), "w") as fd:
-    json.dump(_galaxy_meta_data, fd)
-print("*** Job finished successfully ***")