<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Using ctools from Python
In this notebook you will learn how to use the ctools and cscripts from Python instead of typing the commands in the console.
ctools provides two Python modules that allow using all tools and scripts as Python classes. To use ctools from Python you have to import the ctools and cscripts modules into Python. You should also import the gammalib module, as ctools without GammaLib is generally not very useful.
Warning
Step1:  Simulating events
As an example we will simulate an observation of an hour of the Crab nebula.
Step2:  The first line generates an instance of the ctobssim tool as a Python class. User parameters are then set using the [ ] operator. After setting all parameters the execute() method is called to execute the ctobssim tool. On output the events.fits FITS file is created. Until now everything is analogous to running the tool from the command line, but in Python you can easily combine the different blocks into more complex workflows.
Remember that you can consult the manual of each tool to find out how it works and to discover all the input parameters that you can set.
Step3:  In a Jupyter notebook a code line starting with ! is executed in the shell, so you can do the operation above just from the command line.
One of the advantages of using ctools from Python is that you can run a tool using
Step4:  The main difference to the execute() method is that the run() method will not write the output (i.e., the simulated event list) to disk. Why is this useful? Well, after having typed sim.run() the ctobssim class still exists as an object in memory, including all the simulated events. You can always save to disk later using the save() method.
The ctobssim class has an obs() method that returns an observation container that holds the simulated CTA observation with its associated events. To visualise this container, type
Step5:  There is one CTA observation in the container and to visualise that observation type
Step6:  The observation contains a CTA event list that is implement by the GammaLib class GCTAEventList. You can access the event list using the events() method. To visualise the individual events you can iterate over the events using a for loop. This will show the simulated celestial coordinates (RA, DEC), the coordinate in the camera system [DETX, DETY], the energies and the terrestrial times (TT) of all events. Let's peek at the first events of the list
Step7:  We can use this feature to inspect some of the event properties, for example look at their energy spectrum. For this we will use the Python packages matplotlib. If you do not have matplotlib you can use another plotting package of your choice or skip this step.
Step8:  Fitting a model to the observations
We can use the observation in memory to directly run a likelihood fit.
Step9:  This is very compact. Where do we define the model fit to the data? Where are the user parameters? ctlike doesn’t in fact need any parameters as all the relevant information is already contained in the observation container produced by the ctobssim class. Indeed, you constructed the ctlike instance by using the ctobssim observation container as constructor argument.
An observation container, implemented by the GObservations class of GammaLib, is the fundamental brick of any ctools analysis. An observation container can hold more than events, e.g., in this case it also holds the model that was used to generate the events.
Many tools and scripts handle observation containers, and accept them upon construction and return them after running the tool via the obs() method. Passing observation containers between ctools classes is a very convenient and powerful way of building in-memory analysis pipelines. However, this implies that you need some computing ressources when dealing with large observation containers (for example if you want to analyse a few 100 hours of data at once). Also, if the script crashes the information is lost.
To check how the fit went you can inspect the optimiser used by ctlike by typing
Step10:  You see that the fit converged after 2 iterations. Out of 10 parameters in the model 4 have been fitted (the others were kept fixed). To inspect the fit results you can print the model container that can be access using the models() method of the observation container
Step11:  For example, in this way we can fetch the minimum of the optimized function (the negative of the natural logarithm of the likelihood) to compare different model hypotheses.
Step12:  Suppose you want to repeat the fit by optimising also the position of the point source. This is easy from Python, as you can modify the model and fit interactively. Type the following
Step13:  Now we can quantify the improvement of the model by comparing the new value of the optimized function with the previous one.
Step14:  The test statistic (TS) is expected to be distributed as a $\chi^2_n$ with a number of degrees of freedom $n$ equal to the additional number of degrees of freedom of the (second) more complex model with respect to the (first) more parsimonious one, for our case two degrees of freedom. The chance probability that the likelihood improved that much because of pure statistical fluctuations is the integral from TS to infinity of the chi square with two degrees of freedom. In this case the improvement is negligible (i.e., the chance probability is very high), as expected since in the model the source is already at the true position.
Generating a log file
By default, tools and scripts run from Python will not generate a log file. The reason for this is that Python scripts are often used to build ctools analysis pipelines and workflows, and one generally does not want that such a script pollutes the workspace with log files. You can however instruct a ctool or cscript to generate a log file as follows
<ASSISTANT_TASK:>
Python Code:
import gammalib
import ctools
import cscripts
Explanation: Using ctools from Python
In this notebook you will learn how to use the ctools and cscripts from Python instead of typing the commands in the console.
ctools provides two Python modules that allow using all tools and scripts as Python classes. To use ctools from Python you have to import the ctools and cscripts modules into Python. You should also import the gammalib module, as ctools without GammaLib is generally not very useful.
Warning: Always import gammalib before you import ctools and cscripts.
End of explanation
sim = ctools.ctobssim()
sim['inmodel']   = '${CTOOLS}/share/models/crab.xml'
sim['outevents'] = 'events.fits'
sim['caldb']     = 'prod2'
sim['irf']       = 'South_0.5h'
sim['ra']        = 83.5
sim['dec']       = 22.8
sim['rad']       = 5.0
sim['tmin']      = '2020-01-01T00:00:00'
sim['tmax']      = '2020-01-01T01:00:00'
sim['emin']      = 0.03
sim['emax']      = 150.0
sim.execute()
Explanation: Simulating events
As an example we will simulate an observation of an hour of the Crab nebula.
End of explanation
!ctobssim --help
Explanation: The first line generates an instance of the ctobssim tool as a Python class. User parameters are then set using the [ ] operator. After setting all parameters the execute() method is called to execute the ctobssim tool. On output the events.fits FITS file is created. Until now everything is analogous to running the tool from the command line, but in Python you can easily combine the different blocks into more complex workflows.
Remember that you can consult the manual of each tool to find out how it works and to discover all the input parameters that you can set.
End of explanation
sim.run()
Explanation: In a Jupyter notebook a code line starting with ! is executed in the shell, so you can do the operation above just from the command line.
One of the advantages of using ctools from Python is that you can run a tool using
End of explanation
print(sim.obs())
Explanation: The main difference to the execute() method is that the run() method will not write the output (i.e., the simulated event list) to disk. Why is this useful? Well, after having typed sim.run() the ctobssim class still exists as an object in memory, including all the simulated events. You can always save to disk later using the save() method.
The ctobssim class has an obs() method that returns an observation container that holds the simulated CTA observation with its associated events. To visualise this container, type:
End of explanation
print(sim.obs()[0])
Explanation: There is one CTA observation in the container and to visualise that observation type:
End of explanation
events = sim.obs()[0].events()
for event in events[:1]:
    print(event)
Explanation: The observation contains a CTA event list that is implement by the GammaLib class GCTAEventList. You can access the event list using the events() method. To visualise the individual events you can iterate over the events using a for loop. This will show the simulated celestial coordinates (RA, DEC), the coordinate in the camera system [DETX, DETY], the energies and the terrestrial times (TT) of all events. Let's peek at the first events of the list:
End of explanation
import matplotlib.pyplot as plt
%matplotlib inline
#this will visualize plots inline
ax = plt.subplot()
ax.set_yscale('log')
ax.set_xlabel('Log10(Energy [TeV])')
energies = []
for event in events:
    energies.append(event.energy().log10TeV())
    
n, bins, patches = plt.hist(energies)
Explanation: We can use this feature to inspect some of the event properties, for example look at their energy spectrum. For this we will use the Python packages matplotlib. If you do not have matplotlib you can use another plotting package of your choice or skip this step.
End of explanation
like = ctools.ctlike(sim.obs())
like.run()
Explanation: Fitting a model to the observations
We can use the observation in memory to directly run a likelihood fit.
End of explanation
print(like.opt())
Explanation: This is very compact. Where do we define the model fit to the data? Where are the user parameters? ctlike doesn’t in fact need any parameters as all the relevant information is already contained in the observation container produced by the ctobssim class. Indeed, you constructed the ctlike instance by using the ctobssim observation container as constructor argument.
An observation container, implemented by the GObservations class of GammaLib, is the fundamental brick of any ctools analysis. An observation container can hold more than events, e.g., in this case it also holds the model that was used to generate the events.
Many tools and scripts handle observation containers, and accept them upon construction and return them after running the tool via the obs() method. Passing observation containers between ctools classes is a very convenient and powerful way of building in-memory analysis pipelines. However, this implies that you need some computing ressources when dealing with large observation containers (for example if you want to analyse a few 100 hours of data at once). Also, if the script crashes the information is lost.
To check how the fit went you can inspect the optimiser used by ctlike by typing:
End of explanation
print(like.obs().models())
Explanation: You see that the fit converged after 2 iterations. Out of 10 parameters in the model 4 have been fitted (the others were kept fixed). To inspect the fit results you can print the model container that can be access using the models() method of the observation container:
End of explanation
like1 = like.opt().value()
print(like1)
Explanation: For example, in this way we can fetch the minimum of the optimized function (the negative of the natural logarithm of the likelihood) to compare different model hypotheses.
End of explanation
like.obs().models()['Crab']['RA'].free()
like.obs().models()['Crab']['DEC'].free()
like.run()
print(like.obs().models())
Explanation: Suppose you want to repeat the fit by optimising also the position of the point source. This is easy from Python, as you can modify the model and fit interactively. Type the following:
End of explanation
like2 = like.opt().value()
ts    = -2.0 * (like2 - like1)
print(ts)
Explanation: Now we can quantify the improvement of the model by comparing the new value of the optimized function with the previous one.
End of explanation
like.logFileOpen()
like.run()
Explanation: The test statistic (TS) is expected to be distributed as a $\chi^2_n$ with a number of degrees of freedom $n$ equal to the additional number of degrees of freedom of the (second) more complex model with respect to the (first) more parsimonious one, for our case two degrees of freedom. The chance probability that the likelihood improved that much because of pure statistical fluctuations is the integral from TS to infinity of the chi square with two degrees of freedom. In this case the improvement is negligible (i.e., the chance probability is very high), as expected since in the model the source is already at the true position.
Generating a log file
By default, tools and scripts run from Python will not generate a log file. The reason for this is that Python scripts are often used to build ctools analysis pipelines and workflows, and one generally does not want that such a script pollutes the workspace with log files. You can however instruct a ctool or cscript to generate a log file as follows:
End of explanation
<END_TASK>