<jupyter_start><jupyter_text>This example notebook, uses the provided OGIP files, to plot the Light Curves One can follow the tutorial from gammapy here https://docs.gammapy.org/0.19/tutorials/analysis/time/light_curve.html<jupyter_code>%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec

from pathlib import Path

import astropy.units as u
from astropy.time import Time
from astropy.io import fits
from astropy.table import Table

import numpy as np

from gammapy.modeling import Fit
from gammapy.modeling.models import (
    PowerLawSpectralModel,
    LogParabolaSpectralModel,
    create_crab_spectral_model,
    SkyModel,
)
from gammapy.datasets import Datasets, SpectrumDatasetOnOff

from gammapy.estimators import LightCurveEstimator<jupyter_output><empty_output><jupyter_text>1. Get the OGIP files and make some selections if need be<jupyter_code>base_dir = "/fefs/aswg/workspace/analysis-school-2022/DL3/"
dir_path = "Crab_src_indep/"

ogip_path = Path(base_dir+dir_path+"OGIP/")

# Create the Paths if they do not exist already
ogip_path.mkdir(exist_ok=True)
# read all obs ids of the pha files in the given directory
obs_ids = []
pha_files = list(ogip_path.glob("obs_*_arf.fits.gz"))
    
for p in pha_files:
    run = int(p.name[4:8])
    obs_ids.append(run)
obs_ids = np.sort(np.array(obs_ids))
# Generate the Datasets object back from the OGIP files
datasets = Datasets()
for obs in obs_ids:
    file = ogip_path / f"obs_{obs}.fits.gz"
    datasets.append(SpectrumDatasetOnOff.read(file))<jupyter_output><empty_output><jupyter_text>Read the OGIP files and some metadata stored in it<jupyter_code>fits.open(ogip_path / f"obs_{obs_ids[1]}.fits.gz").info()
Table.read(ogip_path / f"obs_{obs_ids[1]}.fits.gz", hdu="SPECTRUM").meta
Table.read(ogip_path / f"obs_{obs_ids[1]}.fits.gz", hdu="REGION").meta
fits.open(ogip_path / f"obs_{obs_ids[1]}_bkg.fits.gz").info()
Table.read(ogip_path / f"obs_{obs_ids[1]}_bkg.fits.gz", hdu="SPECTRUM").meta
fits.open(ogip_path / f"obs_{obs_ids[1]}_arf.fits.gz").info()
Table.read(ogip_path / f"obs_{obs_ids[1]}_arf.fits.gz", hdu="SPECRESP")
fits.open(ogip_path / f"obs_{obs_ids[1]}_rmf.fits.gz").info()
Table.read(ogip_path / f"obs_{obs_ids[1]}_rmf.fits.gz", hdu="MATRIX")<jupyter_output><empty_output><jupyter_text>2. Get some parameters to use from the datasets<jupyter_code>t_start = []
t_stop = []
tot_time = []
obj_name = []

for obs in obs_ids:
    file = ogip_path / f"obs_{obs}.fits.gz"
    t = Table.read(file, hdu="GTI")
    r = Table.read(file, hdu="REGION").meta
    
    t_start.append(t["START"][0])
    t_stop.append(t["STOP"][0])
    tot_time.append(t["STOP"][0] - t["START"][0])
    
    obj_name.append(r["OBJECT"])
print(obs_ids)

t_start = np.sort(np.array(t_start))
t_stop = np.sort(np.array(t_stop))
tot_time = np.array(tot_time)
print(tot_time.sum()/3600)

obj_name = np.unique(np.array(obj_name))[0]
e_reco_edges = datasets[0].counts.geom.axes["energy"].edges
# Energy range for estimating the Light Curve.
# One can play with different energy bins to check the different LCs

# Using these values to check with Crab reference flux of MAGIC JHEAP 2015
e_fit_min = 0.3 * u.TeV
e_fit_max = 100 * u.TeV<jupyter_output><empty_output><jupyter_text>3. Get Pivot energy to fix the reference energy and define the Spectrum Model<jupyter_code># Find pivot (decorrelation) energy for a Power Law model to get the reference energy for Log Parabola model
def get_pivot_energy(datasets, e_ref, e_edges, obj_name):
    """
    Using Power Law spectral model with the given reference energy and
    get the decorrelation energy of the fit, within the fit energy range, e_edges
    """
    spectral_model = PowerLawSpectralModel(
        index=2, amplitude=2e-11 * u.Unit("cm-2 s-1 TeV-1"), reference=e_ref
    )
    model = SkyModel(spectral_model=spectral_model, name=obj_name)
    model_check = model.copy()

    # Stacked dataset method
    stacked_dataset = Datasets(datasets).stack_reduce()
    stacked_dataset.models = model_check

    fit_stacked = Fit()
    result_stacked = fit_stacked.run(datasets=stacked_dataset)
    
    return model_check.spectral_model.pivot_energy

# Using a reference energy close to the expected decorrelation energy
ref = get_pivot_energy(datasets, 0.4 * u.TeV, e_reco_edges, obj_name)
print(ref.to_value(u.GeV))
# Final spectral model of Log Parabola, to be used for estimating the LC.
# One can try different Spectral Models as well.
# Be careful in the choice of Spectral Model being used for the 2 examples presented here

# Using the same values as we get from plotting SED, from the other notebook.
# Crab
spectral_model_lp = LogParabolaSpectralModel(
        amplitude = 1.75e-10 * u.Unit('cm-2 s-1 TeV-1'),
        reference = ref,
        alpha = 2.325 * u.Unit(''),
        beta = 0.0697 * u.Unit('')
)
model_lp = SkyModel(spectral_model=spectral_model_lp, name=obj_name)
# Use the appropriate models, as per the selection of the source/dataset
params=model_lp.to_dict()['spectral']['parameters']
params<jupyter_output><empty_output><jupyter_text>4. Estimate the LC<jupyter_code>lc_maker_1d = LightCurveEstimator(
    energy_edges=[e_fit_min, e_fit_max], 
    reoptimize=False, # Re-optimizing other free model parameters (not belonging to the source)
    source=obj_name, 
    selection_optional="all" # Estimates asymmetric errors, upper limits and fit statistic profiles
)
# Need the time start and end
t_start = Time(t_start, format='unix')
t_stop = Time(t_stop, format='unix')

t_day = np.unique(np.rint(t_start.mjd))

# To make the range night-wise, keep the MJD range in half integral values
t_range = [Time([t-0.5, t+0.5], format="mjd", scale="utc") for t in t_day]

lc_maker_night_wise = LightCurveEstimator(
    energy_edges=[e_fit_min, e_fit_max], 
    time_intervals=t_range,
    reoptimize=False, 
    source=obj_name,
    selection_optional="all"
)
# Assigning the model for each dataset
for data in datasets:
    data.models = model_lp

lc_1d = lc_maker_1d.run(datasets)
lc_night = lc_maker_night_wise.run(datasets)
# Check the various column data of the Light Curve object
lc_1d.to_table(sed_type="flux", format="lightcurve")
# If there are more than 1 night of data, one can see the integrated light curve for each night
lc_night.to_table(sed_type="flux", format="lightcurve")<jupyter_output><empty_output><jupyter_text>5. Plot the Light Curve<jupyter_code># Calculate & plot Crab reference flux
# https://doi.org/10.1016/j.jheap.2015.01.002
crab = create_crab_spectral_model("magic_lp")
crab.amplitude.error = 0.03e-11 * u.Unit("cm-2 s-1 TeV-1")
crab.alpha.error = 0.01
crab.beta.error = 0.01/np.log(10)

flux_crab, flux_crab_error = crab.integral_error(e_fit_min, e_fit_max)
print(flux_crab, flux_crab_error)
fig_lc = plt.figure(figsize=(8,10))

gs2 = GridSpec(10, 5)

gs2.update(wspace=0.4)
args1 = [gs2[:5,:]]
args2 = [gs2[5:,:]]

fig_gs1 = fig_lc.add_subplot(*args1)
fig_gs2 = fig_lc.add_subplot(*args2, sharey=fig_gs1)

lc_1d.plot(
    ax=fig_gs1,
    sed_type="flux",
    marker="o", label="LST-1"
)
fig_gs1.axhline(
    flux_crab.to_value("cm-2 s-1"), c='red', ls='--', 
    label='Crab (MAGIC, JHEAp 2015)'
)
fig_gs1.axhspan(
    (flux_crab - flux_crab_error).to_value("cm-2 s-1"), 
    (flux_crab + flux_crab_error).to_value("cm-2 s-1"), 
    alpha=0.2, color='tab:orange'
)
fig_gs1.get_xaxis().set_ticklabels([])
fig_gs1.grid(which='both')
fig_gs1.set_title(
    f'LC LST-1 {obj_name}: {e_fit_min} < E < {e_fit_max} \nRun-wise {tot_time.sum()/3600:.2f} hrs, night-wise {len(t_day)} nights'
)
fig_gs1.legend()
fig_gs1.get_yaxis().get_offset_text().set_position((-0.06,1))

lc_night.plot(
    ax=fig_gs2,
    sed_type="flux",
    axis_name="time",
    marker="o", label="LST-1"
)
fig_gs2.axhline(
    flux_crab.to_value("cm-2 s-1"), c='red', ls='--', 
    label='Crab (MAGIC, JHEAp 2015)'
)
fig_gs2.axhspan(
    (flux_crab - flux_crab_error).to_value("cm-2 s-1"), 
    (flux_crab + flux_crab_error).to_value("cm-2 s-1"), 
    alpha=0.2, color='tab:orange'
)

fig_gs2.grid(which='both')
fig_gs2.legend()<jupyter_output>No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
{'fraction_non_alphanumeric': 0.1271999029008375, 'fraction_numerical': 0.021968685520087388, 'mean_word_length': 4.898353614889048, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 28, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '372470', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3248, 'n_tokens_neox': 3114, 'n_words': 695, 'path': 'notebooks/plot_lc_from_DL4.ipynb', 'repo_name': 'mireianievas/cta-lstchain', 'size': '8239'}