Add How-to for joint analysis (#1868)

Jürgen Knödlseder
1 parent 3e80f463
Showing 4 changed files with 242 additions and 23 deletions
cscripts/csspec.py
doc/source/users/tutorials/howto/howto_combine_instruments.png
doc/source/users/tutorials/howto/howto_combine_instruments.rst
doc/source/users/tutorials/howto/index.rst
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Generates a spectrum.
 #
-# Copyright (C) 2014-2017 Michael Mayer
+# Copyright (C) 2014-2018 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
@@ -100,26 +100,29 @@ class csspec(ctools.csobservation):
         n_other = self.obs().size() - n_cta
  
         # If spectrum method is not "NODES" then set spectrum method and
-        # script mode according to type of CTA observations
+        # script mode according to type of observations
         if self._method != 'NODES':
-            if n_unbinned == 0 and n_binned != 0 and n_onoff == 0:
-                self._binned_mode = True
-                self._method      = 'SLICE'
-            elif n_unbinned == 0 and n_binned == 0 and n_onoff != 0:
-                self._onoff_mode = True
-                self._method      = 'SLICE'
-            elif n_unbinned == 0 and n_binned != 0 and n_onoff != 0:
-                msg = 'Mix of binned and On/Off CTA observations found in ' \
-                      'observation container. csscript does not support this ' \
-                      'mix.'
-                raise RuntimeError(msg)
-            elif n_unbinned != 0 and (n_binned != 0 or n_onoff != 0):
-                msg = 'Mix of unbinned and binned or On/Off CTA observations ' \
-                      'found in observation container. csscript does not support ' \
-                      'this mix.'
-                raise RuntimeError(msg)
-            elif n_unbinned != 0:
-                self._method = 'SLICE'
+            if n_other > 0:
+                self._method = 'NODES'
+            else:
+                if n_unbinned == 0 and n_binned != 0 and n_onoff == 0:
+                    self._binned_mode = True
+                    self._method      = 'SLICE'
+                elif n_unbinned == 0 and n_binned == 0 and n_onoff != 0:
+                    self._onoff_mode = True
+                    self._method      = 'SLICE'
+                elif n_unbinned == 0 and n_binned != 0 and n_onoff != 0:
+                    msg = 'Mix of binned and On/Off CTA observations found ' \
+                          'in observation container. csscript does not support ' \
+                          'this mix.'
+                    raise RuntimeError(msg)
+                elif n_unbinned != 0 and (n_binned != 0 or n_onoff != 0):
+                    msg = 'Mix of unbinned and binned or On/Off CTA observations ' \
+                          'found in observation container. csscript does not ' \
+                          'support this mix.'
+                    raise RuntimeError(msg)
+                elif n_unbinned != 0:
+                    self._method = 'SLICE'
  
         # Set ebounds
         self._set_ebounds()
@@ -148,10 +151,9 @@ class csspec(ctools.csobservation):
         # If there are a mix of CTA and non-CTA observations and the method
         # is 'SLICE' then log a warning that non-CTA observations will be
         # ignored
+        warning = False
         if n_cta > 0 and n_other > 0 and self._method == 'SLICE':
-            self._log_string(gammalib.TERSE, '\nWARNING: Only CTA observation '
-                             'can be handled with the "SLICE" method, all '
-                             'non-CTA observation will be ignored.')
+            warning = True
  
         # If there are only non-CTA observations and the method is 'SLICE'
         # then stop now
@@ -167,6 +169,12 @@ class csspec(ctools.csobservation):
         # Log selected spectrum method
         self._log_value(gammalib.TERSE, 'Selected spectrum method', self._method)
  
+        # Signal warning
+        if warning:
+            self._log_string(gammalib.TERSE, ' WARNING: Only CTA observation '
+                             'can be handled with the "SLICE" method, all '
+                             'non-CTA observation will be ignored.')
+
         # Return
         return
  
+.. _howto_combine_instrument:
+
+How to jointly analyse data from different instruments?
+-------------------------------------------------------
+
+  .. admonition:: What you will learn
+
+     You will learn how to **jointly analyse data from different instruments**.
+
+In this tutorial you will learn how to jointly analyse data observed with
+Fermi-LAT and simulated data for CTA for the Vela pulsar and pulsar wind
+nebula.
+
+Prepare Fermi-LAT data
+^^^^^^^^^^^^^^^^^^^^^^
+
+The preparation of Fermi-LAT data will follow the one
+`described in this section <howto_fermi_prepare>`_.
+
+
+Prepare CTA data
+^^^^^^^^^^^^^^^^
+
+For CTA the simulations of the Galactic Plane Scan from the first Data Challenge
+will be used, but if you do not have these data at hand you can create your own
+simulated CTA dataset by
+`following the description here <start_simulating>`_.
+
+First, bin the simulated data into a counts cube
+
+.. code-block:: bash
+
+   $ ctbin
+   Input event list or observation definition XML file [events.fits] obs.xml
+   First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63] 128.838
+   Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.51] -45.178
+   Projection method (AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN) [CAR]
+   Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL]
+   Image scale (in degrees/pixel) [0.02]
+   Size of the X axis in pixels [200]
+   Size of the Y axis in pixels [200]
+   Algorithm for defining energy bins (FILE|LIN|LOG) [LOG]
+   Lower energy limit (TeV) [0.1] 0.03
+   Upper energy limit (TeV) [100.0] 150.0
+   Number of energy bins (1-200) [20] 40
+   Output counts cube file [cntcube.fits]
+
+Then compute the response cubes as follows
+
+.. code-block:: bash
+
+   $ ctexpcube
+   Input event list or observation definition XML file [NONE] obs.xml
+   Input counts cube file to extract exposure cube definition [NONE] cntcube.fits
+   Output exposure cube file [expcube.fits]
+
+.. code-block:: bash
+
+   $ ctpsfcube
+   Input event list or observation definition XML file [NONE] obs.xml
+   Input counts cube file to extract PSF cube definition [NONE]
+   First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63] 128.838
+   Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.51] -45.178
+   Projection method (AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN) [CAR]
+   Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL]
+   Image scale (in degrees/pixel) [1.0]
+   Size of the X axis in pixels [10]
+   Size of the Y axis in pixels [10]
+   Algorithm for defining energy bins (FILE|LIN|LOG) [LOG]
+   Lower energy limit (TeV) [0.1] 0.03
+   Upper energy limit (TeV) [100.0] 150.0
+   Number of energy bins (1-200) [20] 40
+   Output PSF cube file [psfcube.fits]
+
+.. code-block:: bash
+
+   $ ctbkgcube
+   Input event list or observation definition XML file [NONE] obs.xml
+   Input counts cube file to extract background cube definition [NONE] cntcube.fits
+   Input model definition XML file [NONE] $CTOOLS/share/models/bkg_irf.xml
+   Output background cube file [bkgcube.fits]
+   Output model definition XML file [NONE] bkgcube.xml
+
+
+Combining the observations
+^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Now you have all the data and hand. You have to create an
+:ref:`observation definition file <glossary_obsdef>`
+to combine the data for an analysis:
+
+.. code-block:: xml
+
+   <?xml version="1.0" standalone="no"?>
+   <observation_list title="observation library">
+     <observation name="Vela" id="000001" instrument="CTA">
+       <parameter name="CountsCube"   file="cntcube.fits"/>
+       <parameter name="ExposureCube" file="expcube.fits"/>
+       <parameter name="PsfCube"      file="psfcube.fits"/>
+       <parameter name="BkgCube"      file="bkgcube.fits"/>
+     </observation>
+     <observation name="Vela" id="000001" instrument="LAT">
+       <parameter name="CountsMap"    file="srcmaps.fits"/>
+       <parameter name="ExposureMap"  file="expmap.fits"/>
+       <parameter name="LiveTimeCube" file="ltcube.fits"/>
+       <parameter name="IRF"          value="P8R2_SOURCE_V6"/>
+     </observation>
+   </observation_list>
+
+The
+:ref:`observation definition file <glossary_obsdef>`
+contains two observations of the Vela pulsar, a first done with CTA and a
+second done with Fermi-LAT. The ``instrument`` attribute distinguishes between
+both instruments.
+
+
+Generate a spectral energy distribution
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Before being able to generate a spectral energy distribution (SED) you have
+to define a
+:ref:`model definition file <glossary_moddef>`
+that will be used to model the events for both observations.
+The model that will be used in this analysis is shown below:
+
+.. code-block:: xml
+
+   <?xml version="1.0" standalone="no"?>
+   <source_library title="source library">
+     <source type="PointSource" name="Vela">
+       <spectrum type="PowerLaw">
+          <parameter name="Prefactor"   scale="1e-16" value="5.7"  min="1e-07" max="1000.0" free="1"/>
+          <parameter name="Index"       scale="-1"    value="2.48" min="0.0"   max="+5.0"   free="1"/>
+          <parameter name="PivotEnergy" scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="0"/>
+       </spectrum>
+       <spatialModel type="PointSource">
+         <parameter name="RA"  scale="1.0" value="128.84" min="-360" max="360" free="1"/>
+         <parameter name="DEC" scale="1.0" value="-45.18" min="-90"  max="90"  free="1"/>
+       </spatialModel>
+     </source>
+     <source type="DiffuseSource" name="Galactic_diffuse" instrument="LAT">
+       <spectrum type="Constant">
+         <parameter name="Normalization" scale="1.0" value="1.0" min="0.1" max="1000.0" free="1"/>
+       </spectrum>
+       <spatialModel type="DiffuseMapCube" file="gll_iem_v06.fits">
+         <parameter name="Normalization" scale="1.0" value="1.0" min="0.1" max="10.0" free="0"/>
+       </spatialModel>
+     </source>
+     <source type="DiffuseSource" name="Extragalactic_diffuse" instrument="LAT">
+       <spectrum type="FileFunction" file="iso_P8R2_SOURCE_V6_v06.txt">
+         <parameter name="Normalization" scale="1.0" value="1.0" min="0.0" max="1000.0" free="0"/>
+       </spectrum>
+       <spatialModel type="DiffuseIsotropic">
+         <parameter name="Value" scale="1.0" value="1.0" min="0.0" max="10.0" free="0"/>
+       </spatialModel>
+     </source>
+     <source name="Background" type="CTACubeBackground" instrument="CTA">
+       <spectrum type="PowerLaw">
+         <parameter name="Prefactor"   scale="1.0"  value="1.0"  min="1e-3" max="1e+3"   free="1"/>
+         <parameter name="Index"       scale="1.0"  value="0.0"  min="-5.0" max="+5.0"   free="1"/>
+         <parameter name="PivotEnergy" scale="1e6"  value="1.0"  min="0.01" max="1000.0" free="0"/>
+       </spectrum>
+     </source>
+   </source_library>
+
+The model contains a point source located at the position of the Vela pulsar
+with a power law spectrum. In addition, it contains two ``DiffuseSource``
+components that are only applied for Fermi-LAT observations, which is indicated
+by their ``instrument="LAT"`` attribute. Both components model the diffuse
+background that prevails at GeV energies. Finally, the model contains a
+``CTACubeBackground`` component that applies to CTA. You may have noticed that
+the point source is the only component that has no ``instrument`` attribute,
+meaning that this component applies to both instruments.
+
+Now you are ready to generate the spectral energy distribution for the
+combined data set. You do this using the :ref:`csspec` script as follows
+
+.. code-block:: bash
+
+   $ csspec
+   Input event list, counts cube, or observation definition XML file [events.fits] obs.xml
+   Input model definition XML file [$CTOOLS/share/models/crab.xml] models.xml
+   Source name [Crab] Vela
+   Spectrum generation method (SLICE|NODES|AUTO) [AUTO]
+   Algorithm for defining energy bins (FILE|LIN|LOG) [LOG]
+   Start value for first energy bin in TeV [0.1] 0.0001
+   Stop value for last energy bin in TeV [100.0] 150.0
+   Number of energy bins (1-200) [20]
+   Output spectrum file [spectrum.fits]
+
+This will generate a logarithmically spaced spectrum composed of 20 energy bins
+comprised within 100 MeV and 150 TeV. The :ref:`csspec` tool is run in the ``AUTO``
+mode, which for different instruments corresponds to the ``NODES`` method.
+The resulting spectral energy distribution is shown below:
+
+.. figure:: howto_combine_instruments.png
+   :width: 600px
+   :align: center
+
+   *Vela spectrum derived using csspec from real Fermi-LAT and simulated CTA data*
+
+.. note::
+   The figure was created using the ``show_spectrum.py`` script that is
+   located in the ctools example folder. The example script requires the
+   `matplotlib <http://matplotlib.org>`_ Python module for display.
+   You may reproduce the plot by typing
+
+   .. code-block:: bash
+
+      $ $CTOOLS/share/examples/python/show_spectrum.py spectrum.fits
@@ -46,3 +46,4 @@ Analysis of data from other instruments
    How to analyse Fermi-LAT data with ctools? <fermi/index>
    How to analyse COMPTEL data with ctools? <comptel/index>
    How to analyse H.E.S.S. data with ctools? <iact/index>
+   howto_combine_instruments