Update 1DC tutorials

Jürgen Knödlseder
1 parent 593f8844
Showing 26 changed files with 458 additions and 180 deletions
doc/source/users/tutorials/1dc/first_detect.rst
doc/source/users/tutorials/1dc/first_fitting.rst
doc/source/users/tutorials/1dc/first_improving.rst
doc/source/users/tutorials/1dc/first_select_events.rst
doc/source/users/tutorials/1dc/first_select_obs.rst
doc/source/users/tutorials/1dc/first_skymap.rst
doc/source/users/tutorials/1dc/first_stacked.rst
doc/source/users/tutorials/1dc/first_unbinned.rst
doc/source/users/tutorials/1dc/howto.rst
doc/source/users/tutorials/1dc/howto_display.rst
doc/source/users/tutorials/1dc/howto_display_phases.png
doc/source/users/tutorials/1dc/howto_edisp.rst
doc/source/users/tutorials/1dc/howto_exclude.rst
doc/source/users/tutorials/1dc/howto_extent.rst
doc/source/users/tutorials/1dc/howto_filtered_exclude.png
doc/source/users/tutorials/1dc/howto_lightcurve.rst
doc/source/users/tutorials/1dc/howto_phasecurve.rst
doc/source/users/tutorials/1dc/howto_phasecurve_events.png
doc/source/users/tutorials/1dc/howto_phasecurve_mcid.png
doc/source/users/tutorials/1dc/howto_phasecurve_phases.png
@@ -3,15 +3,15 @@
 Detecting candidate sources
 ---------------------------
-.. admonition:: You will learn ...
+  .. admonition:: What you will learn
-   ... how to **detect candidate sources** in a sky map.
+     You will learn how to **detect candidate sources** in a sky map.
-   This will provide you with an initial
-   :ref:`model definition XML file <glossary_moddef>`
-   that can be directly used for maximum likelihood fitting. In subsequent
-   steps you will iteratively improve this model to obtain a reliable
-   description of the observed gamma-ray sky.
+     This will provide you with an initial
+     :ref:`model definition file <glossary_moddef>`
+     in XML format that can be directly used for maximum likelihood fitting.
+     In subsequent steps you will iteratively improve this model to obtain a
+     reliable description of the observed gamma-ray sky.
 As next step you need to identify the location of candidate sources in the
 sky map.
@@ -3,22 +3,20 @@
 Fitting the model components to the counts cube
 -----------------------------------------------
-.. admonition:: You will learn ...
-
-   ... how to **fit a parametric model** to the counts cube using an
-   maximum liklihood algorithm.
-
-   You will also learn how to display the fit results in form of
-   `ds9 <http://ds9.si.edu>`_
-   region files, butterfly diagrams and spectral points.
-
-   Please note that display tools are **not** part of ctools, yet some
-   scripts for result display that use the ``matplotlib`` Python module can be
-   found in the
-   ``$CTOOLS/share/examples/python``
-   folder. For example you may use ``show_butterfly.py`` to display a butterfly
-   diagram or ``show_spectrum.py`` to display spectral points. Do not hesitate to
-   copy and adjust these scripts as needed, they are examples.
+  .. admonition:: What you will learn
+
+     You will learn how to **fit a parametric model** to the counts cube using
+     a maximum likelihood algorithm.
+
+     You will also learn how to display the fit results in form of
+     `ds9 <http://ds9.si.edu>`_
+     region files, butterfly diagrams and spectral points.
+
+     Please note that display tools are **not** part of ctools, yet some
+     scripts for result display that use the ``matplotlib`` Python module can be
+     found in the
+     ``$CTOOLS/share/examples/python``
+     folder (see :ref:`how to display results<1dc_howto_display>`).
 Now you are ready to fit the model to the counts cube and to determine its
 maximum likelihood parameters.
@@ -3,10 +3,10 @@
 Iteratively improving the model
 -------------------------------
-.. admonition:: You will learn ...
+  .. admonition:: What you will learn
-   ... how to **iteratively improve** your source model by inspecting the
-   fit residuals and by adjusting the model components as needed.
+     You will learn how to **iteratively improve** your source model by
+     inspecting the fit residuals and by adjusting the model components as needed.
 After doing the model fit you should investigate the residuals to verify that
 the model components properly describe the observed event distribution.
@@ -3,14 +3,15 @@
 Selecting the relevant events
 -----------------------------
-.. admonition:: You will learn ...
+  .. admonition:: What you will learn
-   ... how to **select a subset of events** from the selected observations for
-   your analysis.
+     You will learn how to **select a subset of events** from the selected
+     observations for your analysis.
-   The event selection is useful to restrict for example the energy range of
-   the events or the time interval covered, but also to select only events
-   within a given radial acceptance region (also called region of interest).
+     The event selection is necessary if you want to restrict for example the
+     energy range of the events or the time interval that should be covered,
+     but also if you want to select only events within a given radial
+     acceptance region (also called region of interest).
 As next step, you have to select the relevant events from the selected
 observations. You do this with the :ref:`ctselect` tool by typing:
@@ -3,13 +3,13 @@
 Selecting the relevant observations
 -----------------------------------
-.. admonition:: You will learn ...
+  .. admonition:: What you will learn
-   ... how to **select a subset of observations** from the Data Challenge data
-   for your analysis.
+     You will learn how to **select a subset of observations** from the Data
+     Challenge data for your analysis.
-   We recommended to analyse the data region by region since the analysis of
-   the Data Challenge data in a single shot will be very very time consuming.
+     We recommended to analyse the data region by region since the analysis of
+     the Data Challenge data in a single shot will be very very time consuming.
 Let's assume that you want to analyse the central region of our Galaxy using
 the data obtained during the Galactic Plane Survey.
@@ -3,15 +3,15 @@
 Generating a sky map from the events
 ------------------------------------
-.. admonition:: You will learn ...
+  .. admonition:: What you will learn
-   ... how to **generate a sky map** from the selected events.
+     You will learn how to **generate a sky map** from the selected events.
-   A background model based on the background template included in the
-   :ref:`instrument response functions <glossary_irf>`
-   may be subtracted from the sky map to display only celestial excess events,
-   but be aware that for the moment the background subtraction is computationally
-   intensive.
+     A background model based on the background template included in the
+     :ref:`instrument response functions <glossary_irf>`
+     may be subtracted from the sky map to display only celestial excess events,
+     but be aware that for the moment the background subtraction is
+     computationally intensive.
 After having selected the observations and the events you can begin your
 analysis.
@@ -3,18 +3,18 @@
 Stacking the data
 -----------------
-.. admonition:: You will learn ...
-
-   ... how to **prepare the data for analysis** by stacking all events into
-   a counts cube and by computing the effective
-   :ref:`instrument response function <glossary_irf>`
-   for this counts cube.
-
-   Please note that for the moment we do not consider the energy dispersion
-   since for many cases the effect of the energy dispersion is negligible,
-   and taking energy dispersion into account is computationally intensive. If
-   you want to learn how to take the energy dispersion into account please
-   read :ref:`this section<1dc_howto_edisp>`.
+  .. admonition:: What you will learn
+
+     You will learn how to **prepare the data for analysis** by stacking all
+     events into a counts cube and by computing the effective
+     :ref:`instrument response function <glossary_irf>`
+     for this counts cube.
+
+     Please note that for the moment we do not consider the energy dispersion
+     since for many cases the effect of the energy dispersion is negligible,
+     and taking energy dispersion into account is computationally intensive. If
+     you want to learn how to take the energy dispersion into account please
+     read :ref:`this section<1dc_howto_edisp>`.
 To analyse the selected observations, we recommend to stack the events into
 a counts cube. You do this using the :ref:`ctbin` tool:
@@ -3,15 +3,15 @@
 Fitting the model components directly to the events
 ---------------------------------------------------
-.. admonition:: You will learn ...
+  .. admonition:: What you will learn
-   ... how to analyse the event data using an **unbinned maximum
-   likelihood analysis**.
+     You will learn how to analyse the event data using an **unbinned maximum
+     likelihood analysis**.
-   Use an unbinned maximum likelihood analysis if you have doubts about the
-   impact of the selected binning on your analysis and if the number of
-   observations is not too large (typically a few tens of 30 minute
-   observations). Or if you simply to not worry about computation time.
+     Use an unbinned maximum likelihood analysis if you have doubts about the
+     impact of the selected binning on your analysis and if the number of
+     observations is not too large (typically a few tens of 30 minute
+     observations). Or if you simply to not worry about computation time.
 Instead of binning the events into a counts cube and computing the effective
 :ref:`instrument response functions <glossary_irf>`
@@ -7,7 +7,7 @@ The following sections explain how you can address some specific analysis
 issues.
 .. toctree::
-   :maxdepth: 2
+   :maxdepth: 1
    howto_display
    howto_ts
@@ -3,10 +3,13 @@
 How to display the results?
 ---------------------------
-.. admonition:: You will learn ...
+  .. admonition:: What you will learn
-   ... how to **display results** using some example scripts written in
-   Python and requiring the ``matplotlib`` module.
+     You will learn how to **display results** using some example scripts
+     written in Python.
+
+     Please not that these scripts need the ``matplotlib`` Python module
+     installed.
 The ctools package does not contain any tools or scripts for graphical
 display of results since results are generally written into standard FITS
@@ -29,10 +32,15 @@ Python module. You can find these scripts in the
   +------------------------+-----------------------------+
   | ``show_pointings.py``  | Display pointing directions |
   +------------------------+-----------------------------+
+  | ``show_phases.py``     | Display event phases        |
+  +------------------------+-----------------------------+
   | ``show_spectrum.py``   | Display a spectrum          |
   +------------------------+-----------------------------+
-Below some usage examples and the expected output:
+Below some usage examples and the expected output.
+
+show_butterfly.py
+^^^^^^^^^^^^^^^^^
 .. code-block:: bash
@@ -44,6 +52,10 @@ Below some usage examples and the expected output:
    *Butterfly diagram displayed with show_butterfly.py*
+
+show_lightcurve.py
+^^^^^^^^^^^^^^^^^^
+
 .. code-block:: bash
    $ $CTOOLS/share/examples/python/show_lightcurve.py lightcurve.fits
@@ -54,6 +66,10 @@ Below some usage examples and the expected output:
    *Light curve displayed with show_lightcurve.py*
+
+show_obs.py
+^^^^^^^^^^^
+
 .. code-block:: bash
    $ $CTOOLS/share/examples/python/show_obs.py obs_selected.xml
@@ -64,6 +80,10 @@ Below some usage examples and the expected output:
    *Observation summary displayed with show_obs.py*
+
+show_pointings.py
+^^^^^^^^^^^^^^^^^
+
 .. code-block:: bash
    $ $CTOOLS/share/examples/python/show_pointings.py obs_selected.xml
@@ -74,6 +94,24 @@ Below some usage examples and the expected output:
    *Observation summary displayed with show_pointings.py (zoomed in)*
+
+show_phases.py
+^^^^^^^^^^^^^^
+
+.. code-block:: bash
+
+   $ $CTOOLS/share/examples/python/show_phases.py -n 50 events_phased.fits
+
+.. figure:: howto_display_phases.png
+   :width: 600px
+   :align: center
+
+   *Event phases displayed with show_phases.py*
+
+
+show_spectrum.py
+^^^^^^^^^^^^^^^^
+
 .. code-block:: bash
    $ $CTOOLS/share/examples/python/show_spectrum.py spectrum_src001.fits
 .. _1dc_howto_edisp:
-How to take the energy dispersion into account
-----------------------------------------------
+How to take the energy dispersion into account?
+-----------------------------------------------
-.. admonition:: You will learn ...
+  .. admonition:: What you will learn
-   ... how to **take into account the energy dispersion** in a maximum
-   likelihood analysis.
+     You will learn how to **take the energy dispersion into account** when
+     fitting a model to the data using a maximum likelihood analysis.
-   Although the effect of the energy dispersion is often neglegible there
-   may be cases where you want to consider energy dispersion in an analysis,
-   for example if you are analysing the data down to very low energies.
-   Energy dispersion is fully implemented in ctools but is computationally
-   intensive. So be aware that the tools and scripts will take longer.
+     Although the effect of the energy dispersion is often neglegible there
+     may be cases where you want to consider energy dispersion in an analysis,
+     for example if you are analysing the data down to very low energies.
+     Energy dispersion is fully implemented in ctools but is computationally
+     intensive. So be aware that the tools and scripts will take longer.
 If you are doing a stacked analysis the first thing you need is an energy
 dispersion cube. You generate this cube using the :ref:`ctedispcube` tool:
@@ -34,7 +34,7 @@ dispersion cube. You generate this cube using the :ref:`ctedispcube` tool:
    Number of energy bins [20] 30
    Output energy dispersion cube file [edispcube.fits]
-This produces a
+This produces an
 :ref:`energy dispersion cube <glossary_edispcube>`
 FITS file that contains the weighted energy dispersion as function of
 sky position and energy.
@@ -53,9 +53,10 @@ during the fit:
    Input event list, counts cube or observation definition XML file [events.fits] cntcube.fits
    Input exposure cube file (only needed for stacked analysis) [NONE] expcube.fits
    Input PSF cube file (only needed for stacked analysis) [NONE] psfcube.fits
+   Input energy dispersion cube file (only needed for stacked analysis) [NONE] edispcube.fits
    Input background cube file (only needed for stacked analysis) [NONE] bkgcube.fits
    Input model definition XML file [$CTOOLS/share/models/crab.xml] stacked_models.xml
-   Output model definition XML file [crab_results.xml] stacked_results.xml
+   Output model definition XML file [crab_results.xml] stacked_results_edisp.xml
 If you are doing an unbinned analysis you do not need to generate an energy
 dispersion cube and you can directly run :ref:`ctlike` with the ``edisp=yes``
 .. _1dc_howto_exclude:
-How to exclude pixels from a counts cube?
------------------------------------------
+How to exclude bins from a counts cube?
+---------------------------------------
-.. admonition:: You will learn ...
+  .. admonition:: What you will learn
-   ... how to **exclude pixels of a counts cube** from an analysis using
-   the :ref:`ctcubemask` tool.
+     You will learn how to **exclude bins from a counts cube** to restrict an
+     analysis to a specific region of interest.
 You can exclude bins in an analysis from a counts cube by using the
 :ref:`ctcubemask` tool.
@@ -42,3 +42,38 @@ in the following image.
    :align: center
    *Filtered counts cube*
+
+In the next step we exclude in addition a circular region of 0.3 degrees
+around ``Src001`` from the counts cube. For this purpose you first have to
+generate a
+`ds9 region file <http://ds9.si.edu/doc/ref/region.html>`_ file:
+
+.. code-block:: bash
+
+   $ nano exclude.reg
+   fk5
+   circle(266.42, -29.01, 0.3)
+
+and use this file as ``Input exclusion region file in ds9 format`` in
+the :ref:`ctcubemask` run:
+
+.. code-block:: bash
+
+   $ ctcubemask
+   Input counts cube or observation definition XML file [cntcube.fits]
+   Input exclusion region file in ds9 format [NONE] exclude.reg
+   RA for ROI centre (degrees) (0-360) [83.63] 266.42
+   Dec for ROI centre (degrees) (-90-90) [22.01] -29.01
+   Radius of ROI (degrees) (0-180) [3.0] 1.0
+   Lower energy limit (TeV) [0.1]
+   Upper energy limit (TeV) [100.0]
+   Output counts cube or observation definition XML file [cntcube_filtered.fits] cntcube_filtered_exclude.fits
+
+The resulting counts cube, summed over the energy axis, is shown
+in the following image.
+
+.. figure:: howto_filtered_exclude.png
+   :width: 400px
+   :align: center
+
+   *Filtered counts cube with exclusion region*
@@ -3,10 +3,10 @@
 How to determine the extension of a source?
 -------------------------------------------
-.. admonition:: You will learn ...
+  .. admonition:: What you will learn
-   ... how you **determine the spatial extent of a source** by fitting an
-   extended source model to the data.
+     You will learn how you **determine the spatial extent of a source** by
+     fitting an extended source model to the data.
 To determine the spatial extension of a source, an extended source model
 has to be used in the
@@ -15,7 +15,7 @@ of the source of interest.
 In the example below, a ``RadialGaussian`` model is used for the source
 ``Src001`` to model the spatial morphology of the source by an axisymmetric
 Gaussian profile with the Gaussian width ``sigma`` as free parameter.
-Note that at the same time the source type has been set to
+Note that at the same time the source type was set to
 ``ExtendedSource``.
 .. code-block:: xml
@@ -41,31 +41,31 @@ log file and will also be written into the output
 :ref:`model definition file <glossary_moddef>`.
 For illustration, an excerpt of the :ref:`ctlike` log file is shown below.
-The Gaussian width of ``Src001`` has been fitted to 0.0218 +/- 0.0006 degrees.
+The Gaussian width of ``Src001`` has been fitted to 0.049 +/- 0.001 degrees.
 .. code-block:: xml
-   2017-03-04T17:44:06: === GModelSky ===
-   2017-03-04T17:44:06:  Name ......................: Src001
-   2017-03-04T17:44:06:  Instruments ...............: all
-   2017-03-04T17:44:06:  Instrument scale factors ..: unity
-   2017-03-04T17:44:06:  Observation identifiers ...: all
-   2017-03-04T17:44:06:  Model type ................: ExtendedSource
-   2017-03-04T17:44:06:  Model components ..........: "RadialGaussian" * "PowerLaw" * "Constant"
-   2017-03-04T17:44:06:  Number of parameters ......: 7
-   2017-03-04T17:44:06:  Number of spatial par's ...: 3
-   2017-03-04T17:44:06:   RA .......................: 266.416716505268 +/- 0.000753191747203773 [-360,360] deg (free,scale=1)
-   2017-03-04T17:44:06:   DEC ......................: -29.0044104159472 +/- 0.000655420291471696 [-90,90] deg (free,scale=1)
-   2017-03-04T17:44:06:   Sigma ....................: 0.02183871432113 +/- 0.000636651479386333 [0.001,10] deg (free,scale=1)
-   2017-03-04T17:44:06:  Number of spectral par's ..: 3
-   2017-03-04T17:44:06:   Prefactor ................: 4.85894838349519e-17 +/- 7.40461052122477e-19 [0,infty[ ph/cm2/s/MeV (free,scale=5.7e-18,gradient)
-   2017-03-04T17:44:06:   Index ....................: -2.29168205335762 +/- 0.00909936860902909 [10,-10]  (free,scale=-2.48,gradient)
-   2017-03-04T17:44:06:   PivotEnergy ..............: 300000 MeV (fixed,scale=300000,gradient)
-   2017-03-04T17:44:06:  Number of temporal par's ..: 1
-   2017-03-04T17:44:06:   Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
+   2017-06-03T16:09:45: === GModelSky ===
+   2017-06-03T16:09:45:  Name ......................: Src001
+   2017-06-03T16:09:45:  Instruments ...............: all
+   2017-06-03T16:09:45:  Instrument scale factors ..: unity
+   2017-06-03T16:09:45:  Observation identifiers ...: all
+   2017-06-03T16:09:45:  Model type ................: ExtendedSource
+   2017-06-03T16:09:45:  Model components ..........: "RadialGaussian" * "ExponentialCutoffPowerLaw" * "Constant"
+   2017-06-03T16:09:45:  Number of parameters ......: 8
+   2017-06-03T16:09:45:  Number of spatial par's ...: 3
+   2017-06-03T16:09:45:   RA .......................: 266.420293012532 +/- 0.00132113631921452 [-360,360] deg (free,scale=1)
+   2017-06-03T16:09:45:   DEC ......................: -29.0055014958646 +/- 0.00116849557635297 [-90,90] deg (free,scale=1)
+   2017-06-03T16:09:45:   Sigma ....................: 0.0488254879944832 +/- 0.000993438752316212 [0.001,10] deg (free,scale=1)
+   2017-06-03T16:09:45:  Number of spectral par's ..: 4
+   2017-06-03T16:09:45:   Prefactor ................: 3.2695646853429e-17 +/- 6.35060270585935e-19 [1e-25,1e-15] ph/cm2/s/MeV (free,scale=1e-18,gradient)
+   2017-06-03T16:09:45:   Index ....................: -2.25022445798271 +/- 0.0188890753731813 [-0,-5]  (free,scale=-1,gradient)
+   2017-06-03T16:09:45:   CutoffEnergy .............: 28286052.9486927 +/- 6883002.00397237 [100000,10000000000] MeV (free,scale=10000000,gradient)
+   2017-06-03T16:09:45:   PivotEnergy ..............: 300000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)
+   2017-06-03T16:09:45:  Number of temporal par's ..: 1
+   2017-06-03T16:09:45:   Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
 .. warning::
-   The errors returned by :ref:`ctlike` are purely statistical. In addition
-   there are systematic uncertainties, such as for example the limits in the
-   knowledge of the point spread function. It remains to be seen whether CTA
-   will be able to reliably determine source extents as small as 0.02 degrees.
+   The parameter errors returned by :ref:`ctlike` are purely statistical. In
+   addition to the statistical errors there are systematic uncertainties, such
+   as for example the limits on the knowledge of the point spread function.
@@ -3,10 +3,10 @@
 How to generate a light curve?
 ------------------------------
-.. admonition:: You will learn ...
+  .. admonition:: What you will learn
-   ... how to **generate a light curve** for a source using the
-   :ref:`cslightcrv` script.
+     You will learn how to **generate a light curve** for a source to study
+     the temporal variation of the flux and spectrum.
 You can generate a light curve using the :ref:`cslightcrv` script which
 takes on input an
@@ -25,7 +25,7 @@ The following example illustrates how you can generate a light curve for
 corresponding to the dates from 2021-01-21 (12:00:00) to 2021-01-22 (00:00:00)
 (to convert dates to MJD and vice versa we suggest to use the
 `xTime Date and Time Conversion Utility <https://heasarc.gsfc.nasa.gov/cgi-bin/Tools/xTime/xTime.pl>`_).
-The time interval is linearly divided into 100 time bins.
+The time interval is linearly divided into 20 time bins.
 .. code-block:: bash
 .. _1dc_howto_phasecurve:
-How to generate a phase curve of a pulsar or a gamma-ray binary?
-----------------------------------------------------------------
-
-So far there are no tools in ctools to generate phase curves of pulsars or
-gamma-ray binaries. A plugin for
-`tempo2 <http://www.atnf.csiro.au/research/pulsar/tempo2/>`_
-is under development by the ctools developers, and you may get in touch
-with Thierry Reposeur (reposeur@in2p3.fr) for information about this
-development.
-
-If you know the phase information of the pulsar or gamma-ray binary,
-you may specify a ``temporal`` model component of type ``PhaseCurve`` in the
+How to generate a phase curve of a pulsar?
+------------------------------------------
+
+  .. admonition:: What you will learn
+
+     You will learn how to **derive a phase curve for a pulsar**.
+
+To determine the phase curve of a pulsar you first need the ephemerides
+of the pulsar of interest. In this example we suppose that you want to
+analyse the Vela Pulsar which for an epoch of 54686.2 (MJD) has the
+following ephemerides:
+
+* :math:`f = 11.19 \, {\rm s}^{-1}`
+* :math:`\dot{f} = -1.55 \, 10^{-11} \, {\rm s}^{-2}`
+* :math:`\ddot{f} = 6.46 \, 10^{-23} \, {\rm s}^{-3}`
+
+You need to enter these ephemerides into a
 :ref:`model definition file <glossary_moddef>`
-for the source and fit the information using an unbinned analysis.
-An example of a
+with the following structure (note that the spectral and spatial information
+is irrelevant for this example):
+
+.. code-block:: xml
+
+   <?xml version="1.0" encoding="UTF-8" standalone="no"?>
+   <source_library title="source library">
+     <source name="VelaPulsar" type="PointSource">
+       <spectrum type="PowerLaw">
+         <parameter name="Prefactor"   scale="1e-14" value="19.4" min="1e-07" max="1000.0" free="1"/>
+         <parameter name="Index"       scale="-1"    value="4.45" min="0.0"   max="+10.0"  free="1"/>
+         <parameter name="PivotEnergy" scale="1e3"   value="20.0" min="0.01"  max="1000.0" free="0"/>
+       </spectrum>
+       <spatialModel type="PointSource">
+         <parameter name="RA"  scale="1" value="128.838" min="-360" max="360" free="0"/>
+         <parameter name="DEC" scale="1" value="-45.178" min="-90"  max="90"  free="0"/>
+       </spatialModel>
+       <temporal type="PhaseCurve" file="$CTADATA/models/phasecurve_J0835-4510.fits">
+         <parameter name="Normalization" scale="1"         value="1.0"     min="0.0"   max="1000.0"   free="0"/>
+         <parameter name="MJD"           scale="1"         value="54686.2" min="0.0"   max="100000.0" free="0"/>
+         <parameter name="Phase"         scale="1"         value="0.0"     min="0.0"   max="1.0"      free="0"/>
+         <parameter name="F0"            scale="1"         value="11.19"   min="0.0"   max="1000.0"   free="0"/>
+         <parameter name="F1"            scale="1e-11"     value="-1.55"   min="-10.0" max="1000.0"   free="0"/>
+         <parameter name="F2"            scale="1e-23"     value="6.46"    min="0.0"   max="1000.0"   free="0"/>
+       </temporal>
+     </source>
+   </source_library>
+
+As next step you need to select all observations from the data with pointing
+directions close to the Vela pulsar. In the example below, all observations
+from the Galactic Plane Scan with pointing directions within 5 degrees of the
+Vela pulsar are selected:
+
+.. code-block:: bash
+
+   $ csobsselect
+   Input event list or observation definition XML file [obs.xml] $CTADATA/obs/obs_gps_baseline.xml
+   Pointing selection region shape (CIRCLE|BOX) [CIRCLE]
+   Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL]
+   Galactic longitude of selection centre (deg) (0-360) [184.56] 128.838
+   Galactic latitude of selection centre (deg) (-90-90) [-5.79] -45.178
+   Radius of selection circle (deg) (0-180) [5.0]
+   Output observation definition XML file [outobs.xml] obs_vela.xml
+
+This results in a selection of 45 observations from the Galactic Plane Survey
+that will be compiled in the ``obs_vela.xml``
+:ref:`observation definition file <glossary_obsdef>`.
+Now you can compute the phase for every event in your selection using the
+ephemerides of the Vela pulsar. You do this with the :ref:`ctphase` tool
+as follows:
+
+.. code-block:: bash
+
+   $ ctphase
+   Input event list or observation definition XML file [events.fits] obs_vela.xml
+   Model definition XML file [$CTOOLS/share/models/crab.xml] models.xml
+   Source in model definition XML file to use for phase computation [Crab] VelaPulsar
+   Output event list or observation definition XML file [phased_events.fits] obs_vela_phased.xml
+
+This produces a new
+:ref:`observation definition file <glossary_obsdef>`
+``obs_vela_phased.xml`` that points to event files to which a ``PHASE``
+column was added. The format of the file is dislayed below. Each event
+has now an associated phase value that was computed from the trigger time of
+the event.
+
+.. figure:: howto_phasecurve_events.png
+   :width: 700px
+   :align: center
+
+   *Event list file with phase information for the Vela pulsar*
+
+To generate now a phase curve of the Vela pulsar you need to select all events
+that are close to the source using the :ref:`ctselect` tool. To select for
+example all events within 0.5 degrees of the Vela pulsar you need to type
+the following:
+
+.. code-block:: bash
+
+   $ ctselect
+   Input event list or observation definition XML file [events.fits] obs_vela_phased.xml
+   RA for ROI centre (degrees) (0-360) [83.63] 128.838
+   Dec for ROI centre (degrees) (-90-90) [22.01] -45.178
+   Radius of ROI (degrees) (0-180) [3.0] 0.5
+   Start time (CTA MET in seconds) [0.0]
+   End time (CTA MET in seconds) [0.0]
+   Lower energy limit (TeV) [0.1] 0.03
+   Upper energy limit (TeV) [100.0] 160.0
+   Output event list or observation definition XML file [selected_events.fits] obs_vela_phased_selected.xml
+
+The selected events as a function of phase are displayed below. There is no
+obvious indication of a phase dependence of the event rate.
+
+.. figure:: howto_phasecurve_phases.png
+   :width: 600px
+   :align: center
+
+   *Phase histogram for events selected within 0.5 degrees of the Vela pulsar*
+
+An alternative method to select events from the
+:ref:`observation definition file <glossary_obsdef>`
+is based on the computation of the probability that the event arises from a
+specific model component. To perform such a computation you have to
+
+* add a background model component, and
+* drop the temporal component from the Vela pulsar
+
+in the
+:ref:`model definition file <glossary_moddef>`. The resulting
 :ref:`model definition file <glossary_moddef>`
-including a ``PhaseCurve`` model component is shown below.
-For more information, please read the following
-:ref:`description of the temporal model components <sec_temporal_models>`.
+should look as follows:
 .. code-block:: xml
-   <?xml version="1.0" standalone="no"?>
+   <?xml version="1.0" encoding="UTF-8" standalone="no"?>
    <source_library title="source library">
-     <source name="Crab" type="PointSource">
+     <source name="VelaPulsar" type="PointSource">
+       <spectrum type="PowerLaw">
+         <parameter name="Prefactor"   scale="1e-14" value="19.4" min="1e-07" max="1000.0" free="1"/>
+         <parameter name="Index"       scale="-1"    value="4.45" min="0.0"   max="+10.0"  free="1"/>
+         <parameter name="PivotEnergy" scale="1e3"   value="20.0" min="0.01"  max="1000.0" free="0"/>
+       </spectrum>
        <spatialModel type="PointSource">
-         <parameter name="RA"  scale="1" value="83.6331" min="-360" max="360" free="0"/>
-         <parameter name="DEC" scale="1" value="22.0145" min="-90"  max="90"  free="0"/>
+         <parameter name="RA"  scale="1" value="128.838" min="-360" max="360" free="0"/>
+         <parameter name="DEC" scale="1" value="-45.178" min="-90"  max="90"  free="0"/>
        </spatialModel>
+     </source>
+     <source name="Background" type="CTAIrfBackground">
        <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"/>
+         <parameter name="Prefactor" value="1" error="0" scale="1" min="0" free="1" />
+         <parameter name="Index" value="0" error="0" scale="1" min="-10" max="10" free="1" />
+         <parameter name="PivotEnergy" value="1" scale="1000000" free="0" />
        </spectrum>
-       <temporal type="PhaseCurve" file="model_temporal_phasecurve.fits">
-         <parameter name="Normalization" scale="1" value="1.0"     min="0.0" max="1000.0"   free="0"/>
-         <parameter name="MJD"           scale="1" value="51544.5" min="0.0" max="100000.0" free="0"/>
-         <parameter name="Phase"         scale="1" value="0.0"     min="0.0" max="1.0"      free="0"/>
-         <parameter name="F0"            scale="1" value="1.0"     min="0.0" max="1000.0"   free="0"/>
-         <parameter name="F1"            scale="1" value="0.1"     min="0.0" max="1000.0"   free="0"/>
-         <parameter name="F2"            scale="1" value="0.01"    min="0.0" max="1000.0"   free="0"/>
-       </temporal>
      </source>
    </source_library>
+You can then run the :ref:`ctprob` tool that computes the probability that an
+event arises either from the source or from the background component:
+
+.. code-block:: bash
+
+   $ ctprob
+   Input event list or observation definition XML file [events.fits] obs_vela_phased.xml
+   Input model definition XML file [$CTOOLS/share/models/crab.xml] models_prob.xml
+   Output event list or observation definition XML file [prob_events.fits] obs_vela_phased_prob.xml
+
+This produces a new
+:ref:`observation definition file <glossary_obsdef>`
+``obs_vela_phased_prob.xml`` that points to event files to which the
+``PROB_VelaPulsar`` and ``PROB_Background`` columns have been added.
+The format of the file is dislayed below. Each event has now an associated
+probability value for all model components.
+
+.. figure:: howto_phasecurve_prob.png
+   :width: 800px
+   :align: center
+
+   *Event list file with phase and probability information for the Vela pulsar and probability information for the background model*
+
+You now can perform an event selection based on the probability that an event
+originates from the Vela pulsar. In the following example we set the
+probability threshold to >20% by adding the attribute
+``expr="PROB_VelaPulsar>0.2"`` after the :ref:`ctselect` command:
+
+.. code-block:: bash
+
+   $ ctselect expr="PROB_VelaPulsar>0.2"
+   Input event list or observation definition XML file [obs_vela_phased.xml] obs_vela_phased_prob.xml
+   RA for ROI centre (degrees) (0-360) [128.838]
+   Dec for ROI centre (degrees) (-90-90) [-45.178]
+   Radius of ROI (degrees) (0-180) [0.5] 5.0
+   Start time (CTA MET in seconds) [0.0]
+   End time (CTA MET in seconds) [0.0]
+   Lower energy limit (TeV) [0.03]
+   Upper energy limit (TeV) [160.0]
+   Output event list or observation definition XML file [obs_vela_phased_selected.xml] obs_vela_phased_prob_selected.xml
+
+The selected events as a function of phase are displayed below. Unfortunately
+there is still no indication of a phase dependence of the event rate.
+
+.. figure:: howto_phasecurve_probphases.png
+   :width: 600px
+   :align: center
+
+   *Phase histogram of events for events having a probability of >20% in originating from the Vela pulsar*
+
+Since every event in the event files of the Data Challenge is tagged with a
+Monte Carlo identifier that specifies from which model component it was
+simulated, you may check how many events from the Vela pulsar are actually
+expected in the data. The correspondance between the model components and
+the Monte Carlo identifiers is given in the header of the event files
+(see :ref:`1dc_data_organisation`), and from there you will find out that
+the Vela pulsar has the Monte Carlo identifier ``86``. You can then select only
+the events with this identifier by adding the attribute ``expr="MC_ID==86"``
+after the :ref:`ctselect` command:
+
+.. code-block:: bash
+
+   $ ctselect expr="MC_ID==86"
+   Input event list or observation definition XML file [obs_vela_phased_prob.xml] obs_vela_phased.xml
+   RA for ROI centre (degrees) (0-360) [128.838]
+   Dec for ROI centre (degrees) (-90-90) [-45.178]
+   Radius of ROI (degrees) (0-180) [5.0]
+   Start time (CTA MET in seconds) [0.0]
+   End time (CTA MET in seconds) [0.0]
+   Lower energy limit (TeV) [0.03]
+   Upper energy limit (TeV) [160.0]
+
+The resulting phase histogram is shown below. It now becomes obvious that there
+are too few events in the simulated data to allow for a significant detection
+of the pulsar. The events from the Vela pulsar are in fact burried into the
+background events.
+
+.. figure:: howto_phasecurve_mcid.png
+   :width: 700px
+   :align: center
+
+   *Distribution of Vela pulsar events as function of event phase*
+
+.. note::
+   The plots were created using the ``show_phases.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 plots by typing
+
+   .. code-block:: bash
+
+      $ $CTOOLS/share/examples/python/show_phases.py obs_vela_phased_selected.xml
+      $ $CTOOLS/share/examples/python/show_phases.py obs_vela_phased_prob_selected.xml
+      $ $CTOOLS/share/examples/python/show_phases.py obs_vela_phased_mcid.xml
@@ -3,15 +3,15 @@
 How to compute the signifiance of a source?
 -------------------------------------------
-.. admonition:: You will learn ...
+  .. admonition:: What you will learn
-   ... how to **compute the significance of a source detection** using the
-   :ref:`ctlike` tool.
+     You will learn how to **compute the significance of a source detection**
+     using the :ref:`ctlike` tool.
 To compute the detection significance of a given source you need to add the
 attribute ``tscalc="1"`` to the ``<source>`` tag in the
 :ref:`model definition file <glossary_moddef>`
-and execute :ref:`ctlike` using this file.
+and execute :ref:`ctlike` using the modified file.
 The structure of the XML file with significance computation requested for
 ``Src001`` is shown below.
@@ -43,7 +43,7 @@ freedom, where :math:`p` is the number of free parameters of ``Src001``.
 For :math:`p=1`, the source significance is the square root of the Test
 Statistic.
-An excerpt of the :ref:`ctlike` log file is shown below where signifance
+An excerpt of the :ref:`ctlike` log file is shown below where significance
 computation was requested for ``Src001`` and ``Src002``. The computation was
 done using a stacked analysis.
@@ -3,10 +3,10 @@
 How to generate a Test Statistic map?
 -------------------------------------
-.. admonition:: You will learn ...
+  .. admonition:: What you will learn
-   ... how to **generate a map of source detection significances** using the
-   :ref:`cttsmap` tool.
+     You will learn how to **generate a map of source detection significances**
+     to localise a gamma-ray source.
 If you want to know how the Test Statistic changes with source position you
 can create a Test Statistic map using the :ref:`cttsmap` tool. A Test
@@ -56,12 +56,12 @@ The figure below shows the resulting Test Statistic map.
    The ``binmin`` and ``binmax`` parameters can be used to compute only a
    subset of the bins in a map (see :ref:`cttsmap`). This functionality is
-   used by the :ref:`cstsmapsplit` to generate a command sequence that is
+   used by the :ref:`cstsmapsplit` script to generate a command sequence that is
    useful for parallelizing the Test Statistic map computation. The
    :ref:`cstsmapmerge` script can then be used to recombine the different
    partial maps into a single Test Statistic map.
-   To produce the same result that has been produced above by the single
+   To produce the same result that was produced above by the single
    :ref:`cttsmap` run, now however split over ten :ref:`cttsmap` jobs that
    can be executed in parallel, type
@@ -86,7 +86,7 @@ The figure below shows the resulting Test Statistic map.
       Compute null hypothesis first? [yes]
       ASCII file containing all commands [commands.dat] tsmap_commands.dat
-   The :ref:`cstsmapmerge` script creates the following ASCII file on output:
+   The :ref:`cstsmapsplit` script creates the following ASCII file on output:
    .. code-block:: bash
@@ -107,9 +107,9 @@ The figure below shows the resulting Test Statistic map.
       ./tsmap_commands.dat
-   This will execute ten :ref:`cttsmap` jobs and parallel that each creates
-   ten pixels of the Test Statistic map. You combine the resulting ten Test
-   Statistic maps into a single map by typing
+   This will execute 10 :ref:`cttsmap` jobs in parallel. Each of the jobs will
+   compute 10 pixels of the Test Statistic map. You combine the resulting 10
+   Test Statistic maps into a single map by typing
    .. code-block:: bash
@@ -3,10 +3,10 @@
 How to compute upper limits?
 ----------------------------
-.. admonition:: You will learn ...
+  .. admonition:: What you will learn
-   ... how to **determine the upper flux limit for a source** using the
-   :ref:`ctulimit` tool.
+     You will learn how to **determine the upper flux limit for a source**
+     using the :ref:`ctulimit` tool.
 You determine upper limits for gamma-ray sources using the :ref:`ctulimit`
 tool.
@@ -34,7 +34,7 @@ for example
      </spatialModel>
    </source>
-and then run the :ref:`ctulimit` as follows:
+and then run the :ref:`ctulimit` tool as follows:
 .. code-block:: bash
@@ -51,12 +51,12 @@ extracted from the :ref:`ctulimit` log file:
 .. code-block:: bash
-   2017-03-04T20:29:08: +=====================+
-   2017-03-04T20:29:08: | Upper limit results |
-   2017-03-04T20:29:08: +=====================+
-   2017-03-04T20:29:08:  Differential flux limit ...: 1.01134913337855e-20 ph/cm2/s/MeV at 1 TeV
-   2017-03-04T20:29:08:  Integral flux limit .......: 6.82594737233183e-15 ph/cm2/s within [1-100] TeV
-   2017-03-04T20:29:08:  Energy flux limit .........: 3.00560602939917e-14 erg/cm2/s within [1-100] TeV
+   2017-06-02T17:00:09: +=====================+
+   2017-06-02T17:00:09: | Upper limit results |
+   2017-06-02T17:00:09: +=====================+
+   2017-06-02T17:00:09:  Differential flux limit ...: 4.50187740703557e-20 ph/cm2/s/MeV at 1 TeV
+   2017-06-02T17:00:09:  Integral flux limit .......: 3.03847378149801e-14 ph/cm2/s within [1-100] TeV
+   2017-06-02T17:00:09:  Energy flux limit .........: 1.33790294880664e-13 erg/cm2/s within [1-100] TeV
 .. warning::
    It is important to fix the spectral index of the test source in the
@@ -3,12 +3,12 @@
 How to use Virtual Observatory tools with ctools?
 -------------------------------------------------
-.. admonition:: You will learn ...
+  .. admonition:: What you will learn
-   ... how to interaction with Virtual Observatory tools such as
-   `Aladin <http://aladin.u-strasbg.fr>`_
-   or
-   `Topcat <http://www.star.bris.ac.uk/~mbt/topcat/>`_.
+     You will learn how to interact with the Virtual Observatory tools
+     `Aladin <http://aladin.u-strasbg.fr>`_ for image display and
+     `Topcat <http://www.star.bris.ac.uk/~mbt/topcat/>`_ for display of tables
+     such as event lists.
 Several ctools can interact with Virtual Observatory tools such as
 `Aladin <http://aladin.u-strasbg.fr>`_
@@ -117,9 +117,7 @@ energy (see below).
    *Display of event list in 3-dimensional representation*
 Below the list of tools and scripts that interact with Virtual Observatory
-tools. Some tools or scripts publish tables instead of images, and the tables
-can for example by displayed and explored using
-`Topcat <http://www.star.bris.ac.uk/~mbt/topcat/>`_.
+tools.
   +-------------------+--------+-------------------+
   | Tool or script    | Object | VO tool (example) |
@@ -146,6 +144,10 @@ can for example by displayed and explored using
   +-------------------+--------+-------------------+
   | :ref:`ctobssim`   | Table  | Topcat            |
   +-------------------+--------+-------------------+
+  | :ref:`ctphase`    | Table  | Topcat            |
+  +-------------------+--------+-------------------+
+  | :ref:`ctprob`     | Table  | Topcat            |
+  +-------------------+--------+-------------------+
   | :ref:`ctselect`   | Table  | Topcat            |
   +-------------------+--------+-------------------+
   | :ref:`ctskymap`   | Image  | Aladin            |