Added code snippets to retrieve analysis results in python

Michael Mayer
1 parent 907508c6
Showing 2 changed files with 115 additions and 0 deletions
doc/source/user_manual/getting_started/python.rst
scripts/csspec.py
@@ -354,6 +354,109 @@ The ``obsutils.sim()`` function then simulates the event data, the
    functions will be added in the future.
+Access common analysis results
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+Here are some examples how to access your analysis results in python. In the following it is
+assumed that you have a ctlike object called 'like' which is setup and runs:
+
+.. code-block:: python
+   
+   import gammalib
+   import ctools
+   like = ctools.ctlike(sim.obs())
+   like.run()
+   ...
+   
+* Best-fit parameters:
+
+  The following command shows how to access the fit parameters and errors:
+  
+  .. code-block:: python
+
+   obs = like.obs() # Get observations object
+   obs.models()["Crab"]["Prefactor"].value() # This returns the actual fit value
+   obs.models()["Crab"]["Prefactor"].error() # This returns the actual fit error
+   ...
+
+* Open XML model file from python and print models on screen:
+
+  .. code-block:: python
+  
+    models = gammalib.GModels("$CTOOLS/share/models/crab.xml")
+    print models
+    ...
+  
+* Likelihood value
+ 
+ .. code-block:: python
+  
+    like.opt().value() # Returns likelihood value
+    ...    
+    
+* Covariance Matrix
+
+ .. code-block:: python
+   
+   covariance = like.opt().curvature().invert() # Return GMatrix object 
+   print covariance
+   ...
+   
+Modify and work on XML models
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+* Fix a certain parameter
+	
+ .. code-block:: python
+
+   models["Crab"]["Prefactor"].fix()
+   ...
+     
+* Release a certain parameter
+
+ .. code-block:: python
+
+   models["Crab"]["Prefactor"].free()
+   ...
+     
+* Check if parameter is free
+
+ .. code-block:: python
+
+   models["Crab"]["Prefactor"].is_free()
+   ...
+     
+* Set parameter value
+
+ .. code-block:: python
+
+   models["Crab"]["Index"].value(-2.5)
+   ...
+     
+* Set parameter boundaries
+
+ .. code-block:: python 
+   
+   models["Crab"]["Index"].min(-5.0)
+   models["Crab"]["Index"].max(-1.0)
+   ...
+     
+* Loop over models and parameters (e.g. fix all parameters)
+
+ .. code-block:: python 
+ 
+   for model in models: # Loop over models
+   	for par in model: # Loop over parameters
+   		par.fix() # fix parameter
+   ...
+     
+* Save modified XML model
+ 
+ .. code-block:: python
+   
+   models.save("myXMLmodel.xml")
+   ...
+
+
 Beyond the first steps
 ~~~~~~~~~~~~~~~~~~~~~~
@@ -294,6 +294,8 @@ class csspec(ctools.cscript):
         TSvalues    = gammalib.GFitsTableDoubleCol("TS", self.m_ebounds.size())
         ulim_values = gammalib.GFitsTableDoubleCol("UpperLimit", self.m_ebounds.size())
         ulim_values.unit("erg/cm2/s")
+        Npred_values = gammalib.GFitsTableDoubleCol("Npred", self.m_ebounds.size())
+        
         # Loop over energy bins
         for i in range(self.m_ebounds.size()):
@@ -459,6 +461,7 @@ class csspec(ctools.cscript):
                 flux_err[i]    = 0.0
                 TSvalues[i]    = 0.0
                 ulim_values[i] = 0.0
+                Npred_values[i] = 0.0
                 continue
             # Get results
@@ -491,6 +494,11 @@ class csspec(ctools.cscript):
             TS = -1.0
             if self.m_calc_ts:
                 TS = source.ts() 
+            
+            # Compute Npred value
+            Npred = 0.0
+            for observation in like.obs():
+                Npred += observation.npred(source)  
             # Get differential flux    
             fitted_flux = source.spectral().eval(elogmean,gammalib.GTime())
@@ -503,6 +511,9 @@ class csspec(ctools.cscript):
             # Set values for storage
             TSvalues[i] = TS
+            # Set npred values 
+            Npred_values[i] = Npred
+            
             # Convert fluxes to nuFnu
             flux[i]     = fitted_flux * elogmean2 * gammalib.MeV2erg
             flux_err[i] = e_flux      * elogmean2 * gammalib.MeV2erg
@@ -528,6 +539,7 @@ class csspec(ctools.cscript):
         table.append(flux_err)
         table.append(TSvalues)
         table.append(ulim_values)
+        table.append(Npred_values)
         # Create the FITS file now
         self.fits = gammalib.GFits()