Merge branch '1051-create-ctbutterfly' into integration

Jürgen Knödlseder
2 parents d975bd1a 6b06b2a7
Showing 6 changed files with 442 additions and 1 deletions
doc/source/fhelp/ctbutterfly.html
doc/source/tools_reference.rst
doc/source/user_manual/getting_started/butterfly.jpg
doc/source/user_manual/getting_started/quickstart.rst
examples/Makefile.am
examples/show_butterfly.py
+<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<html lang="en"><head>
+<meta content="text/html; charset=ISO-8859-1" http-equiv="content-type"><title>ctbutterfly</title>
+
+<link rel="stylesheet" type="text/css" media="screen,print" href="ctools.css">
+<meta content="Michael Mayer" name="author">
+<meta content="Spectral butterfly (confidence band) computation tools" name="description">
+</head>
+<body>
+<p class="Header">ctbutterfly</p>
+<p class="Body">Butterfly calculation tool.</p>
+<br>
+<p class="Title">Synopsis</p>
+<p style="padding-top: 0pt; text-align: left;" class="Body">Calculates
+the confidence band of a spectral model, taking into account the
+covariance matrix of a likelihood fit.
+</p>
+<br>
+<p class="Title">Links</p>
+<ul class="Body">
+<li><a href="ftools_usage.html">Basic FTOOL
+Parameter Usage</a></li>
+</ul>
+<p style="padding-top: 0pt; text-align: left;" class="Body"></p>
+<br>
+<p class="Title">General parameters</p>
+<table id="table1" class="Body" border="0" cellpadding="0" cellspacing="0" width="600">
+<tbody>
+<tr>
+<td width="50">&nbsp;</td>
+<td colspan="2"><span class="irafpar"><strong>infile
+[file]</strong></span></td>
+</tr>
+<tr>
+<td width="50">&nbsp;</td>
+<td width="25">&nbsp;</td>
+<td width="525">XML or fits file that defines the observation</td>
+</tr>
+<tr>
+<td width="50">&nbsp;</td>
+<td>&nbsp;</td>
+<td width="525">&nbsp;</td>
+</tr>
+<tr>
+<td width="50">&nbsp;</td>
+<td colspan="2"><span class="irafpar"><strong>srcmdl
+[file]</strong></span></td>
+</tr>
+<tr>
+<td width="50">&nbsp;</td>
+<td width="25">&nbsp;</td>
+<td width="525">Model XML file containing the source and background definitions</td>
+</tr>
+<tr>
+<td width="50">&nbsp;</td>
+<td>&nbsp;</td>
+<td width="525">&nbsp;</td>
+</tr>
+<td width="50">&nbsp;</td>
+<td colspan="2"><span class="irafpar"><strong>srcname
+    [file]</strong></span></td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td width="25">&nbsp;</td>
+    <td width="525">Name of model for which butterfly should be computed</td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td>&nbsp;</td>
+    <td width="525">&nbsp;</td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td colspan="2"><span class="irafpar"><strong>caldb
+        [string]</strong></span></td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td width="25">&nbsp;</td>
+    <td width="525">Calibration database repository.</td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td>&nbsp;</td>
+    <td width="525">&nbsp;</td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td colspan="2"><span class="irafpar"><strong>irf
+        [string]</strong></span></td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td width="25">&nbsp;</td>
+    <td width="525">Instrumental response function that
+        should be used for event generation.</td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td>&nbsp;</td>
+    <td width="525">&nbsp;</td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td colspan="2"><span class="irafpar"><strong>outfile
+        [file]</strong></span></td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td width="25">&nbsp;</td>
+    <td width="525">Output butterfly file name</td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td>&nbsp;</td>
+    <td width="525">&nbsp;</td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td colspan="2"><span class="irafpar"><strong>(matrix = "NONE")
+        [file]</strong></span></td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td width="25">&nbsp;</td>
+    <td width="525">Input covariance matrix file (not used)</td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td>&nbsp;</td>
+    <td width="525">&nbsp;</td>
+</tr>
+<td width="50">&nbsp;</td>
+<td colspan="2"><span class="irafpar"><strong>emin
+    [double]</strong></span></td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td width="25">&nbsp;</td>
+    <td width="525">Minimum energy of butterfly (in
+        TeV).</td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td>&nbsp;</td>
+    <td width="525">&nbsp;</td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td colspan="2"><span class="irafpar"><strong>emax
+        [double]</strong></span></td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td width="25">&nbsp;</td>
+    <td width="525">Maximum energy of butterfly (in
+        TeV).</td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td>&nbsp;</td>
+    <td width="525">&nbsp;</td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td colspan="2"><span class="irafpar"><strong>(ebinalg
+        = "LOG") &lt;FILE|LIN|LOG&gt;
+        [string]</strong></span></td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td width="25">&nbsp;</td>
+    <td width="525">Butterfly binning algorithm</td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td>&nbsp;</td>
+    <td width="525">&nbsp;</td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td colspan="2"><span class="irafpar"><strong>(enumbins
+        = 100)
+        [string]</strong></span></td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td width="25">&nbsp;</td>
+    <td width="525">Number of energy bins of butterfly</td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td>&nbsp;</td>
+    <td width="525">&nbsp;</td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td colspan="2"><span class="irafpar"><strong>(ebinfile
+        = "NONE")
+        [string]</strong></span></td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td width="25">&nbsp;</td>
+    <td width="525">File for energy binning algorithm</td>
+</tr>
+<tr>
+    <td width="50">&nbsp;</td>
+    <td>&nbsp;</td>
+    <td width="525">&nbsp;</td>
+</tr>
+</tbody>
+</table>
+<p class="Title">Standard parameters</p>
+<table id="table1" class="Body" border="0" cellpadding="0" cellspacing="0" width="600">
+<tbody>
+<tr>
+<td width="50">&nbsp;</td>
+<td colspan="2"><span class="irafpar"><strong>(chatter
+= 2) [integer]</strong></span></td>
+</tr>
+<tr>
+<td width="50">&nbsp;</td>
+<td width="25">&nbsp;</td>
+<td width="525">Verbosity of the executable:<br>
+<code>chatter = 0</code>: no information will be
+logged<br>
+<code>chatter = 1</code>: only errors will be logged<br>
+<code>chatter = 2</code>: errors and actions will be
+logged<br>
+<code>chatter = 3</code>: report about the task
+execution<br>
+<code>chatter = 4</code>: detailed report about the
+task execution</td>
+</tr>
+<tr>
+<td width="50">&nbsp;</td>
+<td>&nbsp;</td>
+<td width="525">&nbsp;</td>
+</tr>
+<tr>
+<td width="50">&nbsp;</td>
+<td colspan="2"><span class="irafpar"><strong>(clobber
+= yes) [boolean]</strong></span></td>
+</tr>
+<tr>
+<td width="50">&nbsp;</td>
+<td width="25">&nbsp;</td>
+<td width="525">Specifies whether an existing output
+event list should be overwritten.</td>
+</tr>
+<tr>
+<td width="50">&nbsp;</td>
+<td>&nbsp;</td>
+<td width="525">&nbsp;</td>
+</tr>
+<tr>
+<td width="50">&nbsp;</td>
+<td colspan="2"><span class="irafpar"><strong>(debug
+= no) [boolean]</strong></span></td>
+</tr>
+<tr>
+<td width="50">&nbsp;</td>
+<td width="25">&nbsp;</td>
+<td width="525">Enables debug mode. In debug mode
+the executable will dump any log file output in the console.</td>
+</tr>
+<tr>
+<td width="50">&nbsp;</td>
+<td>&nbsp;</td>
+<td width="525">&nbsp;</td>
+</tr>
+<tr>
+<td width="50">&nbsp;</td>
+<td colspan="2"><span class="irafpar"><strong>(mode
+= ql) [string]</strong></span></td>
+</tr>
+<tr>
+<td width="50">&nbsp;</td>
+<td width="25">&nbsp;</td>
+<td width="525">Mode of automatic parameters.</td>
+</tr>
+<tr>
+<td width="50">&nbsp;</td>
+<td>&nbsp;</td>
+<td width="525">&nbsp;</td>
+</tr>
+</tbody>
+</table>
+<br>
+<p class="Title">Related FTOOLS</p>
+<br>
+<hr size="1">
+<p class="Body">Last updated by: Michael
+Mayer 30 October
+2014
+</p>
+<br>
+</body></html>
 \ No newline at end of file
@@ -14,6 +14,7 @@ ctools
 * `ctobssim <fhelp/ctobssim.html>`__ --- Observation simulator.
 * `ctselect <fhelp/ctselect.html>`__ --- Performs sub-selections of event data which are needed for unbinned maximum likelihood analysis.
 * `ctskymap <fhelp/ctskymap.html>`__ --- Computes sky map from CTA data.
+* `ctbutterfly <fhelp/ctbutterfly.html>`__ --- Computes a butterfly from a spectral model.
 cscripts
 --------
@@ -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,64 @@ 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. It will write the 
+butterfly information into an ASCII file.
+The following example shows how to compute such a butterfly from the command
+line.
+
+.. 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: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 and 
+that you will have an ASCII file named ``butterfly.txt'' on disk you can
+visualise the butterfly using the script ``show_butterfly.py'' that is 
+in the ctools ``example'' folder. You will need matplotlib on your system
+to make this work. To launch the script, type:
+
+.. code-block:: bash
+	python $CTOOLS/examples/python/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
@@ -7,7 +7,8 @@ scripts = README \
           make_ts_distributions.py \
           make_unbinned_analysis.py \
           show_pull_histogram.py \
-          show_ts_distribution.py
+          show_ts_distribution.py \
+          show_butterfly
 # Set directories
 scriptsdir   = $(datadir)/examples/python
+#! /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 <http://www.gnu.org/licenses/>.
+#
+# ==========================================================================
+
+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