<jupyter_start><jupyter_text>load the VEGAS fits first to get the energy bands<jupyter_code>ff = fits.open(("../VEGAS/CrabDL3{0:s}Point_s6.fits".format(cuts)))
points = ff[1].data

points['Energy'] *= u.TeV
for f in ['Flux', 'FluxError', 'FluxUpperLimit']:
    points[f] = (points[f]*u.Unit('TeV-1 m-2 s-1')).to(u.Unit('TeV-1 cm-2 s-1'))

fit = ff[2].data
cov = ff[3].data

# this trickery allows us to reject points after the first non-significant one
cond = (points['Significance'] > 2) * (points['On'] > 2)
for i,c in enumerate(cond):
    if c == False:
        cond[i:] = False<jupyter_output><empty_output><jupyter_text>Now do the gammapy analysis<jupyter_code>data_store = DataStore.from_dir(("../VEGAS/DL3{0:s}/".format(cuts)))
data_store.info()
obs_id = data_store.obs_table['OBS_ID'].data
print(("Use observations {}".format(obs_id)))

obs_list = data_store.obs_list(obs_id)
crab_pos = SkyCoord(83.633083, 22.0145, unit='deg')

if cuts == "soft":
    radius = 0.173 * u.deg
else:
    radius = 0.1 * u.deg

on_region = CircleSkyRegion(crab_pos, radius)

model = models.LogParabola(
    alpha=2.5,
    beta=0.1,
    amplitude=3e-11 * u.Unit('cm-2 s-1 TeV-1'),
    reference=1 * u.TeV,
)

flux_point_binning = EnergyBounds.equal_log_spacing(
    points['EBinLowEdge'][cond][0], points['EBinHiEdge'][cond][-1],
    len(points[cond]), u.TeV)
config = dict(
    outdir = None,
    background = dict(
        on_region=on_region,
#         exclusion_mask=exclusion_mask,
        min_distance = 0.5 * u.deg,
    ),
    extraction = dict(containment_correction=False),
    fit = dict(
        model=model,
        stat='wstat',
        forward_folded=True,
        fit_range = flux_point_binning[[0, -1]]
    ),
    fp_binning=flux_point_binning
)
ana = SpectrumAnalysisIACT(
    observations=obs_list,
    config=config,
)
ana.run()
print(ana.fit.result[0])
ax = ana.spectrum_result.plot(
    energy_range=ana.fit.fit_range,
    energy_power=2,
    flux_unit='TeV-1 cm-2 s-1',
    fig_kwargs=dict(figsize = (8,8)),
)
ax[1].set_ylim(-2, 2)
ax[0].set_ylim(5e-12, 1e-10)
fig, ax, cbar = ana.background_estimator.plot()
ax.set_xlim(-70, 70)
ax.set_ylim(0, 100)
vegas.parameters.set_parameter_errors?
vegas = gammapy.spectrum.models.LogParabola(
    amplitude=(fit['Normalization'][0] * u.Unit('1 / (m2 s TeV)')).to(
        u.Unit('1 / (cm2 s TeV)')),
    reference=fit['E0'][0] * u.TeV,
    alpha=fit['Index'][0],
    beta=fit['Cutoff'][0])

vegas.parameters.set_parameter_errors({
    "amplitude":
    fit['Normalization'][1] * u.Unit('1 / (m2 s TeV)'),
    "alpha":
    fit['Index'][1] * u.Unit(''),
    "beta":
    fit['Cutoff'][1] * u.Unit('')
})
# we need to convert this into a form that is understood by gammapy, grrrr

cov2 = np.array(cov.copy())
cov2 = np.array(cov2.tolist())
cov2 = np.insert(np.insert(cov2, 1, 0, axis=0), 1, 0, axis=1)
cov2[:, 0] = cov2[:, 0] * 1e-4
cov2[0, :] = cov2[0, :] * 1e-4

vegas.parameters.covariance = cov2
vegas.parameters.to_table()
ax = ana.spectrum_result.plot(
    energy_range=ana.fit.fit_range,
    energy_power=2,
    flux_unit='TeV-1 cm-2 s-1',
    fig_kwargs=dict(figsize=(8, 8)),
    butterfly_kwargs=dict(facecolor="C0"),
    point_kwargs=dict(c="C0"))

vegas.plot(
    energy_range=ana.fit.fit_range,
    ax=ax[0],
    energy_power=2,
    flux_unit='TeV-1 cm-2 s-1',
)
vegas.plot_error(
    energy_range=ana.fit.fit_range,
    ax=ax[0],
    energy_power=2,
    flux_unit='TeV-1 cm-2 s-1',
    facecolor="C1",
)

ax[0].errorbar(
    points['Energy'][cond],
    points['Flux'][cond] * points['Energy'][cond]**2,
    yerr=points['FluxError'][cond] * points['Energy'][cond]**2,
    c="C1",
    marker="d",
    ls="")

ax[1].set_ylim(-2, 2)
ax[0].set_ylim(5e-12, 1e-10)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18723290598290598, 'fraction_numerical': 0.039797008547008544, 'mean_word_length': 3.934123847167325, 'pattern_counts': {'":': 4, '<': 9, '<?xml version=': 0, '>': 11, 'https://': 0, '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': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '81889', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 1611, 'n_tokens_neox': 1561, 'n_words': 257, 'path': 'notebooks/VERITAS_spectrum.ipynb', 'repo_name': 'VERITAS-Observatory/V2DL3', 'size': '3744'}