Merge branch '1501-add-documentation-for-python-interface' into integration

Jürgen Knödlseder
2 parents 370ee4fa 10b8ae6a
Showing 2 changed files with 126 additions and 11 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
 ~~~~~~~~~~~~~~~~~~~~~~
@@ -281,19 +281,21 @@ class csspec(ctools.cscript):
         table.card("TELESCOP", "CTA", "Name of Telescope")
         # Create FITS table columns
-        energy      = gammalib.GFitsTableDoubleCol("Energy", self.m_ebounds.size())
+        energy       = gammalib.GFitsTableDoubleCol("Energy", self.m_ebounds.size())
         energy.unit("TeV")
-        energy_low  = gammalib.GFitsTableDoubleCol("ed_Energy", self.m_ebounds.size())
+        energy_low   = gammalib.GFitsTableDoubleCol("ed_Energy", self.m_ebounds.size())
         energy_low.unit("TeV")
         energy_high  = gammalib.GFitsTableDoubleCol("eu_Energy", self.m_ebounds.size())
         energy_high.unit("TeV")
-        flux        = gammalib.GFitsTableDoubleCol("Flux", self.m_ebounds.size())
+        flux         = gammalib.GFitsTableDoubleCol("Flux", self.m_ebounds.size())
         flux.unit("erg/cm2/s")
-        flux_err    = gammalib.GFitsTableDoubleCol("e_Flux", self.m_ebounds.size())
+        flux_err     = gammalib.GFitsTableDoubleCol("e_Flux", self.m_ebounds.size())
         flux_err.unit("erg/cm2/s")
-        TSvalues    = gammalib.GFitsTableDoubleCol("TS", self.m_ebounds.size())
-        ulim_values = gammalib.GFitsTableDoubleCol("UpperLimit", self.m_ebounds.size())
+        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()):
@@ -455,10 +457,11 @@ class csspec(ctools.cscript):
                     self.log("No event in this bin. Bin is skipped\n")
                 # Set all values to 0
-                flux[i]        = 0.0
-                flux_err[i]    = 0.0
-                TSvalues[i]    = 0.0
-                ulim_values[i] = 0.0
+                flux[i]         = 0.0
+                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
@@ -520,7 +531,7 @@ class csspec(ctools.cscript):
                     self.log(", UL = "+str(ulim_values[i])+" [erg/cm2/s]")
                 self.log("\n")
-        # Append filles columns to fits table    
+        # Append filled columns to fits table    
         table.append(energy)
         table.append(energy_low)
         table.append(energy_high)
@@ -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()