%matplotlib inline
%matplotlib notebook
import numpy as np
from threeML import *

get_available_plugins()

"""
Explanation: 3ML with Fermi GBM TTE  and LAT LLE Data
Purpose
This demonstrates of you can use LLE data directly in 3ML
For more detail on the Fermi TTE plugin, check out its demo.
Let's check it out!
Import 3ML as always to make sure you have the plugin.
End of explanation
"""


# os.path.join is a way to generate system-independent
# paths (good for unix, windows, Mac...)
trigger_number = 'bn080916009'

data_dir_gbm = os.path.join('gbm',trigger_number)
gbm_data = download_GBM_trigger_data(trigger_number,detectors=['n3','n4','b0'],destination_directory=data_dir_gbm,compress_tte=True)


data_dir_lle = os.path.join('lat')

lle_data = download_LLE_trigger_data(trigger_number,destination_directory=data_dir_lle)

src_selection = "0.-71."

"""
Explanation: Obtaining data
Use the LLE and GBM data downloaders to get the data:
End of explanation
"""


lle = FermiLATLLELike("LLE",
                      os.path.join(data_dir_lle, "gll_lle_bn080916009_v10.fit"),
                      os.path.join(data_dir_lle, "gll_pt_bn080916009_v10.fit"),
                      os.path.join(data_dir_lle, "gll_cspec_bn080916009_v10.rsp"),
                      src_selection,
                       "-100-0,100-200")



nai3 = FermiGBMTTELike('NAI3',
                       os.path.join(data_dir_gbm, "glg_tte_n3_bn080916009_v01.fit.gz"),
                       os.path.join(data_dir_gbm, "glg_cspec_n3_bn080916009_v00.rsp2"),
                       src_selection,
                        "-10-0,120-200",
                       poly_order=2)

nai4 = FermiGBMTTELike('NAI4',
                       os.path.join(data_dir_gbm, "glg_tte_n4_bn080916009_v01.fit.gz"),
                       os.path.join(data_dir_gbm, "glg_cspec_n4_bn080916009_v00.rsp2"),
                       src_selection,
                       "-10-0,120-200",
                       poly_order=-1,
                      verbose=False)

bgo0 = FermiGBMTTELike('BGO0',
                       os.path.join(data_dir_gbm, "glg_tte_b0_bn080916009_v01.fit.gz"),
                       os.path.join(data_dir_gbm, "glg_cspec_b0_bn080916009_v00.rsp2"),
                        src_selection,
                       "-10-0,120-200")



"""
Explanation: Plugin building
Build our plugins from our data. The LLE plugin is an EventList like the GBM TTE plugin and operates in a similar way. 
Simply select a source interval, background interval(s), etc. (all of which can be changed later) and the plugin will prepare the data for you.
As with the TTE plugin (and any OGIPLike plugin), the data can be saved out to PHA files for cross-checking with XSPEC.
End of explanation
"""


lle.display()

"""
Explanation: What's in the plugin?
End of explanation
"""


lle.view_lightcurve(-10,100,dt=.5)

"""
Explanation: Let's look at the lightcurve of LLE file to check out background fit:
End of explanation
"""


nai3.view_lightcurve(-10,100.,.5)

"""
Explanation: And an NaI:
End of explanation
"""



lle.set_active_measurements("50000-100000")
lle.view_count_spectrum()


nai3.set_active_measurements("10.0-30.0", "40.0-900.0")
nai4.set_active_measurements("10.0-30.0", "40.0-900.0")
bgo0.set_active_measurements("250-43000")


"""
Explanation: Energy selection
Like any OGIPLike, we need to select the energies we would like to fit. LLE data below 50 MeV has a lot of dispersion and should not be selected as indicated by the BAD (red) region.
End of explanation
"""



ra = 121.8
dec = -61.3

data_list = DataList(nai3,nai4,bgo0,lle )

band= Band()

GRB = PointSource( trigger_number, ra, dec, spectral_shape=band )

model = Model( GRB )


jl = JointLikelihood( model, data_list, verbose=False )

res = jl.fit()

"""
Explanation: Fitting!
We are now ready for the standard 3ML process. We perform and MLE fit here, but Bayesian fitting is also possible.
End of explanation
"""


sp = SpectralPlotter(jl)

_=sp.plot_model(y_unit='erg2/(cm2 s keV)', num_ene=200,x_max=1E5)

"""
Explanation: Examining our fit
End of explanation
"""


_ = display_ogip_model_counts(jl,min_rate=1,step=False)

res = jl.get_errors()

res = jl.get_contours(band.xp,200,900,20)

res = jl.get_contours(band.xp,250,900,50,band.alpha,-1.2,-0.9,50)

"""
Explanation: We can examine our fit with the data:
End of explanation
"""


