From b6640c4fa29999918c7e1e0508227273e6aed5d0 Mon Sep 17 00:00:00 2001
From: Michael Mayer <michael.mayer@desy.de>
Date: Thu, 2 Jul 2015 14:08:38 +0200
Subject: [PATCH] Added code snippets to retrieve analysis results in python

---
 doc/source/user_manual/getting_started/python.rst | 103 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 scripts/csspec.py                                 |  12 ++++++++++++
 2 files changed, 115 insertions(+), 0 deletions(-)

diff --git a/doc/source/user_manual/getting_started/python.rst b/doc/source/user_manual/getting_started/python.rst
index 117b47b..8b7ce51 100644
--- a/doc/source/user_manual/getting_started/python.rst
+++ b/doc/source/user_manual/getting_started/python.rst
@@ -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
 ~~~~~~~~~~~~~~~~~~~~~~
 
diff --git a/scripts/csspec.py b/scripts/csspec.py
index 796ad9b..bfc5b50 100755
--- a/scripts/csspec.py
+++ b/scripts/csspec.py
@@ -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()
--
libgit2 0.21.4