import gammalib

"""
Explanation: How to handle models from Python?
In this tutorial you will learn how to handle models from Python. This is done using the GammaLib classes GModels, GModel, GModelSpatial, GModelSpectral, and GModelTemporal. You can find all the information on the GammaLib classes in the Doxygen documentation.
In GammaLib there are some conventions about units. If not otherwise specified:


energies are in MeV


photon fluxes are in photons/cm2/s


differential photon fluxes are in photons/cm2/s/MeV


energy fluxes are in erg/cm2/s


To start we import the gammalib Python module.
End of explanation
"""


container = gammalib.GModels('$CTOOLS/share/models/crab.xml')

"""
Explanation: Reading and updating an existing XML model
Opening an XML model and parsing its content
An existing XML model can be read into a model container (GModels class).
End of explanation
"""


print(container)

"""
Explanation: The easiest way to inspect a model container is by using the function print that will display in the terminal its entire content.
End of explanation
"""


print('value', container['Crab']['Prefactor'].value())
print('error', container['Crab']['Prefactor'].error())
print('min value', container['Crab']['Prefactor'].min())
print('max value', container['Crab']['Prefactor'].max())
print('is free?', container['Crab']['Prefactor'].is_free())

"""
Explanation: The model container contains two models, for a total of ten parameters. The first model is a sky model named Crab, the second a background model. Each model has a spatial, a spectral, and a temporal component. Each model component has parameters. For each parameter you have a value, an error (in our example it is always zero because no fit to the data was done yet), an allowed value range for fitting the parameter to data, and, whenever relevant, units. In the parenthesis you can se if the parameter is free or fixed, and some more technical information, i.e, if the model has an analytical gradient w.r.t. this parameter implemented (can be fit to data quickly under certain conditions), and if internally the parameter is scaled (parameters should be scaled so that value/scale is of order unity for a more efficient handling by the optimiser during a fit to the data). Furthemore each model can be restricted to a specific instrument or observation. In the example the only restriction is that the CTA background model is used only for CTA data (Instruments ...............: CTA).
Model containers can be parsed using the model and parameter names if they are know to you. Below for example we directly access the properties of the Prefactor for the spectral model of the Crab source.
End of explanation
"""


print(container['Crab'].spectral())

"""
Explanation: You can access specific model components through the spatial, spectral, and temporal methods, for example:
End of explanation
"""


for s in range(container.size()):
    model = container[s]
    print('###### Name:', model.name())
    for i in range(model.spectral().size()):
        param = model.spectral()[i]
        print(param.name(),'=',param.value())
        

"""
Explanation: You can also blindly parse models and model parameters. Below we get names and values for the spectral parameters of all the models in the container.
End of explanation
"""


# calculate Crab spectral model flux in the 1-10 TeV energy range

# define energy bounds
emin = gammalib.GEnergy(1.,'TeV')
emax = gammalib.GEnergy(10.,'TeV')

# extract model from container'
crab = container['Crab']

# photon flux (cm-2 s-1)
flux = crab.spectral().flux(emin,emax)
print('Crab photon flux:', flux, 'cm-2 s-1')

# energy flux (cm-2 s-1)
eflux = crab.spectral().eflux(emin,emax)
print('Crab energy flux', eflux, 'erg cm-2 s-1')

# calculate Crab spatial model flux in 0.2 deg region centred on the source

# region centre = source position
centre = crab.spatial().dir()

# circular region
reg = gammalib.GSkyRegionCircle(centre,0.2)

# flux in the circular region (relative value)
flux = crab.spatial().flux(reg)
print('Crab flux in region:',flux,'(relative value)')

"""
Explanation: Finally there are convenience methods to compute for a model component the flux and energy flux over an energy interval (spectral components), and the flux within a circular region (spatial models).
End of explanation
"""


# change R.A.
crab['RA'].value(crab['RA'].value() + .1)
# free spatial parameters
crab['RA'].free()
crab['DEC'].free()
# restric range
crab['RA'].min(79.)
crab['RA'].max(89.)
crab['DEC'].min(21.5)
crab['DEC'].max(22.5)

# inspect model to verify changes
print(crab)

"""
Explanation: Spatial models can be normalised to unity. This is always the case for analytical models including point sources. Therefore we do expect to have a spatial model flux of 1 for a region encompassing the source.
Updating a model
Once you have a model container you can update its parameters, change model components and attributes, delete existing models or append new ones.
In the first example we change the position of the Crab source by 0.1 deg toward positive R.A. and we free the spatial parameters, but we restrict the allowed range around the known position.
End of explanation
"""


expplaw = gammalib.GModelSpectralExpPlaw()
expplaw['Prefactor'].value(crab['Prefactor'].value())
expplaw['Index'].value(crab['Index'].value())
expplaw['PivotEnergy'].value(crab['PivotEnergy'].value())
expplaw['CutoffEnergy'].value(50.e6) # value in MeV
crab.spectral(expplaw)

# inspect model to verify changes
print(crab)

"""
Explanation: In the second example we change the spectral model from a power law to an exponentially-cutoff power law, keeping the power law parameters to the original values, and adding a cutoff at 50 TeV.
End of explanation
"""


# remove existing background model
container.remove('CTABackgroundModel')

# open example background cube model
bkgcube_container = gammalib.GModels('$CTOOLS/share/models/bkg_cube.xml')

# append first model in new container to old container
container.append(bkgcube_container[0])

# inspect model to verify changes
print(container)

"""
Explanation: In the last example we remove the IRF background model, and append a cube background model from another container.
End of explanation
"""


container.save('my_crab.xml')

"""
Explanation: Save model changes to an XML file
You can use directly the model container in Python to pass it to ctools or cscripts, but you can also write it to an XML file.
End of explanation
"""


my_container = gammalib.GModels()

"""
Explanation: Creating a new model
New models can be created directly in Python. First, create a new empty model container.
End of explanation
"""


# define source direction
srcdir = gammalib.GSkyDir()
# set R.A. and Dec
srcdir.radec_deg(54.,-19.)

# Gaussian spatial component
spatial = gammalib.GModelSpatialRadialGauss(srcdir,0.5) # centre and radius in deg

"""
Explanation: To create a new source we need to define at least spatial and spectral components (if not specified the temporal component is taken to be constants). Refer to the Doxygen documentation for finding out all the models available in GammaLib and how to use them.
We first create a Gaussian spatial component.
End of explanation
"""


# define pivot energy
pivot = gammalib.GEnergy(1,'TeV')

# power law
spectral = gammalib.GModelSpectralPlaw(1.e-18,-2.5,pivot) # differential photon flux, index, pivot energy

"""
Explanation: We create a power-law spectral component.
End of explanation
"""


# create source
source = gammalib.GModelSky(spatial,spectral)
source.name('my_source')

# append to container
my_container.append(source)

"""
Explanation: We create the source and append it to the model container.
End of explanation
"""


# spectral correction
spectral = gammalib.GModelSpectralPlaw(1, 0, gammalib.GEnergy(1, 'TeV'))

# create background model
bkgmodel = gammalib.GCTAModelIrfBackground(spectral)
bkgmodel.name('Background')
bkgmodel.instruments('CTA')

# append to container
my_container.append(bkgmodel)

"""
Explanation: We also add a CTA IRF background model with power-law spectral correction.
End of explanation
"""


print(my_container)

"""
Explanation: We inspect the model we just created.
End of explanation
"""


my_container.save('mymodel.xml')

"""
Explanation: The container can be used in Python or written to disk in an XML file.
End of explanation
"""


