Update ctools User manual

Jürgen Knödlseder
1 parent 03ce9b01
Showing 15 changed files with 1452 additions and 1134 deletions
doc/source/users/user_manual/models.rst
doc/source/users/user_manual/models_concept.rst
doc/source/users/user_manual/models_implementation.rst
doc/source/users/user_manual/models_spatial_bgd.rst
doc/source/users/user_manual/models_spatial_src.rst
doc/source/users/user_manual/models_spectral.rst
doc/source/users/user_manual/models_temporal.rst
doc/source/users/user_manual/plotting.rst
doc/source/users/user_manual/plotting_butterfly.png
doc/source/users/user_manual/plotting_lightcurve.png
doc/source/users/user_manual/plotting_obs.png
doc/source/users/user_manual/plotting_phases.png
doc/source/users/user_manual/plotting_pointings.png
doc/source/users/user_manual/plotting_spectrum.png
doc/source/users/user_manual/response.rst
@@ -3,1066 +3,12 @@
 Modelling CTA data
 ------------------
  
-.. _sec_model_concept:
-
-Overall concept
-~~~~~~~~~~~~~~~
-
-A central element of any data analysis is the modeling of the observed 
-event distribution.
-Generally, the measured events can be separated into two distinct classes:
-events being attributed to gamma rays of celestial origin (the "source" 
-events) or events being attributed to any other kind of trigger (the 
-"background" events).
-In the first case, the :ref:`Instrument Response Functions (IRFs) <sec_response>`
-describe how an incident gamma-ray converts into a measured event.
-In the second case there is often no general prescription, and the 
-distribution of background events is commonly modeled directly in the data 
-space of the observable quantities.
-The main difficulty of gamma-ray astronomy is that not all events can be
-tagged a priori as being source or background, hence their separation can 
-only be done on a statistical basis, based on the different morphologies, 
-spectral characteristics or eventually temporal signatures of both event
-categories.
-
-For this purpose, ctools use a general model that describes the spatial, 
-spectral and temporal properties of the source and background components.
-The model is composed of an arbitrary number of model components that
-add up linearly to provide a prediction of the expected number of events.
-Model components are generally parametric, and model parameters can be 
-adjusted through a maximum likelihood procedure to find the set of 
-parameters that represent best the measured data.
-Model components representing celestial sources are convolved with the 
-:ref:`IRFs <sec_response>` to predict the expected number of source events in 
-the data.
-Background model components will be directly expressed as expected number 
-of background events without any :ref:`IRF <sec_response>` convolution.
-
-
-.. _sec_model_implementation:
-
-Implementation
-~~~~~~~~~~~~~~
-
-The general model is describe in ctools using a model definition XML file. 
-Below is a simple example of such a file comprising one source and one 
-background model.
-Each model is factorised into
-a spectral (tag ``<spectrum>``),
-a spatial (tags ``<spatialModel>`` and ``<radialModel>``), and
-a temporal component (tag ``<temporal>``).
-
-.. math::
-  M(x,y,E,t) = M_{\rm spatial}(x,y|E) \times M_{\rm spectral}(E) \times M_{\rm temporal}(t)
-
-In this specific example, the source component ``Crab`` describes 
-a point source at the location of the Crab nebula with a power law spectral
-shape that is constant in time.
-The background component ``Background`` is modelled as a radial Gaussian 
-function in offset angle squared (with the offset angle being defined as 
-the angle between pointing and measured event direction) and a spectral
-function that is tabulated in an ASCII file.
-
-.. code-block:: xml
-
-  <?xml version="1.0" standalone="no"?>
-  <source_library title="source library">
-    <source name="Crab" type="PointSource">
-      <spatialModel type="PointSource">
-        <parameter name="RA"  scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
-        <parameter name="DEC" scale="1.0" value="22.0145" min="-90"  max="90"  free="1"/>
-      </spatialModel>
-      <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>
-      <temporal type="Constant">
-        <parameter name="Normalization" scale="1.0" value="1.0" min="0.0" max="1000.0" free="0"/>
-      </temporal>
-    </source>
-    <source name="Background" type="RadialAcceptance" instrument="CTA">
-      <spectrum type="FileFunction" file="$CTOOLS/share/models/bkg_dummy.txt">
-        <parameter name="Normalization" scale="1.0" value="1.0" min="0.0" max="1000.0" free="1"/>
-      </spectrum>
-      <radialModel type="Gaussian">
-         <parameter name="Sigma" scale="1.0" value="3.0" min="0.01" max="10.0" free="1"/>
-      </radialModel>
-    </source>
-  </source_library>
-
-Model parameters are specified by a ``<parameter>`` tag with a certain 
-number of attributes.
-The ``name`` attribute gives the (case-sensitive) parameter name that 
-needs to be unique within a given model component.
-The ``scale`` attribute gives a scaling factor that will be multiplied by 
-the ``value`` attribute to provide the real (physical) parameter value.
-The ``min`` and ``max`` attributes specify boundaries for the ``value``
-term of the parameter.
-And the ``free`` attribute specifies whether a parameter should be fitted 
-(``free="1"``) or kept fixed (``free="0"``) in a maximum likelihood 
-analysis.
-After a maximum likelihood fit using :ref:`ctlike`, an
-``error`` attribute giving the statisical uncertainty of the ``value``
-term will be appended to each ``<parameter>`` tag.
-
-.. note::
-   For compatibility reasons with the Fermi/LAT XML format the ``<temporal>``
-   tag can be omitted for models that are constant in time.
-
-.. note::
-
-   XML files are ASCII files and can be edited by hand using any text 
-   editor.
-   The indentation of the XML elements is not mandatory.
-
-.. note::
-
-   The splitting of parameter values into a ``value`` and ``scale`` factor 
-   is mainly for numerical purposes.
-   Parameter fitting algorithms can be ill-conditioned if several 
-   parameters of vastly different orders of magnitudes need to be 
-   optimised simultaneously.
-   Splitting a value into two components allows to "prescale" the 
-   parameters so that the effective parameters to be optimised (the ``value`` terms) 
-   are all of about unity.
-
-.. note::
-
-   The syntax of the model definition XML file has been inspired from the
-   syntax used by the Fermi/LAT ScienceTools, but for reasons of clarity and
-   homogenity of the various model and parameter names we have made some
-   modifications.
-   Nevertheless, the format used by the Fermi/LAT ScienceTools is also
-   supported.
-
-
-.. _sec_spatial_src_models:
-
-Spatial source model components
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-The following sections present the spatial model components that are available 
-in ctools for gamma-ray sources.
-
-Point source
-^^^^^^^^^^^^
-
-  .. code-block:: xml
-
-    <source name="Crab" type="PointSource">
-      <spatialModel type="PointSource">
-        <parameter name="RA"  scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
-        <parameter name="DEC" scale="1.0" value="22.0145" min="-90"  max="90"  free="1"/>
-      </spatialModel>
-      <spectrum type="...">
-        ...
-      </spectrum>
-    </source>
-
-  .. note::
-
-    For compatibility with the Fermi/LAT ScienceTools the model type
-    ``PointSource`` can be replaced by ``SkyDirFunction``.
-
-
-Radial source
-^^^^^^^^^^^^^
-
-  .. code-block:: xml
-
-    <source name="Crab" type="ExtendedSource">
-      <spatialModel type="RadialDisk">
-        <parameter name="RA"     scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
-        <parameter name="DEC"    scale="1.0" value="22.0145" min="-90"  max="90"  free="1"/>
-        <parameter name="Radius" scale="1.0" value="0.20"    min="0.01" max="10"  free="1"/>
-      </spatialModel>
-      <spectrum type="...">
-        ...
-      </spectrum>
-    </source>
-
-  .. code-block:: xml
-
-    <source name="Crab" type="ExtendedSource">
-      <spatialModel type="RadialGaussian">
-        <parameter name="RA"    scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
-        <parameter name="DEC"   scale="1.0" value="22.0145" min="-90"  max="90"  free="1"/>
-        <parameter name="Sigma" scale="1.0" value="0.20"    min="0.01" max="10"  free="1"/>
-      </spatialModel>
-      <spectrum type="...">
-        ...
-      </spectrum>
-    </source>
-
-  .. code-block:: xml
-
-    <source name="Crab" type="ExtendedSource">
-      <spatialModel type="RadialShell">
-        <parameter name="RA"     scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
-        <parameter name="DEC"    scale="1.0" value="22.0145" min="-90"  max="90"  free="1"/>
-        <parameter name="Radius" scale="1.0" value="0.30"    min="0.01" max="10"  free="1"/>
-        <parameter name="Width"  scale="1.0" value="0.10"    min="0.01" max="10"  free="1"/>
-      </spatialModel>
-      <spectrum type="...">
-        ...
-      </spectrum>
-    </source>
-
-
-Elliptical source
-^^^^^^^^^^^^^^^^^
-
-  .. code-block:: xml
-
-    <source name="Crab" type="ExtendedSource">
-      <spatialModel type="EllipticalDisk">
-        <parameter name="RA"          scale="1.0" value="83.6331" min="-360"  max="360" free="1"/>
-        <parameter name="DEC"         scale="1.0" value="22.0145" min="-90"   max="90"  free="1"/>
-        <parameter name="PA"          scale="1.0" value="45.0"    min="-360"  max="360" free="1"/>
-        <parameter name="MinorRadius" scale="1.0" value="0.5"     min="0.001" max="10"  free="1"/>
-        <parameter name="MajorRadius" scale="1.0" value="2.0"     min="0.001" max="10"  free="1"/>
-      </spatialModel>
-      <spectrum type="...">
-        ...
-      </spectrum>
-    </source>
-
-  .. code-block:: xml
-
-    <source name="Crab" type="ExtendedSource">
-      <spatialModel type="EllipticalGaussian">
-        <parameter name="RA"          scale="1.0" value="83.6331" min="-360"  max="360" free="1"/>
-        <parameter name="DEC"         scale="1.0" value="22.0145" min="-90"   max="90"  free="1"/>
-        <parameter name="PA"          scale="1.0" value="45.0"    min="-360"  max="360" free="1"/>
-        <parameter name="MinorRadius" scale="1.0" value="0.5"     min="0.001" max="10"  free="1"/>
-        <parameter name="MajorRadius" scale="1.0" value="2.0"     min="0.001" max="10"  free="1"/>
-      </spatialModel>
-      <spectrum type="...">
-        ...
-      </spectrum>
-    </source>
-
-
-Diffuse source
-^^^^^^^^^^^^^^
-
-  .. code-block:: xml
-
-    <source name="Crab" type="DiffuseSource">
-      <spatialModel type="DiffuseIsotropic">
-         <parameter name="Value" scale="1" value="1" min="1"  max="1" free="0"/>
-      </spatialModel>
-      <spectrum type="...">
-        ...
-      </spectrum>
-    </source>
-
-  .. note::
-
-    For compatibility with the Fermi/LAT ScienceTools the model type
-    ``DiffuseIsotropic`` can be replaced by ``ConstantValue``.
-
-  .. code-block:: xml
-
-    <source name="Crab" type="DiffuseSource">
-      <spatialModel type="DiffuseMap" file="map.fits">
-         <parameter name="Normalization" scale="1" value="1" min="0.001" max="1000.0" free="0"/>
-      </spatialModel>
-      <spectrum type="...">
-        ...
-      </spectrum>
-    </source>
-
-  .. note::
-
-    For compatibility with the Fermi/LAT ScienceTools the model type
-    ``DiffuseMap`` can be replaced by ``SpatialMap`` and the parameter
-    ``Normalization`` can be replaced by ``Prefactor``.
-
-  .. code-block:: xml
-
-    <source name="Crab" type="DiffuseSource">
-      <spatialModel type="DiffuseMapCube" file="map_cube.fits">
-        <parameter name="Normalization" scale="1" value="1" min="0.001" max="1000.0" free="0"/>
-      </spatialModel>
-      <spectrum type="...">
-        ...
-      </spectrum>
-    </source>
-
-  .. note::
-
-    For compatibility with the Fermi/LAT ScienceTools the model type
-    ``DiffuseMapCube`` can be replaced by ``MapCubeFunction`` and the parameter
-    ``Normalization`` can be replaced by ``Value``.
-
-
-Composite model
-^^^^^^^^^^^^^^^
-
-  .. code-block:: xml
-
-     <source name="Crab" type="CompositeSource">
-       <spatialModel type="Composite">
-         <spatialModel type="PointSource" component="PointSource">
-           <parameter name="RA"    scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
-           <parameter name="DEC"   scale="1.0" value="22.0145" min="-90"  max="90"  free="1"/>
-         </spatialModel>
-         <spatialModel type="RadialGaussian">
-           <parameter name="RA"    scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
-           <parameter name="DEC"   scale="1.0" value="22.0145" min="-90"  max="90"  free="1"/>
-           <parameter name="Sigma" scale="1.0" value="0.20"    min="0.01" max="10"  free="1"/>
-         </spatialModel>
-       </spatialModel>
-       <spectrum type="...">
-         ...
-       </spectrum>
-     </source>
-
-  This spatial model component implements a composite model that is the
-  sum of an arbitrary number of spatial models, computed using
-
-  .. math::
-     M_{\rm spatial}(x,y|E) = \frac{1}{N} \sum_{i=0}^{N-1} M_{\rm spatial}^{(i)}(x,y|E)
-
-  where :math:`M_{\rm spatial}^{(i)}(x,y|E)` is any spatial model component
-  (including another composite model), and :math:`N` is the number of
-  model components that are combined.
-
-
-.. _sec_spatial_bgd_models:
-
-Spatial background model components
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-The following sections present the spatial model components that are available 
-in ctools for background modelling.
-
-CTA radial background
-^^^^^^^^^^^^^^^^^^^^^
-
-  .. code-block:: xml
-
-    <source name="Background" type="RadialAcceptance" instrument="CTA">
-      <radialModel type="Gaussian">
-        <parameter name="Sigma" scale="1.0" value="3.0" min="0.01" max="10.0" free="1"/>
-      </radialModel>
-      <spectrum type="...">
-        ...
-      </spectrum>
-    </source>
-
-  .. code-block:: xml
-
-    <source name="Background" type="RadialAcceptance" instrument="CTA">
-      <radialModel type="Profile">
-        <parameter name="Width" scale="1.0" value="1.5" min="0.1" max="1000.0" free="1"/>
-        <parameter name="Core"  scale="1.0" value="3.0" min="0.1" max="1000.0" free="1"/>
-        <parameter name="Tail"  scale="1.0" value="5.0" min="0.1" max="1000.0" free="1"/>
-      </radialModel>
-      <spectrum type="...">
-        ...
-      </spectrum>
-    </source>
-
-  .. code-block:: xml
-
-    <source name="Background" type="RadialAcceptance" instrument="CTA">
-      <radialModel type="Polynom">
-        <parameter name="Coeff0" scale="1.0" value="+1.00000"   min="-10.0" max="10.0" free="0"/>
-        <parameter name="Coeff1" scale="1.0" value="-0.1239176" min="-10.0" max="10.0" free="1"/>
-        <parameter name="Coeff2" scale="1.0" value="+0.9751791" min="-10.0" max="10.0" free="1"/>
-        <parameter name="Coeff3" scale="1.0" value="-3.0584577" min="-10.0" max="10.0" free="1"/>
-        <parameter name="Coeff4" scale="1.0" value="+2.9089535" min="-10.0" max="10.0" free="1"/>
-        <parameter name="Coeff5" scale="1.0" value="-1.3535372" min="-10.0" max="10.0" free="1"/>
-        <parameter name="Coeff6" scale="1.0" value="+0.3413752" min="-10.0" max="10.0" free="1"/>
-        <parameter name="Coeff7" scale="1.0" value="-0.0449642" min="-10.0" max="10.0" free="1"/>
-        <parameter name="Coeff8" scale="1.0" value="+0.0024321" min="-10.0" max="10.0" free="1"/>
-      </radialModel>
-      <spectrum type="...">
-        ...
-      </spectrum>
-    </source>
-
-
-CTA IRF background
-^^^^^^^^^^^^^^^^^^
-
-  .. code-block:: xml
-
-    <source name="Background" type="CTAIrfBackground" instrument="CTA">
-      <spectrum type="...">
-        ...
-      </spectrum>
-    </source>
-
-
-CTA cube background
-^^^^^^^^^^^^^^^^^^^
-
-  .. code-block:: xml
-
-    <source name="Background" type="CTACubeBackground" instrument="CTA">
-      <spectrum type="...">
-        ...
-      </spectrum>
-    </source>
-
-
-
-.. _sec_spectral_models:
-
-Spectral model components
-~~~~~~~~~~~~~~~~~~~~~~~~~
-
-The following sections present the spectral model components that are available 
-in ctools.
-
-.. warning::
-
-   Source intensities are generally given in units of
-   :math:`{\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm MeV}^{-1}`.
-
-   An exception to this rule exists for the ``DiffuseMapCube`` spatial
-   model where intensities are unitless and the spectral model presents a
-   relative scaling of the diffuse model cube values.
-
-   If spectral models are used for a background model component, intensity 
-   units are generally given in
-   :math:`{\rm counts}\,\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm MeV}^{-1}\,{\rm sr}^{-1}`
-   and correspond to the on-axis count rate.
-
-   Exceptions to this rule exist for the ``CTAIrfBackground`` and
-   ``CTACubeBackground`` models where intensities are unitless and the spectral
-   model presents a relative scaling of the background model values.
-
-
-Constant
-^^^^^^^^
-
-  .. code-block:: xml
-
-   <spectrum type="Constant">
-     <parameter name="Normalization" scale="1e-16" value="5.7" min="1e-07" max="1000.0" free="1"/>
-   </spectrum>
-
-  This spectral model component implements the constant function
-
-  .. math::
-    M_{\rm spectral}(E) = N_0
-
-  where
-
-  * :math:`N_0` = ``Normalization``
-    :math:`({\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm MeV}^{-1})`
-
-  .. note::
-
-    For compatibility with the Fermi/LAT ScienceTools the model type
-    ``Constant`` can be replaced by ``ConstantValue`` and the parameter
-    ``Normalization`` by ``Value``.
-
-
-Power law
-^^^^^^^^^
-
-  .. code-block:: xml
-
-   <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>
-
-  This spectral model component implements the power law function
-
-  .. math::
-    M_{\rm spectral}(E) = k_0 \left( \frac{E}{E_0} \right)^{\gamma}
-
-  where
-
-  * :math:`k_0` = ``Prefactor``
-    :math:`({\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm MeV}^{-1})`
-  * :math:`\gamma` = ``Index``
-  * :math:`E_0` = ``PivotEnergy``
-    :math:`({\rm MeV})`
-
-  .. warning::
-
-    The ``PivotEnergy`` parameter is not intended to be fitted.
-
-  .. note::
-
-    For compatibility with the Fermi/LAT ScienceTools the parameter
-    ``PivotEnergy`` can be replaced by ``Scale``.
-
-  An alternative power law function that uses the integral photon flux as
-  parameter rather than the Prefactor is specified by
-
-  .. code-block:: xml
-
-   <spectrum type="PowerLaw">
-     <parameter scale="1e-07" name="PhotonFlux" min="1e-07" max="1000.0"    value="1.0" free="1"/>
-     <parameter scale="1.0"   name="Index"      min="-5.0"  max="+5.0"      value="-2.0" free="1"/>
-     <parameter scale="1.0"   name="LowerLimit" min="10.0"  max="1000000.0" value="100.0" free="0"/>
-     <parameter scale="1.0"   name="UpperLimit" min="10.0"  max="1000000.0" value="500000.0" free="0"/>
-   </spectrum>
-
-  This spectral model component implements the power law function
-
-  .. math::
-    M_{\rm spectral}(E) = \frac{N(\gamma+1)E^{\gamma}}
-                               {E_{\rm max}^{\gamma+1} - E_{\rm min}^{\gamma+1}}
-
-  where
-
-  * :math:`N` = ``PhotonFlux``
-    :math:`({\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1})`
-  * :math:`\gamma` = ``Index``
-  * :math:`E_{\rm min}` = ``LowerLimit``
-    :math:`({\rm MeV})`
-  * :math:`E_{\rm max}` = ``UpperLimit``
-    :math:`({\rm MeV})`
-
-  .. warning::
-
-    The ``LowerLimit`` and ``UpperLimit`` parameters are always treated as fixed
-    and the flux given by the ``PhotonFlux`` parameter is computed over the
-    range set by these two parameters.
-    Use of this model allows the errors on the integral flux to be evaluated directly
-    by :ref:`ctlike`.
-
-  .. note::
-
-    For compatibility with the Fermi/LAT ScienceTools the model type
-    ``PowerLaw`` can be replaced by ``PowerLaw2`` and the parameter
-    ``PhotonFlux`` by ``Integral``.
-
-
-Exponentially cut-off power law
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-
-  .. code-block:: xml
-
-   <spectrum type="ExponentialCutoffPowerLaw">
-     <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="CutoffEnergy" scale="1e6"   value="1.0"  min="0.01"  max="1000.0" free="1"/>
-     <parameter name="PivotEnergy"  scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="0"/>
-   </spectrum>
-
-  This spectral model component implements the exponentially cut-off power law
-  function
-
-  .. math::
-    M_{\rm spectral}(E) = k_0 \left( \frac{E}{E_0} \right)^{\gamma}
-                          \exp \left( \frac{-E}{E_{\rm cut}} \right)
-
-  where
-
-  * :math:`k_0` = ``Prefactor``
-    :math:`({\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm MeV}^{-1})`
-  * :math:`\gamma` = ``Index``
-  * :math:`E_0` = ``PivotEnergy``
-    :math:`({\rm MeV})`
-  * :math:`E_{\rm cut}` = ``CutoffEnergy``
-    :math:`({\rm MeV})`
-
-  .. warning::
-
-    The ``PivotEnergy`` parameter is not intended to be fitted.
-
-  .. note::
-
-    For compatibility with the Fermi/LAT ScienceTools the model type
-    ``ExponentialCutoffPowerLaw`` can be replaced by ``ExpCutoff`` and
-    the parameters ``CutoffEnergy`` by ``Cutoff`` and ``PivotEnergy``
-    by ``Scale``.
-
-
-Super exponentially cut-off power law
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-
-  .. code-block:: xml
-
-   <spectrum type="SuperExponentialCutoffPowerLaw">
-    <parameter name="Prefactor"    scale="1e-16" value="1.0" min="1e-07" max="1000.0" free="1"/>
-    <parameter name="Index1"       scale="-1"    value="2.0" min="0.0"   max="+5.0"   free="1"/>
-    <parameter name="CutoffEnergy" scale="1e6"   value="1.0" min="0.01"  max="1000.0" free="1"/>
-    <parameter name="Index2"       scale="1.0"   value="1.5" min="0.1"   max="5.0"    free="1"/>
-    <parameter name="PivotEnergy"  scale="1e6"   value="1.0" min="0.01"  max="1000.0" free="0"/>
-   </spectrum>
-
-  This spectral model component implements the super exponentially cut-off power
-  law function
-
-  .. math::
-     M_{\rm spectral}(E) = k_0 \left( \frac{E}{E_0} \right)^{\gamma}
-                           \exp \left(
-                           -\left( \frac{E}{E_{\rm cut}} \right)^{\alpha}
-                           \right)
-
-  where
-
-  * :math:`k_0` = ``Prefactor``
-    :math:`({\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm MeV}^{-1})`
-  * :math:`\gamma` = ``Index1``
-  * :math:`\alpha` = ``Index2``
-  * :math:`E_0` = ``PivotEnergy``
-    :math:`({\rm MeV})`
-  * :math:`E_{\rm cut}` = ``CutoffEnergy``
-    :math:`({\rm MeV})`
-
-  .. warning::
-
-    The ``PivotEnergy`` parameter is not intended to be fitted.
-
-  An alternative XML format is supported for compatibility with the Fermi/LAT
-  XML format:
-
-  .. code-block:: xml
-
-   <spectrum type="PLSuperExpCutoff">
-    <parameter name="Prefactor"   scale="1e-16" value="1.0" min="1e-07" max="1000.0" free="1"/>
-    <parameter name="Index1"      scale="-1"    value="2.0" min="0.0"   max="+5.0"   free="1"/>
-    <parameter name="Cutoff"      scale="1e6"   value="1.0" min="0.01"  max="1000.0" free="1"/>
-    <parameter name="Index2"      scale="1.0"   value="1.5" min="0.1"   max="5.0"    free="1"/>
-    <parameter name="Scale"       scale="1e6"   value="1.0" min="0.01"  max="1000.0" free="0"/>
-   </spectrum>
-
-
-Broken power law
-^^^^^^^^^^^^^^^^
-
-  .. code-block:: xml
-
-   <spectrum type="BrokenPowerLaw">
-     <parameter name="Prefactor"   scale="1e-16" value="5.7"  min="1e-07" max="1000.0" free="1"/>
-     <parameter name="Index1"      scale="-1"    value="2.48" min="0.0"   max="+5.0"   free="1"/>
-     <parameter name="BreakEnergy" scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="1"/>
-     <parameter name="Index2"      scale="-1"    value="2.70" min="0.01"  max="1000.0" free="1"/>
-   </spectrum>
-
-  This spectral model component implements the broken power law function
-
-  .. math::
-
-     M_{\rm spectral}(E) = k_0 \times \left \{
-     \begin{eqnarray}
-       \left( \frac{E}{E_b} \right)^{\gamma_1} & {\rm if\,\,} E < E_b \\
-       \left( \frac{E}{E_b} \right)^{\gamma_2} & {\rm otherwise}
-     \end{eqnarray}
-     \right .
-
-  where
-
-  * :math:`k_0` = ``Prefactor``
-    :math:`({\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm MeV}^{-1})`
-  * :math:`\gamma_1` = ``Index1``
-  * :math:`\gamma_2` = ``Index2``
-  * :math:`E_b` = ``BreakEnergy``
-    :math:`({\rm MeV})`
-
-  .. warning::
-
-     Note that the ``BreakEnergy`` parameter may be poorly constrained if
-     there is no clear spectral cut-off in the spectrum.
-     This model may lead to complications in the maximum likelihood fitting.
-
-  .. note::
-
-     For compatibility with the Fermi/LAT ScienceTools the parameters
-     ``BreakEnergy`` can be replaced by ``BreakValue``.
-
-
-Smoothly broken power law
-^^^^^^^^^^^^^^^^^^^^^^^^^
-
-  .. code-block:: xml
-
-     <spectrum type="SmoothBrokenPowerLaw">
-       <parameter name="Prefactor"       scale="1e-16" value="5.7"  min="1e-07" max="1000.0" free="1"/>
-       <parameter name="Index1"          scale="-1"    value="2.48" min="0.0"   max="+5.0"   free="1"/>
-       <parameter name="PivotEnergy"     scale="1e6"   value="1.0"  min="0.01"  max="1000.0" free="0"/>
-       <parameter name="Index2"          scale="-1"    value="2.70" min="0.01"  max="+5.0"   free="1"/>
-       <parameter name="BreakEnergy"     scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="1"/>
-       <parameter name="BreakSmoothness" scale="1.0"   value="0.2"  min="0.01"  max="10.0"   free="0"/>
-     </spectrum>
-
-  This spectral model component implements the smoothly broken power law function
-
-  .. math::
-   
-     M_{\rm spectral}(E) = k_0 \left( \frac{E}{E_0} \right)^{\gamma_1}
-                           \left[ 1 +
-                           \left( \frac{E}{E_b} \right)^{\frac{\gamma_1 - \gamma_2}{\beta}}
-                           \right]^{-\beta}
-
-  where
-
-  * :math:`k_0` = ``Prefactor``
-    :math:`({\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm MeV}^{-1})`
-  * :math:`\gamma_1` = ``Index1``
-  * :math:`E_0` = ``PivotEnergy``
-  * :math:`\gamma_2` = ``Index2``
-  * :math:`E_b` = ``BreakEnergy``
-    :math:`({\rm MeV})`
-  * :math:`\beta` = ``BreakSmoothness``
-
-  .. warning::
-
-     The pivot energy should be set far away from the expected break energy
-     value.
-
-  .. warning::
-
-     When the two indices are close together, the :math:`\beta` parameter
-     becomes poorly constrained. Since the :math:`\beta` parameter also scales
-     the indices, this can cause very large errors in the estimates of the
-     various spectral parameters. In this case, consider fixing :math:`\beta`.
-
-  .. note::
-
-     For compatibility with the Fermi/LAT ScienceTools the parameters
-     ``PivotEnergy`` can be replaced by ``Scale``,
-     ``BreakEnergy`` by ``BreakValue`` and
-     ``BreakSmoothness`` by  ``Beta``.
-
-
-Log parabola
-^^^^^^^^^^^^
-
-  .. code-block:: xml
-
-   <spectrum type="LogParabola">
-     <parameter name="Prefactor"   scale="1e-17" value="5.878"   min="1e-07" max="1000.0" free="1"/>
-     <parameter name="Index"       scale="-1"    value="2.32473" min="0.0"   max="+5.0"   free="1"/>
-     <parameter name="Curvature"   scale="-1"    value="0.074"   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>
-
-  This spectral model component implements the log parabola function
-
-  .. math::
-     M_{\rm spectral}(E) = k_0 \left( \frac{E}{E_0} \right)^{\gamma+\eta \ln(E/E_0)}
-
-  where
-
-  * :math:`k_0` = ``Prefactor``
-    :math:`({\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm MeV}^{-1})`
-  * :math:`\gamma` = ``Index``
-  * :math:`\eta` = ``Curvature``
-  * :math:`E_0` = ``PivotEnergy``
-    :math:`({\rm MeV})`
-
-  .. warning::
-
-     The ``PivotEnergy`` parameter is not intended to be fitted.
-
-  An alternative XML format is supported for compatibility with the Fermi/LAT
-  XML format:
-
-  .. code-block:: xml
-
-   <spectrum type="LogParabola">
-     <parameter name="norm"  scale="1e-17" value="5.878"   min="1e-07" max="1000.0" free="1"/>
-     <parameter name="alpha" scale="1"     value="2.32473" min="0.0"   max="+5.0"   free="1"/>
-     <parameter name="beta"  scale="1"     value="0.074"   min="-5.0"  max="+5.0"   free="1"/>
-     <parameter name="Eb"    scale="1e6"   value="1.0"     min="0.01"  max="1000.0" free="0"/>
-   </spectrum>
-
-  where
-
-  * ``alpha`` = -``Index``
-  * ``beta`` = -``Curvature``
-
-
-Gaussian
-^^^^^^^^
-
-  .. code-block:: xml
-
-   <spectrum type="Gaussian">
-     <parameter name="Normalization" scale="1e-10" value="1.0"  min="1e-07" max="1000.0" free="1"/>
-     <parameter name="Mean"          scale="1e6"   value="5.0"  min="0.01"  max="100.0"  free="1"/>
-     <parameter name="Sigma"         scale="1e6"   value="1.0"  min="0.01"  max="100.0"  free="1"/>
-   </spectrum>
-
-  This spectral model component implements the gaussian function
-
-  .. math::
-     M_{\rm spectral}(E) = \frac{N_0}{\sqrt{2\pi}\sigma}
-                           \exp \left( \frac{-(E-\bar{E})^2}{2 \sigma^2} \right)
-
-  where
-
-  * :math:`N_0` = ``Normalization``
-    :math:`({\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1})`
-  * :math:`\bar{E}` = ``Mean``
-    :math:`({\rm MeV})`
-  * :math:`\sigma` = ``Sigma``
-    :math:`({\rm MeV})`
-
-
-File function
-^^^^^^^^^^^^^
-
-  .. code-block:: xml
-
-   <spectrum type="FileFunction" file="data/filefunction.txt">
-     <parameter scale="1.0" name="Normalization" min="0.0" max="1000.0" value="1.0" free="1"/>
-   </spectrum>
-
-  This spectral model component implements an arbitrary function
-  that is defined by intensity values at specific energies.
-  The energy and intensity values are defined using an ASCII file with
-  columns of energy and differential flux values.
-  Energies are given in units of
-  :math:`{\rm MeV}`,
-  intensities are given in units of
-  :math:`{\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm MeV}^{-1}`.
-  The only parameter is a multiplicative normalization:
-
-  .. math::
-     M_{\rm spectral}(E) = N_0 \left. \frac{dN}{dE} \right\rvert_{\rm file}
-
-  where
-
-  * :math:`N_0` = ``Normalization``
-
-  .. warning::
-
-     If the file name is given without a path it is expected that the file
-     resides in the same directory than the XML file.
-     If the file resides in a different directory, an absolute path name should
-     be specified.
-     Any environment variable present in the path name will be expanded.
-
-
-Node function
-^^^^^^^^^^^^^
-
-  .. code-block:: xml
-
-     <spectrum type="NodeFunction">
-       <node>
-         <parameter name="Energy"    scale="1.0"   value="1.0" min="0.1"   max="1.0e20" free="0"/>
-         <parameter name="Intensity" scale="1e-07" value="1.0" min="1e-07" max="1000.0" free="1"/>
-       </node>
-       <node>
-         <parameter name="Energy"    scale="10.0"  value="1.0" min="0.1"   max="1.0e20" free="0"/>
-         <parameter name="Intensity" scale="1e-08" value="1.0" min="1e-07" max="1000.0" free="1"/>
-       </node>
-     </spectrum>
-
-  This spectral model component implements a generalised broken 
-  power law which is defined by a set of energy and intensity values
-  (the so called nodes) that are piecewise connected by power laws.
-  Energies are given in units of
-  :math:`{\rm MeV}`,
-  intensities are given in units of
-  :math:`{\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm MeV}^{-1}`.
-
-  .. warning::
-
-     An arbitrary number of energy-intensity nodes can be defined in a node
-     function.
-     The nodes need to be sorted by increasing energy.
-     Although the fitting of the ``Energy`` parameters is formally possible
-     it may lead to numerical complications.
-     If ``Energy`` parameters are to be fitted make sure that the ``min``
-     and ``max`` attributes are set in a way that avoids inversion of the energy
-     ordering.
-
-
-Composite model
-^^^^^^^^^^^^^^^
-
-  .. code-block:: xml
-
-     <spectrum type="Composite">
-       <spectrum type="PowerLaw" component="SoftComponent">
-         <parameter name="Prefactor"   scale="1e-17" value="3"  min="1e-07" max="1000.0" free="1"/>
-         <parameter name="Index"       scale="-1"    value="3.5" min="0.0"   max="+5.0"   free="1"/>
-         <parameter name="PivotEnergy" scale="1e6"   value="1"  min="0.01"  max="1000.0" free="0"/>
-       </spectrum>
-       <spectrum type="PowerLaw" component="HardComponent">
-         <parameter name="Prefactor"   scale="1e-17" value="5"  min="1e-07" max="1000.0" free="1"/>
-         <parameter name="Index"       scale="-1"    value="2.0" min="0.0"   max="+5.0"   free="1"/>
-         <parameter name="PivotEnergy" scale="1e6"   value="1"  min="0.01"  max="1000.0" free="0"/>
-       </spectrum>
-     </spectrum>
-
-  This spectral model component implements a composite model that is the
-  sum of an arbitrary number of spectral models, computed using
-
-  .. math::
-     M_{\rm spectral}(E) = \sum_{i=0}^{N-1} M_{\rm spectral}^{(i)}(E)
-
-  where :math:`M_{\rm spectral}^{(i)}(E)` is any spectral model component
-  (including another composite model), and :math:`N` is the number of
-  model components that are combined.
-
-
-Multiplicative model
-^^^^^^^^^^^^^^^^^^^^
-
-  .. code-block:: xml
-
-     <spectrum type="Multiplicative">
-       <spectrum type="PowerLaw" component="PowerLawComponent">
-         <parameter name="Prefactor"   scale="1e-17" value="1.0"  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="1.0"  min="0.01"  max="1000.0" free="0"/>
-       </spectrum>
-       <spectrum type="ExponentialCutoffPowerLaw" component="CutoffComponent">
-         <parameter name="Prefactor"    scale="1.0" value="1.0" min="1e-07" max="1000.0" free="0"/>
-         <parameter name="Index"        scale="1.0" value="0.0" min="-2.0"  max="+2.0"   free="0"/>
-         <parameter name="CutoffEnergy" scale="1e6" value="1.0" min="0.01"  max="1000.0" free="1"/>
-         <parameter name="PivotEnergy"  scale="1e6" value="1.0" min="0.01"  max="1000.0" free="0"/>
-       </spectrum>
-     </spectrum>
-
-  This spectral model component implements a composite model that is the
-  product of an arbitrary number of spectral models, computed using
-
-  .. math::
-     M_{\rm spectral}(E) = \prod_{i=0}^{N-1} M_{\rm spectral}^{(i)}(E)
-
-  where :math:`M_{\rm spectral}^{(i)}(E)` is any spectral model component
-  (including another composite model), and :math:`N` is the number of
-  model components that are multiplied.
-
-
-.. _sec_temporal_models:
-
-Temporal model components
-~~~~~~~~~~~~~~~~~~~~~~~~~
-
-The following sections present the temporal model components that are available 
-in ctools.
-
-Constant
-^^^^^^^^
-
-  .. code-block:: xml
-
-     <temporal type="Constant">
-       <parameter name="Normalization" scale="1.0" value="1.0" min="0.1" max="10.0" free="0"/>
-     </temporal>
-
-  This temporal model component implements a constant source
-
-  .. math::
-
-     M_{\rm temporal}(t) = N_0
-
-  where
-
-  * :math:`N_0` = ``Normalization``
-
-
-Light Curve
-^^^^^^^^^^^
-
-  .. code-block:: xml
-
-     <temporal type="LightCurve" file="model_temporal_lightcurve.fits">
-       <parameter name="Normalization" scale="1" value="1.0" min="0.0" max="1000.0" free="0"/>
-     </temporal>
-
-  This temporal model component implements a light curve :math:`r(t)`
-
-  .. math::
-     M_{\rm temporal}(t) = N_0 \times r(t)
-
-  where
-
-  * :math:`N_0` = ``Normalization``
-
-  The light curve is defined by nodes in a FITS file that specify the relative
-  flux normalization as function of time (file ``model_temporal_lightcurve.fits``
-  in the example above). The structure of the light curve FITS
-  file is shown in the figure below. The light curve is defined in the first
-  extension of the FITS file and consists of a binary table with the columns
-  ``TIME`` and ``NORM``. Times in the ``TIME`` columns are given in seconds
-  and are counted with respect to a time reference that is defined in the
-  header of the binary table. Times need to be specified in ascending order.
-  The values in the ``NORM`` column specify :math:`r(t)` at times :math:`t`,
-  and should be comprised between 0 and 1.
-
-  .. _fig_model_lightcurve:
-
-  .. figure:: model_lightcurve.png
-     :align: center
-     :width: 100%
-
-     *Structure of light curve FITS file*
-
-  .. warning::
-     Fitting of light curves only makes sense for an unbinned maximum likelihood
-     analysis, since in a binned or stacked analysis the times of individual
-     events are dropped.
-
-
-Phase Curve
-^^^^^^^^^^^
-
-  .. code-block:: xml
-
-     <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>
-
-  This temporal model component implements a phase curve :math:`r(\Phi(t))`
-
-  .. math::
-     M_{\rm temporal}(t) = N_0 \times r(\Phi(t))
-
-  where the phase as function of time is computed using
-
-  .. math::
-     \Phi(t) = \Phi_0 + f(t-t_0) + \frac{1}{2}\dot{f} (t-t_0)^2 +
-                                   \frac{1}{6}\ddot{f} (t-t_0)^3
-
-  and
-
-  * :math:`N_0` = ``Normalization``
-  * :math:`t_0` = ``MJD``
-  * :math:`\Phi_0` = ``Phase``
-  * :math:`f` = ``F0``
-  * :math:`\dot{f}` = ``F1``
-  * :math:`\ddot{f}` = ``F2``
-
-  The phase curve is defined by nodes in a FITS file that specify the relative
-  flux normalization as function of phase (file ``model_temporal_phasecurve.fits``
-  in the example above). The structure of the phase curve
-  FITS file is shown in the figure below. The phase curve is defined in the
-  first extension of the FITS file and consists of a binary table with the
-  columns ``PHASE`` and ``NORM``. Phase values in the ``PHASE`` column need to
-  be comprised between 0 and 1 and need to be given in ascending order. The
-  values in the ``NORM`` column specify :math:`r(\Phi(t))` at phases
-  :math:`\Phi(t)`, and should be comprised between 0 and 1.
-
-  .. _fig_model_phasecurve:
-
-  .. figure:: model_phasecurve.png
-     :align: center
-     :width: 40%
-
-     *Structure of phase curve FITS file*
-
-  .. warning::
-     Fitting of phase curves only makes sense for an unbinned maximum likelihood
-     analysis, since in a binned or stacked analysis the times of individual
-     events are dropped.
-
-  .. warning::
-     Fitting of phase curve parameters may not properly work for pulsar
-     frequencies.
+.. toctree::
+   :maxdepth: 1
+
+   models_concept
+   models_implementation
+   models_spatial_src
+   models_spatial_bgd
+   models_spectral
+   models_temporal
+.. _sec_models_concept:
+
+Overall concept
+------------------
+
+A central element of any data analysis is the modeling of the observed
+event distribution.
+Generally, the measured events can be separated into two distinct classes:
+events being attributed to gamma rays of celestial origin (the "source" 
+events) or events being attributed to any other kind of trigger (the 
+"background" events).
+In the first case, the :ref:`Instrument Response Functions (IRFs) <sec_response>`
+describe how an incident gamma-ray converts into a measured event.
+In the second case there is often no general prescription, and the 
+distribution of background events is commonly modeled directly in the data 
+space of the observable quantities.
+The main difficulty of gamma-ray astronomy is that not all events can be
+tagged a priori as being source or background, hence their separation can 
+only be done on a statistical basis, based on the different morphologies, 
+spectral characteristics or eventually temporal signatures of both event
+categories.
+
+For this purpose, ctools use a general model that describes the spatial, 
+spectral and temporal properties of the source and background components.
+The model is composed of an arbitrary number of model components that
+add up linearly to provide a prediction of the expected number of events.
+Model components are generally parametric, and model parameters can be 
+adjusted through a maximum likelihood procedure to find the set of 
+parameters that represent best the measured data.
+Model components representing celestial sources are convolved with the 
+:ref:`IRFs <sec_response>` to predict the expected number of source events in 
+the data.
+Background model components will be directly expressed as expected number 
+of background events without any :ref:`IRF <sec_response>` convolution.
+.. _sec_models_implementation:
+
+Implementation
+--------------
+
+The general model is describe in ctools using a model definition XML file. 
+Below is a simple example of such a file comprising one source and one 
+background model.
+Each model is factorised into
+a spectral (tag ``<spectrum>``),
+a spatial (tags ``<spatialModel>`` and ``<radialModel>``), and
+a temporal component (tag ``<temporal>``).
+
+.. math::
+  M(x,y,E,t) = M_{\rm spatial}(x,y|E) \times M_{\rm spectral}(E) \times M_{\rm temporal}(t)
+
+In this specific example, the source component ``Crab`` describes 
+a point source at the location of the Crab nebula with a power law spectral
+shape that is constant in time.
+The background component ``Background`` is modelled using the template that is
+provided with the
+:ref:`Instrument Response Functions (IRFs) <sec_response>`
+and that is multiplied by a spectral power law function.
+
+.. code-block:: xml
+
+  <?xml version="1.0" standalone="no"?>
+  <source_library title="source library">
+    <source name="Crab" type="PointSource">
+      <spatialModel type="PointSource">
+        <parameter name="RA"  scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
+        <parameter name="DEC" scale="1.0" value="22.0145" min="-90"  max="90"  free="1"/>
+      </spatialModel>
+      <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>
+      <temporal type="Constant">
+        <parameter name="Normalization" scale="1.0" value="1.0" min="0.0" max="1000.0" free="0"/>
+      </temporal>
+    </source>
+    <source name="CTABackgroundModel" type="CTAIrfBackground" 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>
+
+Model parameters are specified by a ``<parameter>`` tag with a certain 
+number of attributes.
+The ``name`` attribute gives the (case-sensitive) parameter name that 
+needs to be unique within a given model component.
+The ``scale`` attribute gives a scaling factor that will be multiplied by 
+the ``value`` attribute to provide the real (physical) parameter value.
+The ``min`` and ``max`` attributes specify boundaries for the ``value``
+term of the parameter.
+And the ``free`` attribute specifies whether a parameter should be fitted 
+(``free="1"``) or kept fixed (``free="0"``) in a maximum likelihood 
+analysis.
+After a maximum likelihood fit using :ref:`ctlike`, an
+``error`` attribute giving the statisical uncertainty of the ``value``
+term will be appended to each ``<parameter>`` tag.
+
+.. note::
+   For compatibility reasons with the Fermi/LAT XML format the ``<temporal>``
+   tag can be omitted for models that are constant in time.
+
+.. note::
+
+   XML files are ASCII files and can be edited by hand using any text 
+   editor.
+   The indentation of the XML elements is not mandatory.
+
+.. note::
+
+   The splitting of parameter values into a ``value`` and ``scale`` factor 
+   is mainly for numerical purposes.
+   Parameter fitting algorithms can be ill-conditioned if several 
+   parameters of vastly different orders of magnitudes need to be 
+   optimised simultaneously.
+   Splitting a value into two components allows to "prescale" the 
+   parameters so that the effective parameters to be optimised (the ``value`` terms) 
+   are all of about unity.
+
+.. note::
+
+   The syntax of the model definition XML file has been inspired from the
+   syntax used by the Fermi/LAT ScienceTools, but for reasons of clarity and
+   homogenity of the various model and parameter names we have made some
+   modifications.
+   Nevertheless, the format used by the Fermi/LAT ScienceTools is also
+   supported.
+.. _sec_models_spatial_bgd:
+
+Spatial background model components
+-----------------------------------
+
+The following sections present the spatial model components that are available 
+in ctools for instrumental background modelling for CTA.
+
+CTA radial background
+^^^^^^^^^^^^^^^^^^^^^
+
+  There exist a number of radial CTA background models that factorise into a
+  spatial and spectral component, and that model the spatial component as a
+  function that only dependends on the offset angle :math:`\theta` from the
+  pointing direction. These models have the type ``RadialAcceptance`` and
+  require a ``<radialModel>`` tag as the spatial component.
+
+  The first radial model is of type ``Gaussian``
+
+  .. code-block:: xml
+
+     <source name="Background" type="RadialAcceptance" instrument="CTA">
+       <radialModel type="Gaussian">
+         <parameter name="Sigma" scale="1.0" value="3.0" min="0.01" max="10.0" free="1"/>
+       </radialModel>
+       <spectrum type="PowerLaw">
+         <parameter name="Prefactor"   scale="1e-6" value="61.8" min="0.0"  max="1000.0" free="1"/>
+         <parameter name="Index"       scale="-1"   value="1.85" min="0.0"  max="+5.0"   free="1"/>
+         <parameter name="PivotEnergy" scale="1e6"  value="1.0"  min="0.01" max="1000.0" free="0"/>
+       </spectrum>
+     </source>
+
+  and implements
+
+  .. math::
+     M_{\rm spatial}(\theta) = \exp \left(-\frac{1}{2}
+                               \left( \frac{\theta^2}{\sigma} \right)^2 \right)
+
+  where
+
+  * :math:`\sigma` = ``Sigma`` (degrees)
+
+  The second radial model is of type ``Profile``
+
+  .. code-block:: xml
+
+     <source name="Background" type="RadialAcceptance" instrument="CTA">
+       <radialModel type="Profile">
+         <parameter name="Width" scale="1.0" value="1.5" min="0.1" max="1000.0" free="1"/>
+         <parameter name="Core"  scale="1.0" value="3.0" min="0.1" max="1000.0" free="1"/>
+         <parameter name="Tail"  scale="1.0" value="5.0" min="0.1" max="1000.0" free="1"/>
+       </radialModel>
+        <spectrum type="PowerLaw">
+          <parameter name="Prefactor"   scale="1e-6" value="61.8" min="0.0"  max="1000.0" free="1"/>
+          <parameter name="Index"       scale="-1"   value="1.85" min="0.0"  max="+5.0"   free="1"/>
+          <parameter name="PivotEnergy" scale="1e6"  value="1.0"  min="0.01" max="1000.0" free="0"/>
+        </spectrum>
+     </source>
+
+  and implements
+
+  .. math::
+     M_{\rm spatial}(\theta) = (1 + (\theta/c_0)^{c_1})^{-c_2/c_1}
+
+  where
+
+  * :math:`c_0` = ``Width`` (degrees)
+  * :math:`c_1` = ``Core``
+  * :math:`c_2` = ``Tail``
+
+  The third radial model is of type ``Polynom``
+
+  .. code-block:: xml
+
+     <source name="Background" type="RadialAcceptance" instrument="CTA">
+       <radialModel type="Polynom">
+         <parameter name="Coeff0" scale="1.0" value="+1.00000"   min="-10.0" max="10.0" free="0"/>
+         <parameter name="Coeff1" scale="1.0" value="-0.1239176" min="-10.0" max="10.0" free="1"/>
+         <parameter name="Coeff2" scale="1.0" value="+0.9751791" min="-10.0" max="10.0" free="1"/>
+         <parameter name="Coeff3" scale="1.0" value="-3.0584577" min="-10.0" max="10.0" free="1"/>
+         <parameter name="Coeff4" scale="1.0" value="+2.9089535" min="-10.0" max="10.0" free="1"/>
+         <parameter name="Coeff5" scale="1.0" value="-1.3535372" min="-10.0" max="10.0" free="1"/>
+         <parameter name="Coeff6" scale="1.0" value="+0.3413752" min="-10.0" max="10.0" free="1"/>
+         <parameter name="Coeff7" scale="1.0" value="-0.0449642" min="-10.0" max="10.0" free="1"/>
+         <parameter name="Coeff8" scale="1.0" value="+0.0024321" min="-10.0" max="10.0" free="1"/>
+       </radialModel>
+       <spectrum type="PowerLaw">
+         <parameter name="Prefactor"   scale="1e-6" value="61.8" min="0.0"  max="1000.0" free="1"/>
+         <parameter name="Index"       scale="-1"   value="1.85" min="0.0"  max="+5.0"   free="1"/>
+         <parameter name="PivotEnergy" scale="1e6"  value="1.0"  min="0.01" max="1000.0" free="0"/>
+       </spectrum>
+     </source>
+
+  and implements
+
+  .. math::
+     M_{\rm spatial}(\theta) = \sum_{i=0}^m c_i \theta^i
+
+  where
+
+  * :math:`c_0` = ``Coeff0``
+  * :math:`c_1` = ``Coeff1``
+  * ...
+
+  (the number of polynomial coefficients is arbitrary).
+
+
+CTA IRF background
+^^^^^^^^^^^^^^^^^^
+
+  The :ref:`Instrument Response Functions (IRFs) <sec_response>` contain a template
+  that predicts the background rate as function of position in the field of view
+  and measured energy in units of
+  :math:`{\rm events} \, {\rm s}^{-1} {\rm MeV}^{-1} {\rm sr}^{-1}`. This template
+  can be used by specifying a model of type ``CTAIrfBackground``. No spatial component
+  will be specified explicitly since the spatial (and spectral) information is
+  already contained in the template. The model will be multiplied by a spectral law.
+  In the example below, the template is multiplied by a power law with normalization
+  of 1 and slope 0 (i.e. the template is taken as is, but a fit can adjust the model
+  to compensate for inaccuracies):
+
+  .. code-block:: xml
+
+     <source name="Background" type="CTAIrfBackground" instrument="CTA">
+       <spectrum type="PowerLaw">
+         <parameter name="Prefactor"   scale="1.0" value="1.0" min="1e-3" max="1e3"    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>
+
+  If the observation is an On/Off observation, do not forget to switch the instrument
+  to ``CTAOnOff``:
+
+  .. code-block:: xml
+
+     <source name="Background" type="CTAIrfBackground" instrument="CTAOnOff">
+       <spectrum type="PowerLaw">
+         <parameter name="Prefactor"   scale="1.0" value="1.0" min="1e-3" max="1e3"    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>
+
+
+CTA effective area background
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+  Instead of using the background template the effective area for gamma rays can
+  also be used to model the instrumental background. Note that in this case the
+  effective area has to be scaled to a reasonable background rate. An example is
+  given below:
+
+  .. code-block:: xml
+
+     <source name="Background" type="CTAAeffBackground" instrument="CTA">
+       <spectrum type="PowerLaw">
+         <parameter name="Prefactor"   scale="1e-14" value="1.0"  min="1e-3" max="1e3"    free="1"/>
+         <parameter name="Index"       scale="1.0"   value="-2.4" 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>
+
+
+CTA cube background
+^^^^^^^^^^^^^^^^^^^
+
+  For a stacked analysis, the background rate is predicted by a so called background
+  cube. The background cube is used by specifying a model of type ``CTACubeBackground``.
+  Similar to the ``CTAIrfBackground`` model, the background cube is multplied by
+  a spectral model.
+
+  .. code-block:: xml
+
+     <source name="Background" type="CTACubeBackground" instrument="CTA">
+       <spectrum type="PowerLaw">
+         <parameter name="Prefactor"   scale="1.0" value="1.0" min="1e-3" max="1e3"    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>
+.. _sec_models_spatial_src:
+
+Spatial source model components
+-------------------------------
+
+The following sections present the spatial model components that are available 
+in ctools for gamma-ray sources.
+
+.. note::
+   Except of the ``DiffuseMapCube`` model, all spatial models are normalized
+   so that when integrated over the sphere the result is unity. For clarify,
+   the spatial normalisation is omitted in the formulae below.
+
+Point source
+^^^^^^^^^^^^
+
+  The ``PointSource`` model describes a point source
+
+  .. code-block:: xml
+
+    <source name="Crab" type="PointSource">
+      <spatialModel type="PointSource">
+        <parameter name="RA"  scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
+        <parameter name="DEC" scale="1.0" value="22.0145" min="-90"  max="90"  free="1"/>
+      </spatialModel>
+      <spectrum type="...">
+        ...
+      </spectrum>
+    </source>
+
+  where
+
+  * ``RA`` is the Right Ascension (degrees)
+  * ``DEC`` is the Declination (degrees)
+
+  The ``RA`` and ``DEC`` parameters can be remplaced by ``GLON`` and ``GLAT``
+  for Galactic coordinates
+
+  .. code-block:: xml
+
+    <source name="Crab" type="PointSource">
+      <spatialModel type="PointSource">
+        <parameter name="GLON" scale="1.0" value="184.5575" min="-360" max="360" free="1"/>
+        <parameter name="GLAT" scale="1.0" value="-5.7844"  min="-90"  max="90"  free="1"/>
+      </spatialModel>
+      <spectrum type="...">
+        ...
+      </spectrum>
+    </source>
+
+  where
+
+  * ``GLON`` is the Galactic longitude (degrees)
+  * ``GLAT`` is the Galactic latitude (degrees)
+
+  .. note::
+     For compatibility with the Fermi/LAT ScienceTools the model type
+     ``PointSource`` can be replaced by ``SkyDirFunction``.
+
+
+Radial source
+^^^^^^^^^^^^^
+
+  The ``RadialDisk`` model describes a uniform intensity distribution within
+  a given radius
+
+  .. code-block:: xml
+
+     <source name="Crab" type="ExtendedSource">
+       <spatialModel type="RadialDisk">
+         <parameter name="RA"     scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
+         <parameter name="DEC"    scale="1.0" value="22.0145" min="-90"  max="90"  free="1"/>
+         <parameter name="Radius" scale="1.0" value="0.20"    min="0.01" max="10"  free="1"/>
+       </spatialModel>
+       <spectrum type="...">
+         ...
+       </spectrum>
+     </source>
+
+  where
+
+  * ``RA`` is the Right Ascension (degrees)
+  * ``DEC`` is the Declination (degrees)
+  * ``Radius`` is the disk radius (degrees)
+
+  The ``RadialGaussian`` model describes a Gaussian intensity distribution
+
+  .. code-block:: xml
+
+     <source name="Crab" type="ExtendedSource">
+       <spatialModel type="RadialGaussian">
+         <parameter name="RA"    scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
+         <parameter name="DEC"   scale="1.0" value="22.0145" min="-90"  max="90"  free="1"/>
+         <parameter name="Sigma" scale="1.0" value="0.20"    min="0.01" max="10"  free="1"/>
+       </spatialModel>
+       <spectrum type="...">
+         ...
+       </spectrum>
+     </source>
+
+  and implements
+
+  .. math::
+     M_{\rm spatial}(\theta) = \frac{1}{2 \pi \sigma^2} \exp
+                    \left(-\frac{1}{2}\frac{\theta^2}{\sigma^2} \right),
+
+  where
+
+  * ``RA`` is the Right Ascension (degrees)
+  * ``DEC`` is the Declination (degrees)
+  * :math:`\sigma` == ``Sigma`` (degrees)
+
+  The ``RadialShell`` model describes a spherical shell projected on the sky
+
+  .. code-block:: xml
+
+     <source name="Crab" type="ExtendedSource">
+       <spatialModel type="RadialShell">
+         <parameter name="RA"     scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
+         <parameter name="DEC"    scale="1.0" value="22.0145" min="-90"  max="90"  free="1"/>
+         <parameter name="Radius" scale="1.0" value="0.30"    min="0.01" max="10"  free="1"/>
+         <parameter name="Width"  scale="1.0" value="0.10"    min="0.01" max="10"  free="1"/>
+       </spatialModel>
+       <spectrum type="...">
+         ...
+       </spectrum>
+     </source>
+
+  and implements
+
+  .. math::
+     M_{\rm spatial}(\theta) =  n_0 \left \{
+     \begin{array}{l l}
+        \displaystyle
+        \sqrt{ \theta_{\rm out}^2 - \theta^2 } - \sqrt{ \theta_{\rm in}^2 - \theta^2 }
+        & \mbox{if $\theta \le \theta_{\rm in}$} \\
+        \\
+       \displaystyle
+        \sqrt{ \theta_{\rm out}^2 - \theta^2 }
+        & \mbox{if $\theta_{\rm in} < \theta \le \theta_{\rm out}$} \\
+        \\
+       \displaystyle
+       0 & \mbox{if $\theta > \theta_{\rm out}$}
+     \end{array}
+     \right .
+
+  where
+
+  * ``RA`` is the Right Ascension (degrees)
+  * ``DEC`` is the Declination (degrees)
+  * :math:`\theta_{\rm out}` = ``Radius`` + ``Width`` (degrees)
+  * :math:`\theta_{\rm in}` = ``Radius`` (degrees)
+
+  .. note::
+     For all radial source models the ``RA`` and ``DEC`` parameters can be
+     remplaced by ``GLON`` and ``GLAT`` for Galactic coordinates.
+
+
+Elliptical source
+^^^^^^^^^^^^^^^^^
+
+  The ``EllipticalDisk`` model describes a uniform intensity distribution within
+  an elliptical circumference:
+
+  .. code-block:: xml
+
+     <source name="Crab" type="ExtendedSource">
+       <spatialModel type="EllipticalDisk">
+         <parameter name="RA"          scale="1.0" value="83.6331" min="-360"  max="360" free="1"/>
+         <parameter name="DEC"         scale="1.0" value="22.0145" min="-90"   max="90"  free="1"/>
+         <parameter name="PA"          scale="1.0" value="45.0"    min="-360"  max="360" free="1"/>
+         <parameter name="MinorRadius" scale="1.0" value="0.5"     min="0.001" max="10"  free="1"/>
+         <parameter name="MajorRadius" scale="1.0" value="2.0"     min="0.001" max="10"  free="1"/>
+       </spatialModel>
+       <spectrum type="...">
+         ...
+       </spectrum>
+     </source>
+
+  where
+
+  * ``RA`` is the Right Ascension (degrees)
+  * ``DEC`` is the Declination (degrees)
+  * ``PA`` is the position angle, counted counterclockwise from North (degrees)
+  * ``MinorRadius`` is the minor radius of the ellipse (degrees)
+  * ``MajorRadius`` is the major radius of the ellipse (degrees)
+
+  The ``EllipticalGaussian`` model describes a Gaussian intensity distribution
+
+  .. code-block:: xml
+
+    <source name="Crab" type="ExtendedSource">
+      <spatialModel type="EllipticalGaussian">
+        <parameter name="RA"          scale="1.0" value="83.6331" min="-360"  max="360" free="1"/>
+        <parameter name="DEC"         scale="1.0" value="22.0145" min="-90"   max="90"  free="1"/>
+        <parameter name="PA"          scale="1.0" value="45.0"    min="-360"  max="360" free="1"/>
+        <parameter name="MinorRadius" scale="1.0" value="0.5"     min="0.001" max="10"  free="1"/>
+        <parameter name="MajorRadius" scale="1.0" value="2.0"     min="0.001" max="10"  free="1"/>
+      </spatialModel>
+      <spectrum type="...">
+        ...
+      </spectrum>
+    </source>
+
+  and implements
+
+  .. math::
+     M_{\rm spatial}(\theta, \phi) = \exp \left( -\frac{\theta^2}{2 r_\mathrm{eff}^2} \right),
+
+  with
+
+  .. math::
+     r_\mathrm{eff} = \frac{ab} {\sqrt{\left( a \sin (\phi - \phi_0) \right)^2 +
+                      \sqrt{\left( b \cos (\phi - \phi_0) \right)^2}}}
+
+  where
+
+  * ``RA`` is the Right Ascension (degrees)
+  * ``DEC`` is the Declination (degrees)
+  * ``PA`` is the position angle, counted counterclockwise from North (degrees)
+  * :math:`a` = ``MinorRadius`` (degrees)
+  * :math:`b` = ``MajorRadius`` (degrees)
+  * :math:`\phi_0` is the position angle of the ellipse, counted counterclockwise
+    from North
+  * :math:`\phi` is the azimuth angle with respect to North.
+
+  .. note::
+     For all elliptical source models the ``RA`` and ``DEC`` parameters can be
+     remplaced by ``GLON`` and ``GLAT`` for Galactic coordinates.
+
+
+Diffuse source
+^^^^^^^^^^^^^^
+
+  The ``DiffuseIsotropic`` model describes an isotropic intensity distribution
+
+  .. code-block:: xml
+
+     <source name="Crab" type="DiffuseSource">
+       <spatialModel type="DiffuseIsotropic">
+         <parameter name="Value" scale="1" value="1" min="1"  max="1" free="0"/>
+       </spatialModel>
+       <spectrum type="...">
+         ...
+       </spectrum>
+     </source>
+
+  where
+
+  * ``Value`` is isotropic intensity
+
+  .. note::
+     For compatibility with the Fermi/LAT ScienceTools the model type
+     ``DiffuseIsotropic`` can be replaced by ``ConstantValue``.
+
+  The ``DiffuseMap`` model describes an arbitrary intensity distribution in
+  form of a sky map
+
+  .. code-block:: xml
+
+     <source name="Crab" type="DiffuseSource">
+       <spatialModel type="DiffuseMap" file="map.fits">
+         <parameter name="Normalization" scale="1" value="1" min="0.001" max="1000.0" free="0"/>
+       </spatialModel>
+       <spectrum type="...">
+         ...
+       </spectrum>
+     </source>
+
+  where
+
+  * ``Normalization`` is a normalization value
+
+  and the ``file`` attribute specifies a sky map FITS file name. If a file name
+  without path is specified it is assumed that the FITS file resides in the same
+  directory as the model definition XML file.
+
+  .. note::
+     For compatibility with the Fermi/LAT ScienceTools the model type
+     ``DiffuseMap`` can be replaced by ``SpatialMap`` and the parameter
+     ``Normalization`` can be replaced by ``Prefactor``.
+
+  .. note::
+     You may add the ``normalize="0"``to the spatial model tag to avoid the
+     normalization of the diffuse map to an integral of unity
+
+     .. code-block:: xml
+
+        <source name="Crab" type="DiffuseSource">
+          <spatialModel type="DiffuseMap" file="map.fits" normalize="0">
+            <parameter name="Normalization" scale="1" value="1" min="0.001" max="1000.0" free="0"/>
+          </spatialModel>
+          <spectrum type="...">
+            ...
+          </spectrum>
+        </source>
+
+  The ``DiffuseMapCube`` model describes an arbitrary energy-dependent intensity
+  distribution in form of a map cube
+
+  .. code-block:: xml
+
+     <source name="Crab" type="DiffuseSource">
+       <spatialModel type="DiffuseMapCube" file="map_cube.fits">
+         <parameter name="Normalization" scale="1" value="1" min="0.001" max="1000.0" free="0"/>
+       </spatialModel>
+       <spectrum type="...">
+         ...
+       </spectrum>
+     </source>
+
+  where
+
+  * ``Normalization`` is a normalization value
+
+  .. note::
+     For compatibility with the Fermi/LAT ScienceTools the model type
+     ``DiffuseMapCube`` can be replaced by ``MapCubeFunction`` and the parameter
+     ``Normalization`` can be replaced by ``Value``.
+
+
+Composite model
+^^^^^^^^^^^^^^^
+
+  The ``Composite`` model implements a composite model that is the sum of an
+  arbitrary number of spatial models
+
+  .. code-block:: xml
+
+     <source name="Crab" type="CompositeSource">
+       <spatialModel type="Composite">
+         <spatialModel type="PointSource" component="PointSource">
+           <parameter name="RA"    scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
+           <parameter name="DEC"   scale="1.0" value="22.0145" min="-90"  max="90"  free="1"/>
+         </spatialModel>
+         <spatialModel type="RadialGaussian">
+           <parameter name="RA"    scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
+           <parameter name="DEC"   scale="1.0" value="22.0145" min="-90"  max="90"  free="1"/>
+           <parameter name="Sigma" scale="1.0" value="0.20"    min="0.01" max="10"  free="1"/>
+         </spatialModel>
+       </spatialModel>
+       <spectrum type="...">
+         ...
+       </spectrum>
+     </source>
+
+  which implements
+
+  .. math::
+     M_{\rm spatial}(x,y|E) = \frac{1}{N} \sum_{i=0}^{N-1} M_{\rm spatial}^{(i)}(x,y|E)
+
+  where :math:`M_{\rm spatial}^{(i)}(x,y|E)` is any spatial model component
+  (including another composite model), and :math:`N` is the number of
+  model components that are combined.
+.. _sec_models_spectral:
+
+Spectral model components
+-------------------------
+
+The following sections present the spectral model components that are available 
+in ctools.
+
+.. warning::
+
+   Source intensities are generally given in units of
+   :math:`{\rm photons}\,\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm MeV}^{-1}`.
+
+   **An exception to this rule exists for the** ``DiffuseMapCube`` **spatial
+   model where intensities are unitless** and the spectral model presents a
+   relative scaling of the diffuse model cube values.
+
+   If spectral models are used in combination with ``RadialAcceptance`` CTA
+   background models, intensity units are given in
+   :math:`{\rm events}\,\,{\rm s}^{-1}\,{\rm MeV}^{-1}\,{\rm sr}^{-1}`
+   and correspond to the on-axis count rate.
+
+   For other CTA background models intensities are unitless and the spectral
+   model presents a relative scaling of the background model values.
+
+
+Constant
+^^^^^^^^
+
+  .. code-block:: xml
+
+     <spectrum type="Constant">
+       <parameter name="Normalization" scale="1e-16" value="5.7" min="1e-07" max="1000.0" free="1"/>
+     </spectrum>
+
+  This spectral model component implements the constant function
+
+  .. math::
+     M_{\rm spectral}(E) = N_0
+
+  where
+
+  * :math:`N_0` = ``Normalization``
+    :math:`({\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm MeV}^{-1})`
+
+  .. note::
+     For compatibility with the Fermi/LAT ScienceTools the model type
+     ``Constant`` can be replaced by ``ConstantValue`` and the parameter
+     ``Normalization`` by ``Value``.
+
+
+Power law
+^^^^^^^^^
+
+  .. code-block:: xml
+
+    <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>
+
+  This spectral model component implements the power law function
+
+  .. math::
+     M_{\rm spectral}(E) = k_0 \left( \frac{E}{E_0} \right)^{\gamma}
+
+  where
+
+  * :math:`k_0` = ``Prefactor``
+    :math:`({\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm MeV}^{-1})`
+  * :math:`\gamma` = ``Index``
+  * :math:`E_0` = ``PivotEnergy``
+    :math:`({\rm MeV})`
+
+  .. warning::
+     The ``PivotEnergy`` parameter is not intended to be fitted.
+
+  .. note::
+     For compatibility with the Fermi/LAT ScienceTools the parameter
+     ``PivotEnergy`` can be replaced by ``Scale``.
+
+  An alternative power law function that uses the integral photon flux as
+  parameter rather than the Prefactor is specified by
+
+  .. code-block:: xml
+
+    <spectrum type="PowerLaw">
+      <parameter scale="1e-07" name="PhotonFlux" min="1e-07" max="1000.0"    value="1.0" free="1"/>
+      <parameter scale="1.0"   name="Index"      min="-5.0"  max="+5.0"      value="-2.0" free="1"/>
+      <parameter scale="1.0"   name="LowerLimit" min="10.0"  max="1000000.0" value="100.0" free="0"/>
+      <parameter scale="1.0"   name="UpperLimit" min="10.0"  max="1000000.0" value="500000.0" free="0"/>
+    </spectrum>
+
+  This spectral model component implements the power law function
+
+  .. math::
+     M_{\rm spectral}(E) = \frac{N(\gamma+1)E^{\gamma}}
+                                {E_{\rm max}^{\gamma+1} - E_{\rm min}^{\gamma+1}}
+
+  where
+
+  * :math:`N` = ``PhotonFlux``
+    :math:`({\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1})`
+  * :math:`\gamma` = ``Index``
+  * :math:`E_{\rm min}` = ``LowerLimit``
+    :math:`({\rm MeV})`
+  * :math:`E_{\rm max}` = ``UpperLimit``
+    :math:`({\rm MeV})`
+
+  .. warning::
+     The ``LowerLimit`` and ``UpperLimit`` parameters are always treated as fixed
+     and the flux given by the ``PhotonFlux`` parameter is computed over the
+     range set by these two parameters.
+     Use of this model allows the errors on the integral flux to be evaluated directly
+     by :ref:`ctlike`.
+
+  .. note::
+     For compatibility with the Fermi/LAT ScienceTools the model type
+     ``PowerLaw`` can be replaced by ``PowerLaw2`` and the parameter
+     ``PhotonFlux`` by ``Integral``.
+
+
+Exponentially cut-off power law
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+  .. code-block:: xml
+
+    <spectrum type="ExponentialCutoffPowerLaw">
+      <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="CutoffEnergy" scale="1e6"   value="1.0"  min="0.01"  max="1000.0" free="1"/>
+      <parameter name="PivotEnergy"  scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="0"/>
+    </spectrum>
+
+  This spectral model component implements the exponentially cut-off power law
+  function
+
+  .. math::
+     M_{\rm spectral}(E) = k_0 \left( \frac{E}{E_0} \right)^{\gamma}
+                           \exp \left( \frac{-E}{E_{\rm cut}} \right)
+
+  where
+
+  * :math:`k_0` = ``Prefactor``
+    :math:`({\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm MeV}^{-1})`
+  * :math:`\gamma` = ``Index``
+  * :math:`E_0` = ``PivotEnergy``
+    :math:`({\rm MeV})`
+  * :math:`E_{\rm cut}` = ``CutoffEnergy``
+    :math:`({\rm MeV})`
+
+  .. warning::
+     The ``PivotEnergy`` parameter is not intended to be fitted.
+
+  .. note::
+     For compatibility with the Fermi/LAT ScienceTools the model type
+     ``ExponentialCutoffPowerLaw`` can be replaced by ``ExpCutoff`` and
+     the parameters ``CutoffEnergy`` by ``Cutoff`` and ``PivotEnergy``
+     by ``Scale``.
+
+
+Super exponentially cut-off power law
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+  .. code-block:: xml
+
+    <spectrum type="SuperExponentialCutoffPowerLaw">
+      <parameter name="Prefactor"    scale="1e-16" value="1.0" min="1e-07" max="1000.0" free="1"/>
+      <parameter name="Index1"       scale="-1"    value="2.0" min="0.0"   max="+5.0"   free="1"/>
+      <parameter name="CutoffEnergy" scale="1e6"   value="1.0" min="0.01"  max="1000.0" free="1"/>
+      <parameter name="Index2"       scale="1.0"   value="1.5" min="0.1"   max="5.0"    free="1"/>
+      <parameter name="PivotEnergy"  scale="1e6"   value="1.0" min="0.01"  max="1000.0" free="0"/>
+    </spectrum>
+
+  This spectral model component implements the super exponentially cut-off power
+  law function
+
+  .. math::
+     M_{\rm spectral}(E) = k_0 \left( \frac{E}{E_0} \right)^{\gamma}
+                           \exp \left(
+                           -\left( \frac{E}{E_{\rm cut}} \right)^{\alpha}
+                           \right)
+
+  where
+
+  * :math:`k_0` = ``Prefactor``
+    :math:`({\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm MeV}^{-1})`
+  * :math:`\gamma` = ``Index1``
+  * :math:`\alpha` = ``Index2``
+  * :math:`E_0` = ``PivotEnergy``
+    :math:`({\rm MeV})`
+  * :math:`E_{\rm cut}` = ``CutoffEnergy``
+    :math:`({\rm MeV})`
+
+  .. warning::
+     The ``PivotEnergy`` parameter is not intended to be fitted.
+
+  An alternative XML format is supported for compatibility with the Fermi/LAT
+  XML format:
+
+  .. code-block:: xml
+
+    <spectrum type="PLSuperExpCutoff">
+      <parameter name="Prefactor"   scale="1e-16" value="1.0" min="1e-07" max="1000.0" free="1"/>
+      <parameter name="Index1"      scale="-1"    value="2.0" min="0.0"   max="+5.0"   free="1"/>
+      <parameter name="Cutoff"      scale="1e6"   value="1.0" min="0.01"  max="1000.0" free="1"/>
+      <parameter name="Index2"      scale="1.0"   value="1.5" min="0.1"   max="5.0"    free="1"/>
+      <parameter name="Scale"       scale="1e6"   value="1.0" min="0.01"  max="1000.0" free="0"/>
+    </spectrum>
+
+
+Broken power law
+^^^^^^^^^^^^^^^^
+
+  .. code-block:: xml
+
+    <spectrum type="BrokenPowerLaw">
+      <parameter name="Prefactor"   scale="1e-16" value="5.7"  min="1e-07" max="1000.0" free="1"/>
+      <parameter name="Index1"      scale="-1"    value="2.48" min="0.0"   max="+5.0"   free="1"/>
+      <parameter name="BreakEnergy" scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="1"/>
+      <parameter name="Index2"      scale="-1"    value="2.70" min="0.01"  max="1000.0" free="1"/>
+    </spectrum>
+
+  This spectral model component implements the broken power law function
+
+  .. math::
+     M_{\rm spectral}(E) = k_0 \times \left \{
+     \begin{eqnarray}
+       \left( \frac{E}{E_b} \right)^{\gamma_1} & {\rm if\,\,} E < E_b \\
+       \left( \frac{E}{E_b} \right)^{\gamma_2} & {\rm otherwise}
+     \end{eqnarray}
+     \right .
+
+  where
+
+  * :math:`k_0` = ``Prefactor``
+    :math:`({\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm MeV}^{-1})`
+  * :math:`\gamma_1` = ``Index1``
+  * :math:`\gamma_2` = ``Index2``
+  * :math:`E_b` = ``BreakEnergy``
+    :math:`({\rm MeV})`
+
+  .. warning::
+     Note that the ``BreakEnergy`` parameter may be poorly constrained if
+     there is no clear spectral cut-off in the spectrum.
+     This model may lead to complications in the maximum likelihood fitting.
+
+  .. note::
+     For compatibility with the Fermi/LAT ScienceTools the parameters
+     ``BreakEnergy`` can be replaced by ``BreakValue``.
+
+
+Smoothly broken power law
+^^^^^^^^^^^^^^^^^^^^^^^^^
+
+  .. code-block:: xml
+
+     <spectrum type="SmoothBrokenPowerLaw">
+       <parameter name="Prefactor"       scale="1e-16" value="5.7"  min="1e-07" max="1000.0" free="1"/>
+       <parameter name="Index1"          scale="-1"    value="2.48" min="0.0"   max="+5.0"   free="1"/>
+       <parameter name="PivotEnergy"     scale="1e6"   value="1.0"  min="0.01"  max="1000.0" free="0"/>
+       <parameter name="Index2"          scale="-1"    value="2.70" min="0.01"  max="+5.0"   free="1"/>
+       <parameter name="BreakEnergy"     scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="1"/>
+       <parameter name="BreakSmoothness" scale="1.0"   value="0.2"  min="0.01"  max="10.0"   free="0"/>
+     </spectrum>
+
+  This spectral model component implements the smoothly broken power law function
+
+  .. math::
+     M_{\rm spectral}(E) = k_0 \left( \frac{E}{E_0} \right)^{\gamma_1}
+                           \left[ 1 +
+                           \left( \frac{E}{E_b} \right)^{\frac{\gamma_1 - \gamma_2}{\beta}}
+                           \right]^{-\beta}
+
+  where
+
+  * :math:`k_0` = ``Prefactor``
+    :math:`({\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm MeV}^{-1})`
+  * :math:`\gamma_1` = ``Index1``
+  * :math:`E_0` = ``PivotEnergy``
+  * :math:`\gamma_2` = ``Index2``
+  * :math:`E_b` = ``BreakEnergy``
+    :math:`({\rm MeV})`
+  * :math:`\beta` = ``BreakSmoothness``
+
+  .. warning::
+     The pivot energy should be set far away from the expected break energy
+     value.
+
+  .. warning::
+     When the two indices are close together, the :math:`\beta` parameter
+     becomes poorly constrained. Since the :math:`\beta` parameter also scales
+     the indices, this can cause very large errors in the estimates of the
+     various spectral parameters. In this case, consider fixing :math:`\beta`.
+
+  .. note::
+     For compatibility with the Fermi/LAT ScienceTools the parameters
+     ``PivotEnergy`` can be replaced by ``Scale``,
+     ``BreakEnergy`` by ``BreakValue`` and
+     ``BreakSmoothness`` by  ``Beta``.
+
+
+Log parabola
+^^^^^^^^^^^^
+
+  .. code-block:: xml
+
+    <spectrum type="LogParabola">
+      <parameter name="Prefactor"   scale="1e-17" value="5.878"   min="1e-07" max="1000.0" free="1"/>
+      <parameter name="Index"       scale="-1"    value="2.32473" min="0.0"   max="+5.0"   free="1"/>
+      <parameter name="Curvature"   scale="-1"    value="0.074"   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>
+
+  This spectral model component implements the log parabola function
+
+  .. math::
+     M_{\rm spectral}(E) = k_0 \left( \frac{E}{E_0} \right)^{\gamma+\eta \ln(E/E_0)}
+
+  where
+
+  * :math:`k_0` = ``Prefactor``
+    :math:`({\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm MeV}^{-1})`
+  * :math:`\gamma` = ``Index``
+  * :math:`\eta` = ``Curvature``
+  * :math:`E_0` = ``PivotEnergy``
+    :math:`({\rm MeV})`
+
+  .. warning::
+     The ``PivotEnergy`` parameter is not intended to be fitted.
+
+  An alternative XML format is supported for compatibility with the Fermi/LAT
+  XML format:
+
+  .. code-block:: xml
+
+     <spectrum type="LogParabola">
+       <parameter name="norm"  scale="1e-17" value="5.878"   min="1e-07" max="1000.0" free="1"/>
+       <parameter name="alpha" scale="1"     value="2.32473" min="0.0"   max="+5.0"   free="1"/>
+       <parameter name="beta"  scale="1"     value="0.074"   min="-5.0"  max="+5.0"   free="1"/>
+       <parameter name="Eb"    scale="1e6"   value="1.0"     min="0.01"  max="1000.0" free="0"/>
+     </spectrum>
+
+  where
+
+  * ``alpha`` = -``Index``
+  * ``beta`` = -``Curvature``
+
+
+Gaussian
+^^^^^^^^
+
+  .. code-block:: xml
+
+     <spectrum type="Gaussian">
+       <parameter name="Normalization" scale="1e-10" value="1.0"  min="1e-07" max="1000.0" free="1"/>
+       <parameter name="Mean"          scale="1e6"   value="5.0"  min="0.01"  max="100.0"  free="1"/>
+       <parameter name="Sigma"         scale="1e6"   value="1.0"  min="0.01"  max="100.0"  free="1"/>
+     </spectrum>
+
+  This spectral model component implements the gaussian function
+
+  .. math::
+     M_{\rm spectral}(E) = \frac{N_0}{\sqrt{2\pi}\sigma}
+                           \exp \left( \frac{-(E-\bar{E})^2}{2 \sigma^2} \right)
+
+  where
+
+  * :math:`N_0` = ``Normalization``
+    :math:`({\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1})`
+  * :math:`\bar{E}` = ``Mean``
+    :math:`({\rm MeV})`
+  * :math:`\sigma` = ``Sigma``
+    :math:`({\rm MeV})`
+
+
+File function
+^^^^^^^^^^^^^
+
+  .. code-block:: xml
+
+     <spectrum type="FileFunction" file="data/filefunction.txt">
+       <parameter scale="1.0" name="Normalization" min="0.0" max="1000.0" value="1.0" free="1"/>
+     </spectrum>
+
+  This spectral model component implements an arbitrary function
+  that is defined by intensity values at specific energies.
+  The energy and intensity values are defined using an ASCII file with
+  columns of energy and differential flux values.
+  Energies are given in units of
+  :math:`{\rm MeV}`,
+  intensities are given in units of
+  :math:`{\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm MeV}^{-1}`.
+  The only parameter is a multiplicative normalization:
+
+  .. math::
+     M_{\rm spectral}(E) = N_0 \left. \frac{dN}{dE} \right\rvert_{\rm file}
+
+  where
+
+  * :math:`N_0` = ``Normalization``
+
+  .. warning::
+     If the file name is given without a path it is expected that the file
+     resides in the same directory than the XML file.
+     If the file resides in a different directory, an absolute path name should
+     be specified.
+     Any environment variable present in the path name will be expanded.
+
+
+Node function
+^^^^^^^^^^^^^
+
+  .. code-block:: xml
+
+     <spectrum type="NodeFunction">
+       <node>
+         <parameter name="Energy"    scale="1.0"   value="1.0" min="0.1"   max="1.0e20" free="0"/>
+         <parameter name="Intensity" scale="1e-07" value="1.0" min="1e-07" max="1000.0" free="1"/>
+       </node>
+       <node>
+         <parameter name="Energy"    scale="10.0"  value="1.0" min="0.1"   max="1.0e20" free="0"/>
+         <parameter name="Intensity" scale="1e-08" value="1.0" min="1e-07" max="1000.0" free="1"/>
+       </node>
+     </spectrum>
+
+  This spectral model component implements a generalised broken 
+  power law which is defined by a set of energy and intensity values
+  (the so called nodes) that are piecewise connected by power laws.
+  Energies are given in units of
+  :math:`{\rm MeV}`,
+  intensities are given in units of
+  :math:`{\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm MeV}^{-1}`.
+
+  .. warning::
+     An arbitrary number of energy-intensity nodes can be defined in a node
+     function.
+     The nodes need to be sorted by increasing energy.
+     Although the fitting of the ``Energy`` parameters is formally possible
+     it may lead to numerical complications.
+     If ``Energy`` parameters are to be fitted make sure that the ``min``
+     and ``max`` attributes are set in a way that avoids inversion of the energy
+     ordering.
+
+
+Composite model
+^^^^^^^^^^^^^^^
+
+  .. code-block:: xml
+
+     <spectrum type="Composite">
+       <spectrum type="PowerLaw" component="SoftComponent">
+         <parameter name="Prefactor"   scale="1e-17" value="3"  min="1e-07" max="1000.0" free="1"/>
+         <parameter name="Index"       scale="-1"    value="3.5" min="0.0"   max="+5.0"   free="1"/>
+         <parameter name="PivotEnergy" scale="1e6"   value="1"  min="0.01"  max="1000.0" free="0"/>
+       </spectrum>
+       <spectrum type="PowerLaw" component="HardComponent">
+         <parameter name="Prefactor"   scale="1e-17" value="5"  min="1e-07" max="1000.0" free="1"/>
+         <parameter name="Index"       scale="-1"    value="2.0" min="0.0"   max="+5.0"   free="1"/>
+         <parameter name="PivotEnergy" scale="1e6"   value="1"  min="0.01"  max="1000.0" free="0"/>
+       </spectrum>
+     </spectrum>
+
+  This spectral model component implements a composite model that is the
+  sum of an arbitrary number of spectral models, computed using
+
+  .. math::
+     M_{\rm spectral}(E) = \sum_{i=0}^{N-1} M_{\rm spectral}^{(i)}(E)
+
+  where :math:`M_{\rm spectral}^{(i)}(E)` is any spectral model component
+  (including another composite model), and :math:`N` is the number of
+  model components that are combined.
+
+
+Multiplicative model
+^^^^^^^^^^^^^^^^^^^^
+
+  .. code-block:: xml
+
+     <spectrum type="Multiplicative">
+       <spectrum type="PowerLaw" component="PowerLawComponent">
+         <parameter name="Prefactor"   scale="1e-17" value="1.0"  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="1.0"  min="0.01"  max="1000.0" free="0"/>
+       </spectrum>
+       <spectrum type="ExponentialCutoffPowerLaw" component="CutoffComponent">
+         <parameter name="Prefactor"    scale="1.0" value="1.0" min="1e-07" max="1000.0" free="0"/>
+         <parameter name="Index"        scale="1.0" value="0.0" min="-2.0"  max="+2.0"   free="0"/>
+         <parameter name="CutoffEnergy" scale="1e6" value="1.0" min="0.01"  max="1000.0" free="1"/>
+         <parameter name="PivotEnergy"  scale="1e6" value="1.0" min="0.01"  max="1000.0" free="0"/>
+       </spectrum>
+     </spectrum>
+
+  This spectral model component implements a composite model that is the
+  product of an arbitrary number of spectral models, computed using
+
+  .. math::
+     M_{\rm spectral}(E) = \prod_{i=0}^{N-1} M_{\rm spectral}^{(i)}(E)
+
+  where :math:`M_{\rm spectral}^{(i)}(E)` is any spectral model component
+  (including another composite model), and :math:`N` is the number of
+  model components that are multiplied.
+.. _sec_models_temporal:
+
+Temporal model components
+-------------------------
+
+The following sections present the temporal model components that are available 
+in ctools.
+
+Constant
+^^^^^^^^
+
+  .. code-block:: xml
+
+     <temporal type="Constant">
+       <parameter name="Normalization" scale="1.0" value="1.0" min="0.1" max="10.0" free="0"/>
+     </temporal>
+
+  This temporal model component implements a constant source
+
+  .. math::
+     M_{\rm temporal}(t) = N_0
+
+  where
+
+  * :math:`N_0` = ``Normalization``
+
+
+Light Curve
+^^^^^^^^^^^
+
+  .. code-block:: xml
+
+     <temporal type="LightCurve" file="model_temporal_lightcurve.fits">
+       <parameter name="Normalization" scale="1" value="1.0" min="0.0" max="1000.0" free="0"/>
+     </temporal>
+
+  This temporal model component implements a light curve :math:`r(t)`
+
+  .. math::
+     M_{\rm temporal}(t) = N_0 \times r(t)
+
+  where
+
+  * :math:`N_0` = ``Normalization``
+
+  The light curve is defined by nodes in a FITS file that specify the relative
+  flux normalization as function of time (file ``model_temporal_lightcurve.fits``
+  in the example above). The structure of the light curve FITS
+  file is shown in the figure below. The light curve is defined in the first
+  extension of the FITS file and consists of a binary table with the columns
+  ``TIME`` and ``NORM``. Times in the ``TIME`` columns are given in seconds
+  and are counted with respect to a time reference that is defined in the
+  header of the binary table. Times need to be specified in ascending order.
+  The values in the ``NORM`` column specify :math:`r(t)` at times :math:`t`,
+  and should be comprised between 0 and 1.
+
+  .. _fig_model_lightcurve:
+
+  .. figure:: model_lightcurve.png
+     :align: center
+     :width: 100%
+
+     *Structure of light curve FITS file*
+
+  .. warning::
+     Fitting of light curves only makes sense for an unbinned maximum likelihood
+     analysis, since in a binned or stacked analysis the times of individual
+     events are dropped.
+
+
+Phase Curve
+^^^^^^^^^^^
+
+  .. code-block:: xml
+
+     <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>
+
+  This temporal model component implements a phase curve :math:`r(\Phi(t))`
+
+  .. math::
+     M_{\rm temporal}(t) = N_0 \times r(\Phi(t))
+
+  where the phase as function of time is computed using
+
+  .. math::
+     \Phi(t) = \Phi_0 + f(t-t_0) + \frac{1}{2}\dot{f} (t-t_0)^2 +
+                                   \frac{1}{6}\ddot{f} (t-t_0)^3
+
+  and
+
+  * :math:`N_0` = ``Normalization``
+  * :math:`t_0` = ``MJD``
+  * :math:`\Phi_0` = ``Phase``
+  * :math:`f` = ``F0``
+  * :math:`\dot{f}` = ``F1``
+  * :math:`\ddot{f}` = ``F2``
+
+  The phase curve is defined by nodes in a FITS file that specify the relative
+  flux normalization as function of phase (file ``model_temporal_phasecurve.fits``
+  in the example above). The structure of the phase curve
+  FITS file is shown in the figure below. The phase curve is defined in the
+  first extension of the FITS file and consists of a binary table with the
+  columns ``PHASE`` and ``NORM``. Phase values in the ``PHASE`` column need to
+  be comprised between 0 and 1 and need to be given in ascending order. The
+  values in the ``NORM`` column specify :math:`r(\Phi(t))` at phases
+  :math:`\Phi(t)`, and should be comprised between 0 and 1.
+
+  .. _fig_model_phasecurve:
+
+  .. figure:: model_phasecurve.png
+     :align: center
+     :width: 40%
+
+     *Structure of phase curve FITS file*
+
+  .. warning::
+     Fitting of phase curves only makes sense for an unbinned maximum likelihood
+     analysis, since in a binned or stacked analysis the times of individual
+     events are dropped.
+
+  .. warning::
+     Fitting of phase curve parameters may not properly work for pulsar
+     frequencies.
 Plotting
 ========
  
-Gammalib and ctools do not include classes and scripts to make plots.
+ctools and GammaLib do not include classes and scripts to make plots.
  
-This is intentional because creating and maintaining plotting tools that
-work on most user's machines and cover most user's needs is a big project
-in itself, and deciding to support a plotting package that is popular now
-as part of the CTA science tools could become very problematic in the
-future.
+This is intentional because creating and maintaining plotting tools that work
+on all platforms and cover all needs is a big project in itself, and deciding
+to support a plotting package that is popular now as part of the CTA science
+tools could become very problematic in the future.
  
-But we realise that most ctools users need to create plots for their
-high-level science results, so on this page we provide a few pointers 
-to external packages and tools.
+But we realise that most ctools users need to create plots for their high-level
+science results, so on this page we provide a few pointers to external packages
+and tools.
  
+Nevertheless, for your convenience several scripts for graphical display are
+included in the ctools package that rely on the
+`matplotlib <http://matplotlib.org>`_
+Python module. You can find these scripts in the
+``$CTOOLS/share/examples/python`` folder. The following scripts are available:
  
-General
--------
+  +------------------------+-----------------------------+
+  | Script                 | Usage                       |
+  +========================+=============================+
+  | ``show_butterfly.py``  | Display a butterfly diagram |
+  +------------------------+-----------------------------+
+  | ``show_lightcurve.py`` | Display a light curve       |
+  +------------------------+-----------------------------+
+  | ``show_obs.py``        | Display observation summary |
+  +------------------------+-----------------------------+
+  | ``show_pointings.py``  | Display pointing directions |
+  +------------------------+-----------------------------+
+  | ``show_phases.py``     | Display event phases        |
+  +------------------------+-----------------------------+
+  | ``show_spectrum.py``   | Display a spectrum          |
+  +------------------------+-----------------------------+
  
-The most common way to create publication-quality astronomy plots at the
-moment is to use the `matplotlib`_ Python plotting package. As illustrated
-by the `matplotlib gallery`_, virtually any plot can be constructed,
-but matplotlib doesn't have built-in easy-to-use classes or functions to
-create the common plots for gamma-ray astronomy, which are sky maps and
-spectra, so in the sections below we give some specific pointers for those.
+Below some usage examples and the expected output.
  
+show_butterfly.py
+^^^^^^^^^^^^^^^^^
  
-Sky maps
---------
+.. code-block:: bash
  
-For creating publication quality sky map plots via a Python script a good
-choice is to use `APLPy`_, `WCSAxes`_ or `kapteyn`_, which are all based on 
-matplotlib.
+   $ $CTOOLS/share/examples/python/show_butterfly.py butterfly_src001.txt
  
-You can also use interactive astronomy image viewers like `DS9`_, `Aladin`_
-or `Ginga`_ which allow you to save the image to a PNG or PDF file after
-customising it interactively, but also contain command line options to create
-the plot, i.e. you can write a bash or Python script to make the plot
-reproducible (DS9 and Ginga also feature direct communication with Python,
-but that is a bit harder to use and usually used for interactive image
-analysis, not plotting scripts).
+.. figure:: plotting_butterfly.png
+   :width: 600px
+   :align: center
  
+   *Butterfly diagram displayed with show_butterfly.py*
  
-Spectra and light curves
-------------------------
  
-To plot spectra or spectral energy distributions (including flux model curves,
-flux model error bands, flux points and residuals) and light curves, you can
-use matplotlib directly, or use the classes and scripts in `Gammapy`_.
+show_lightcurve.py
+^^^^^^^^^^^^^^^^^^
  
-Another option is to use virtual observatory interactive spectrum tools
-like `IRIS`_, `CASSIS`_ or `VOSpec`_ which allow you to save the plot to a
-PNG or PDF file after creating it interactively.
+.. code-block:: bash
  
-.. _DS9: http://ds9.si.edu/
-.. _Aladin: http://aladin.u-strasbg.fr/
-.. _Ginga: http://ejeschke.github.io/ginga/
-.. _matplotlib: http://matplotlib.org/
-.. _matplotlib gallery: http://matplotlib.org/gallery.html
-.. _WCSAxes: http://wcsaxes.readthedocs.org/
-.. _APLPy: http://aplpy.github.io/
-.. _kapteyn: https://www.astro.rug.nl/software/kapteyn/
-.. _Gammapy: https://gammapy.readthedocs.org/
-.. _IRIS: http://cxc.cfa.harvard.edu/iris/
-.. _CASSIS: http://cassis.irap.omp.eu/
-.. _VOSpec: http://www.sciops.esa.int/index.php?project=SAT&page=vospec
+   $ $CTOOLS/share/examples/python/show_lightcurve.py lightcurve.fits
  
-.. raw:: html
+.. figure:: plotting_lightcurve.png
+   :width: 600px
+   :align: center
  
-   <iframe src="http://nbviewer.jupyter.org/github/davidchall/random-direction/blob/master/Generate_Beam_Directions.ipynb" width="100%" height="800"></iframe>
+   *Light curve displayed with show_lightcurve.py*
+
+
+show_obs.py
+^^^^^^^^^^^
+
+.. code-block:: bash
+
+   $ $CTOOLS/share/examples/python/show_obs.py obs_selected.xml
+
+.. figure:: plotting_obs.png
+   :width: 600px
+   :align: center
+
+   *Observation summary displayed with show_obs.py*
+
+
+show_pointings.py
+^^^^^^^^^^^^^^^^^
+
+.. code-block:: bash
+
+   $ $CTOOLS/share/examples/python/show_pointings.py obs_selected.xml
+
+.. figure:: plotting_pointings.png
+   :width: 600px
+   :align: center
+
+   *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:: plotting_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
+
+.. figure:: plotting_spectrum.png
+   :width: 600px
+   :align: center
+
+   *Observation summary displayed with show_spectrum.py*
@@ -48,13 +48,38 @@ following:
     PSF(p' | p, E, t) \times
     E_{\rm disp}(E' | p, E, t)
  
-ctools are shipped with response functions for the northern and southern
-arrays, and variants are available that have been optimised for exposure
-times of 0.5 hours, 5 hours and 50 hours.
+ctools are shipped with so called ``prod2`` response functions for the northern
+and southern arrays, and variants are available that have been optimised for
+exposure times of 0.5 hours, 5 hours and 50 hours.
 In total, the following six instrument response functions are available:
 ``North_0.5h``, ``North_5h``, ``North_50h``, ``South_0.5h``,
 ``South_5h``, and ``South_50h``.
  
+.. note::
+   In cases for which the CTA instrument response functions are used in a
+   research project, we ask to add the following acknowledgement in any
+   resulting publication:
+
+   *This research has made use of the CTA instrument response functions
+   provided by the CTA Consortium and Observatory.*
+
+   The latest CTA reponse function (prod3b-v1) can be downloaded from
+   `this link <http://www.cta-observatory.org/wp-content/uploads/2017/12/CTA-Performance-prod3b-v1-FITS1.tar.gz>`_. To install the response function, type
+
+   .. code-block:: bash
+
+      $ tar xvfz CTA-Performance-prod3b-v1-FITS1.tar.gz -C $CTOOLS/share
+
+   or
+
+   .. code-block:: bash
+
+      $ sudo tar xvfz CTA-Performance-prod3b-v1-FITS1.tar.gz -C $CTOOLS/share
+
+   if you need root priviledges to install in your ``$CTOOLS`` directory.
+   To use the prod3b-v1 response functions, specify ``prod3b-v1`` if a tool
+   or script queries for the calibration database.
+
 Each response is stored in a single FITS file, and each component of
 the response factorisation is stored in a binary table of that FITS
 file.
@@ -103,18 +128,18 @@ calibration database and the ``South_0.5h`` response function:
  
 .. code-block:: bash
  
-  $ ctobssim
-  RA of pointing (degrees) (0-360) [83.63] 
-  Dec of pointing (degrees) (-90-90) [22.51] 
-  Radius of FOV (degrees) (0-180) [5.0] 
-  Start time (MET in s) [0.0] 
-  End time (MET in s) [1800.0] 
-  Lower energy limit (TeV) [0.1] 
-  Upper energy limit (TeV) [100.0] 
-  Calibration database [prod2] 
-  Instrument response function [South_0.5h] 
-  Input model XML file [$CTOOLS/share/models/crab.xml] 
-  Output event data file or observation definition XML file [events.fits]
+   $ ctobssim
+   RA of pointing (degrees) (0-360) [83.63]
+   Dec of pointing (degrees) (-90-90) [22.51]
+   Radius of FOV (degrees) (0-180) [5.0]
+   Start time (MET in s) [0.0]
+   End time (MET in s) [1800.0]
+   Lower energy limit (TeV) [0.1]
+   Upper energy limit (TeV) [100.0]
+   Calibration database [prod2]
+   Instrument response function [South_0.5h]
+   Input model XML file [$CTOOLS/share/models/crab.xml]
+   Output event data file or observation definition XML file [events.fits]
  
 Running the other tools is equivalent.
  
@@ -179,11 +204,11 @@ and response name from within Python:
  
 .. code-block:: python
  
-  import gammalib
-  obs   = gammalib.GCTAObservation()
-  caldb = gammalib.GCaldb("cta", "prod2")
-  irf   = "South_0.5h"
-  obs.response(irf, caldb)
+   >>> import gammalib
+   >>> obs   = gammalib.GCTAObservation()
+   >>> caldb = gammalib.GCaldb("cta", "prod2")
+   >>> irf   = "South_0.5h"
+   >>> obs.response(irf, caldb)
  
 The calibration database is set by creating a ``GCaldb`` object. The
 constructor takes as argument the mission (always ``cta``) and the