import gammalib
import ctools
import cscripts

"""
Explanation: Performing a classical On/Off analysis
In this tutorial you will learn to perform a classical On/Off analysis of the data.
We start by importing the gammalib, ctools, and cscripts Python modules.
End of explanation
"""


%matplotlib inline
import matplotlib.pyplot as plt

"""
Explanation: We will also use matplotlib to display the results.
End of explanation
"""


import sys
import os
sys.path.append(os.environ['CTOOLS']+'/share/examples/python/')

"""
Explanation: And finally we add to our path the directory containing the example plotting scripts provided with the ctools installation.
End of explanation
"""


obsfile = 'obs_crab_selected.xml'
phagen  = cscripts.csphagen()
phagen['inobs']         = obsfile
phagen['inmodel']       = 'NONE' # assume that the source is pointlike
phagen['ebinalg']       = 'LOG'
phagen['emin']          = 0.66
phagen['emax']          = 30.0
phagen['enumbins']      = 20
phagen['coordsys']      = 'CEL'
phagen['ra']            = 83.63
phagen['dec']           = 22.01
phagen['rad']           = 0.2
phagen['bkgmethod']     = 'REFLECTED'
phagen['use_model_bkg'] = False
phagen['maxoffset']     = 2.0
phagen['stack']         = True
phagen['outobs']        = 'obs_crab_onoff.xml'
phagen['outmodel']      = 'onoff_model.xml'
phagen.execute()

"""
Explanation: Classical analysis without background model
Preparation of the On/Off binned data
The first step for a classical analysis is to bin the events in bins of energy for an On (signal) region, and one or several Off (background regions).
We will use the event selection performed in the previous tutorial.
End of explanation
"""


print(phagen.obs()[0])

"""
Explanation: Note that we have used the hidden parameter use_model_bkg to specify that no background model should be used in the computation of the background normalisation factor, which is simply assumed to be the ratio of the solid angles of On/Off regions.
Let us peek at the resulting set of On/Off observations.
End of explanation
"""


print(phagen.obs().models())

"""
Explanation: Note that the instrument is now HESSOnOff, to signify that we are dealing with On/Off observations. Since we have no background model the statistic has been set to WSTAT. The script has also created a model for the On region.
End of explanation
"""


phagen.obs().models()['Dummy'].name('Crab')

"""
Explanation: Since we do not have specified any input model the script has placed a point source named Dummy, at the center of the On region with a power law spectrum. We will rename this source as Crab. At this point you can also modify the spectral model, e.g., to use more appropriate parameter values or a different spectral shape.
End of explanation
"""


like = ctools.ctlike(phagen.obs())
like.run()

"""
Explanation: Likelihood fit and residuals
We can now fit this model to the data.
End of explanation
"""


print(like.opt())
print(like.obs().models())

"""
Explanation: Let's look at the results from the optimisation and the fitted model.
End of explanation
"""


residuals = 'residuals_classical_nobkg.fits'
res = cscripts.csresspec(like.obs())
res['components'] = True
res['algorithm']  = 'SIGNIFICANCE'
res['outfile']    = residuals
res.execute()
from show_residuals import plot_residuals
plot_residuals(residuals,'',0)

"""
Explanation: The fit has properly converged and the results are consistent with those obtained with the unbinned analysis.
We will also check the spectral residuals.
End of explanation
"""


srcmodel   = os.environ['CTOOLS']+'/share/models/crab_nobkg.xml'
bkgmodel   = 'bkgmodel.xml'
startmodel = 'model_classical_bkg.xml'
merge = cscripts.csmodelmerge()
merge['inmodels'] = srcmodel + ' ' + bkgmodel
merge['outmodel'] = startmodel
merge.execute()

"""
Explanation: Classical analysis with a background model
Preparation of the On/Off binned data
If a background model is available this can be used to compute the expected background rates for different observations as a function of energy and normalise the number of background counts between On and Off regions. Here we will use the background model derived using csbkgmodel.
End of explanation
"""


phagen = cscripts.csphagen()
phagen['inobs']         = obsfile
phagen['inmodel']       = startmodel
phagen['srcname']       = 'Crab'
phagen['ebinalg']       = 'LOG'
phagen['emin']          = 0.66
phagen['emax']          = 30.0
phagen['enumbins']      = 20
phagen['coordsys']      = 'CEL'
phagen['ra']            = 83.63
phagen['dec']           = 22.01
phagen['rad']           = 0.2
phagen['bkgmethod']     = 'REFLECTED'
phagen['use_model_bkg'] = True
phagen['maxoffset']     = 2.0
phagen['stack']         = False # treat observations separately in joint likelihood fit
phagen.run()

"""
Explanation: Now we can run csphagen using this model including the background in input. Note that since the background model is defined per observation we will not stack the observations together.
End of explanation
"""


like = ctools.ctlike(phagen.obs())
like.run()
print(like.opt())
print(like.obs().models())

"""
Explanation: Likelihood fit and residuals
We will now fit the model to the data and look at the results.
End of explanation
"""


residuals = 'residuals_classical_bkgmod.fits'
res = cscripts.csresspec(like.obs())
res['components'] = True
res['algorithm']  = 'SIGNIFICANCE'
res['outfile']    = residuals
res['stack']      = False
res.execute()
from show_residuals import plot_residuals
plot_residuals(residuals,'',0)

"""
Explanation: Results for the spectrum of the Crab nebula are once again consistent with those obtained before. Note how the spatial parameters of the background model, although formally free, have not been optimised (their errors are 0), since we are running a spectral analysis only.
Since we have provided a background model to csphagen the statistic was set by default to CSTAT, i.e. the background model was also used in the fit. You may also specify as statistic in ctlike WSTAT. In this case the spectral model for the background would be derived from the data during the fit, and the impact of the model used above would extend only to the computation of the background normalisation factor between On and Off regions. Using an a-priori background model on the other hand can be advantageous to avoid intrinsic biases of WSTAT when the counting statistics are limited.
Let's check the residuals for the first of the 4 observations.
End of explanation
"""


