From df8dd6d936cca8ebe00aeb35a46b4cc32cc8de86 Mon Sep 17 00:00:00 2001 From: Michael Mayer Date: Thu, 30 Oct 2014 19:13:54 +0100 Subject: [PATCH] Added documentation and script to plot butterfly --- doc/source/fhelp/ctbutterfly.html | 300 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/tools_reference.rst | 1 + doc/source/user_manual/getting_started/butterfly.jpg | Bin 0 -> 60362 bytes doc/source/user_manual/getting_started/quickstart.rst | 129 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ examples/show_butterfly.py | 76 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 506 insertions(+), 0 deletions(-) create mode 100644 doc/source/fhelp/ctbutterfly.html create mode 100644 doc/source/user_manual/getting_started/butterfly.jpg create mode 100644 examples/show_butterfly.py diff --git a/doc/source/fhelp/ctbutterfly.html b/doc/source/fhelp/ctbutterfly.html new file mode 100644 index 0000000..f170cd6 --- /dev/null +++ b/doc/source/fhelp/ctbutterfly.html @@ -0,0 +1,300 @@ + + +ctbutterfly + + + + + + +

ctbutterfly

+

Butterfly calculation tool.

+
+

Synopsis

+

Calculates +the confidence band of a spectral model, taking into account the +covariance matrix of a likelihood fit. +

+
+

Links

+ +

+
+

General parameters

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
 infile +[file]
  XML or fits file that defines the observation
   
 srcmdl +[file]
  Model XML file containing the source and background definitions
   
 srcname + [file]
  Name of model for which butterfly should be computed
   
 caldb + [string]
  Calibration database repository.
   
 irf + [string]
  Instrumental response function that + should be used for event generation.
   
 outfile + [file]
  Output butterfly file name
   
 (matrix = "NONE") + [file]
  Input covariance matrix file (not used)
   
 emin + [double]
  Minimum energy of butterfly (in + TeV).
   
 emax + [double]
  Maximum energy of butterfly (in + TeV).
   
 (ebinalg + = "LOG") <FILE|LIN|LOG> + [string]
  Butterfly binning algorithm
   
 (enumbins + = 100) + [string]
  Number of energy bins of butterfly
   
 (ebinfile + = "NONE") + [string]
  File for energy binning algorithm
   
+

Standard parameters

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
 (chatter += 2) [integer]
  Verbosity of the executable:
+chatter = 0: no information will be +logged
+chatter = 1: only errors will be logged
+chatter = 2: errors and actions will be +logged
+chatter = 3: report about the task +execution
+chatter = 4: detailed report about the +task execution
   
 (clobber += yes) [boolean]
  Specifies whether an existing output +event list should be overwritten.
   
 (debug += no) [boolean]
  Enables debug mode. In debug mode +the executable will dump any log file output in the console.
   
 (mode += ql) [string]
  Mode of automatic parameters.
   
+
+

Related FTOOLS

+
+
+

Last updated by: Michael +Mayer 30 October +2014 +

+
+ \ No newline at end of file diff --git a/doc/source/tools_reference.rst b/doc/source/tools_reference.rst index 6a9903b..082a6e8 100644 --- a/doc/source/tools_reference.rst +++ b/doc/source/tools_reference.rst @@ -14,6 +14,7 @@ ctools * `ctobssim `__ --- Observation simulator. * `ctselect `__ --- Performs sub-selections of event data which are needed for unbinned maximum likelihood analysis. * `ctskymap `__ --- Computes sky map from CTA data. +* `ctbutterfly `__ --- Computes a butterfly from a spectral model. cscripts -------- diff --git a/doc/source/user_manual/getting_started/butterfly.jpg b/doc/source/user_manual/getting_started/butterfly.jpg new file mode 100644 index 0000000..814edc1 Binary files /dev/null and b/doc/source/user_manual/getting_started/butterfly.jpg differ diff --git a/doc/source/user_manual/getting_started/quickstart.rst b/doc/source/user_manual/getting_started/quickstart.rst index edc0d2e..824d82e 100644 --- a/doc/source/user_manual/getting_started/quickstart.rst +++ b/doc/source/user_manual/getting_started/quickstart.rst @@ -19,6 +19,8 @@ The example consists of the following analysis steps: - :ref:`sec_unbinned_cta` +- :ref:`sec_butterfly_cta` + .. _sec_simulating_cta: @@ -638,3 +640,130 @@ the same events for both analyses. 2014-01-09T21:47:54: +==============+ 2014-01-09T21:47:54: 2014-01-09T21:47:54: Application "ctlike" terminated after 6 wall clock seconds, consuming 0.163384 seconds of CPU time. + + .. _sec_butterfly_cta: + +Calculate and visualise butterfly +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +To visualise the analysis results retrieved above, one can calculate the +confidence band of the spectral fit. The tool ctbutterfly takes the optimised +source model as input. It takes the covariance matrix from the fit to conduct a +Gaussian error propagation for each energy value. +The following example shows how to compute such a butterfly from the command line +and subsequently shows how to visualise the spectral fit using matplotlib. + +.. code-block:: bash + + $ ctbutterfly + Input event list, cube or observation definition file [$CTOOLS/test/data/crab_events.fits.gz] + Calibration database [share/caldb/data/cta/dummy//cta_dummy_irf] share/caldb/data/cta/dummy/bcf/ + Instrument response function [cta_dummy_irf] + Source model [$CTOOLS/test/data/crab.xml] + Source of interest [CrabNebula] Crab + Start value for first energy bin in TeV [0.1] + Stop value for last energy bin in TeV [100.] + Output ascii file [butterfly.txt] + +Below some lines of the ``ctbutterfly.log``: + +.. code-block:: xml + + 2014-10-30T17:51:24: ******************************************************************************** + 2014-10-30T17:51:24: * ctbutterfly * + 2014-10-30T17:51:24: * ---------------------------------------------------------------------------- * + 2014-10-30T17:51:24: * Version: 00-01-00 * + 2014-10-30T17:51:24: ******************************************************************************** + 2014-10-30T17:51:39: +============+ + 2014-10-30T17:51:39: | Parameters | + 2014-10-30T17:51:39: +============+ + 2014-10-30T17:51:39: infile ....................: $CTOOLS/test/data/crab_events.fits.gz + 2014-10-30T17:51:39: srcmdl ....................: $CTOOLS/test/data/crab.xml + 2014-10-30T17:51:39: srcname ...................: Crab + 2014-10-30T17:51:39: caldb .....................: share/caldb/data/cta/dummy/bcf/ + 2014-10-30T17:51:39: irf .......................: cta_dummy_irf + 2014-10-30T17:51:39: outfile ...................: butterfly.txt + 2014-10-30T17:51:39: matrix ....................: NONE + 2014-10-30T17:51:39: ebinalg ...................: LOG + 2014-10-30T17:51:39: emin ......................: 0.1 + 2014-10-30T17:51:39: emax ......................: 100. + 2014-10-30T17:51:39: enumbins ..................: 100 + 2014-10-30T17:51:39: ebinfile ..................: NONE + 2014-10-30T17:51:39: chatter ...................: 2 + 2014-10-30T17:51:39: clobber ...................: yes + 2014-10-30T17:51:39: debug .....................: no + 2014-10-30T17:51:39: mode ......................: ql + 2014-10-30T17:51:39: logfile ...................: ctbutterfly.log + 2014-10-30T17:51:39: + 2014-10-30T17:51:39: +=============+ + 2014-10-30T17:51:39: | Observation | + 2014-10-30T17:51:39: +=============+ + 2014-10-30T17:51:39: === GObservations === + 2014-10-30T17:51:39: Number of observations ....: 1 + 2014-10-30T17:51:39: Number of predicted events : 0 + 2014-10-30T17:51:39: === GCTAObservation === + 2014-10-30T17:51:39: Name ......................: + 2014-10-30T17:51:39: Identifier ................: + 2014-10-30T17:51:39: Instrument ................: CTA + 2014-10-30T17:51:39: Event file ................: test/data/crab_events.fits.gz + 2014-10-30T17:51:39: Event type ................: EventList + 2014-10-30T17:51:39: Statistics ................: Poisson + 2014-10-30T17:51:39: Ontime ....................: 1800 s + 2014-10-30T17:51:39: Livetime ..................: 1710 s + 2014-10-30T17:51:39: Deadtime correction .......: 0.95 + 2014-10-30T17:51:39: === GCTAPointing === + 2014-10-30T17:51:39: Pointing direction ........: (RA,Dec)=(83.63,22.01) + 2014-10-30T17:51:39: === GCTAResponseIrf === + 2014-10-30T17:51:39: Response name .............: cta_dummy_irf + 2014-10-30T17:51:39: Energy dispersion .........: Not used + 2014-10-30T17:51:39: === GCaldb === + 2014-10-30T17:51:39: Database root .............: share/caldb/data/cta/dummy/bcf/ + 2014-10-30T17:51:39: === GCTAAeffPerfTable === + 2014-10-30T17:51:39: Filename ..................: share/caldb/data/cta/dummy/bcf//cta_dummy_irf.dat + 2014-10-30T17:51:39: Number of energy bins .....: 20 + 2014-10-30T17:51:39: Log10(Energy) range .......: 0.0199526 - 125.893 TeV + 2014-10-30T17:51:39: Offset angle dependence ...: Fixed sigma=3 + 2014-10-30T17:51:39: === GCTAPsfPerfTable === + 2014-10-30T17:51:39: Filename ..................: share/caldb/data/cta/dummy/bcf//cta_dummy_irf.dat + 2014-10-30T17:51:39: Number of energy bins .....: 20 + 2014-10-30T17:51:39: Log10(Energy) range .......: 0.0199526 - 125.893 TeV + 2014-10-30T17:51:39: === GCTAEdispPerfTable === + 2014-10-30T17:51:39: Filename ..................: share/caldb/data/cta/dummy/bcf//cta_dummy_irf.dat + 2014-10-30T17:51:39: Number of energy bins .....: 20 + 2014-10-30T17:51:39: Log10(Energy) range .......: 0.0199526 - 125.893 TeV + 2014-10-30T17:51:39: === GCTABackgroundPerfTable === + 2014-10-30T17:51:39: Filename ..................: share/caldb/data/cta/dummy/bcf//cta_dummy_irf.dat + 2014-10-30T17:51:39: Number of energy bins .....: 20 + 2014-10-30T17:51:39: Log10(Energy) range .......: 0.0199526 - 125.893 TeV + 2014-10-30T17:51:39: Offset angle dependence ...: Fixed sigma=3 + 2014-10-30T17:51:39: === GCTAEventList === + 2014-10-30T17:51:39: Number of events ..........: 6141 + 2014-10-30T17:51:39: Time interval .............: 51544.5 - 51544.5 days + 2014-10-30T17:51:39: + 2014-10-30T17:51:39: +===========================+ + 2014-10-30T17:51:39: | Compute covariance matrix | + 2014-10-30T17:51:39: +===========================+ + 2014-10-30T17:51:39: + 2014-10-30T17:51:39: +====================+ + 2014-10-30T17:51:39: | Generate butterfly | + 2014-10-30T17:51:39: +====================+ + 2014-10-30T17:51:39: + 2014-10-30T17:51:39: +========================+ + 2014-10-30T17:51:39: | Save Butterfly to file | + 2014-10-30T17:51:39: +========================+ + 2014-10-30T17:51:39: + 2014-10-30T17:51:39: Application "ctbutterfly" terminated after 15 wall clock seconds, consuming 0.051253 seconds of CPU time. + +Now that you have computed the confidence band of the spectral fit, you can visualise +file with an example script. The script is located in the "examples" folder. Make sure you +have installed matplotlib on your system. Then run the following script: + +.. code-block:: bash + python $CTOOLS/examples/show_butterfly.py butterfly.txt + +This will result in a canvas which should look like the following: +.. figure:: butterfly.jpg + :width: 60% + + *Confidence band of the fit* + \ No newline at end of file diff --git a/examples/show_butterfly.py b/examples/show_butterfly.py new file mode 100644 index 0000000..1f522a7 --- /dev/null +++ b/examples/show_butterfly.py @@ -0,0 +1,76 @@ +#! /usr/bin/env python +# ========================================================================== +# This script shows how to plot a butterfly created with ctbutterfly +# +# Copyright (C) 2014 Michael Mayer +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# ========================================================================== + +import gammalib +import sys +try: + import matplotlib.pyplot as plt +except: + sys.exit("This script needs matplotlib") + +# =============== # +# Script entry # +# =============== # +if __name__ == "__main__": + + # check for given butterfly file + if not len(sys.argv) == 2: + sys.exit("Usage: show_butterfly.py butterfly.txt") + + # Read given butterfly file + filename = sys.argv[1] + csv = gammalib.GCsv(filename) + + # initialise arrays to be filled + butterfly_x = [] + butterfly_y = [] + + line_x = [] + line_y = [] + + # loop over rows of the file + nrows = csv.nrows() + for row in range(nrows): + # Compute upper edge of confidence band + butterfly_x.append(csv.real(row,0)) + butterfly_y.append(csv.real(row,1)+csv.real(row,2)) + + # Set line values + line_x.append(csv.real(row,0)) + line_y.append(csv.real(row,1)) + + # Loop over the rows backwards to compute the lower edge + # of the confidence band + for row in range(nrows): + index = nrows - 1 - row + + butterfly_x.append(csv.real(index,0)) + butterfly_y.append(csv.real(index,1)-csv.real(index,2)) + + # plot the butterfly and spectral line + plt.figure() + plt.loglog() + plt.grid() + plt.plot(line_x,line_y,color='black',ls='-') + plt.fill(butterfly_x,butterfly_y,color='green',alpha=0.5) + plt.xlabel("Energy [MeV]") + plt.ylabel(r"dN/dE [MeV$^{-1}$ s$^{-1}$ cm$^{-2}$]") + plt.show() \ No newline at end of file -- libgit2 0.21.4