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..65ad310 100755 --- a/scripts/csspec.py +++ b/scripts/csspec.py @@ -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()