From c8fd2ec62f2702e2f369ccf8ecaf38becc1a0594 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Wed, 24 Jan 2018 14:38:12 +0100 Subject: [PATCH] Update ctools User manual --- doc/source/users/user_manual/models.rst | 1072 +++++++++------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- doc/source/users/user_manual/models_concept.rst | 34 ++++++++++++++++++++++++++++++++++ doc/source/users/user_manual/models_implementation.rst | 95 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/users/user_manual/models_spatial_bgd.rst | 181 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/users/user_manual/models_spatial_src.rst | 354 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/users/user_manual/models_spectral.rst | 503 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/users/user_manual/models_temporal.rst | 130 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/users/user_manual/plotting.rst | 152 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------------------------------- doc/source/users/user_manual/plotting_butterfly.png | Bin 0 -> 33989 bytes doc/source/users/user_manual/plotting_lightcurve.png | Bin 0 -> 24732 bytes doc/source/users/user_manual/plotting_obs.png | Bin 0 -> 73814 bytes doc/source/users/user_manual/plotting_phases.png | Bin 0 -> 18403 bytes doc/source/users/user_manual/plotting_pointings.png | Bin 0 -> 27281 bytes doc/source/users/user_manual/plotting_spectrum.png | Bin 0 -> 23006 bytes doc/source/users/user_manual/response.rst | 65 +++++++++++++++++++++++++++++++++++++++++++++-------------------- 15 files changed, 1452 insertions(+), 1134 deletions(-) create mode 100644 doc/source/users/user_manual/models_concept.rst create mode 100644 doc/source/users/user_manual/models_implementation.rst create mode 100644 doc/source/users/user_manual/models_spatial_bgd.rst create mode 100644 doc/source/users/user_manual/models_spatial_src.rst create mode 100644 doc/source/users/user_manual/models_spectral.rst create mode 100644 doc/source/users/user_manual/models_temporal.rst create mode 100644 doc/source/users/user_manual/plotting_butterfly.png create mode 100644 doc/source/users/user_manual/plotting_lightcurve.png create mode 100644 doc/source/users/user_manual/plotting_obs.png create mode 100644 doc/source/users/user_manual/plotting_phases.png create mode 100644 doc/source/users/user_manual/plotting_pointings.png create mode 100644 doc/source/users/user_manual/plotting_spectrum.png diff --git a/doc/source/users/user_manual/models.rst b/doc/source/users/user_manual/models.rst index f9ebdd2..7d1eda1 100644 --- a/doc/source/users/user_manual/models.rst +++ b/doc/source/users/user_manual/models.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) ` -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 ` 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 ` 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 ````), -a spatial (tags ```` and ````), and -a temporal component (tag ````). - -.. 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 - - - - - - - - - - - - - - - - - - - - - - - - - - - -Model parameters are specified by a ```` 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 ```` tag. - -.. note:: - For compatibility reasons with the Fermi/LAT XML format the ```` - 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 - - - - - - - - ... - - - - .. note:: - - For compatibility with the Fermi/LAT ScienceTools the model type - ``PointSource`` can be replaced by ``SkyDirFunction``. - - -Radial source -^^^^^^^^^^^^^ - - .. code-block:: xml - - - - - - - - - ... - - - - .. code-block:: xml - - - - - - - - - ... - - - - .. code-block:: xml - - - - - - - - - - ... - - - - -Elliptical source -^^^^^^^^^^^^^^^^^ - - .. code-block:: xml - - - - - - - - - - - ... - - - - .. code-block:: xml - - - - - - - - - - - ... - - - - -Diffuse source -^^^^^^^^^^^^^^ - - .. code-block:: xml - - - - - - - ... - - - - .. note:: - - For compatibility with the Fermi/LAT ScienceTools the model type - ``DiffuseIsotropic`` can be replaced by ``ConstantValue``. - - .. code-block:: xml - - - - - - - ... - - - - .. 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 - - - - - - - ... - - - - .. 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 - - - - - - - - - - - - - - - ... - - - - 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 - - - - - - - ... - - - - .. code-block:: xml - - - - - - - - - ... - - - - .. code-block:: xml - - - - - - - - - - - - - - - ... - - - - -CTA IRF background -^^^^^^^^^^^^^^^^^^ - - .. code-block:: xml - - - - ... - - - - -CTA cube background -^^^^^^^^^^^^^^^^^^^ - - .. code-block:: xml - - - - ... - - - - - -.. _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 - - - - - - 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 - - - - - - - - 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 - - - - - - - - - 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 - - - - - - - - - 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 - - - - - - - - - - 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 - - - - - - - - - - -Broken power law -^^^^^^^^^^^^^^^^ - - .. code-block:: xml - - - - - - - - - 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 - - - - - - - - - - - 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 - - - - - - - - - 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 - - - - - - - - - where - - * ``alpha`` = -``Index`` - * ``beta`` = -``Curvature`` - - -Gaussian -^^^^^^^^ - - .. code-block:: xml - - - - - - - - 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 - - - - - - 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 - - - - - - - - - - - - - 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 - - - - - - - - - - - - - - - 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 - - - - - - - - - - - - - - - - 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 - - - - - - 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 - - - - - - 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 - - - - - - - - - - - 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 diff --git a/doc/source/users/user_manual/models_concept.rst b/doc/source/users/user_manual/models_concept.rst new file mode 100644 index 0000000..c3d148d --- /dev/null +++ b/doc/source/users/user_manual/models_concept.rst @@ -0,0 +1,34 @@ +.. _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) ` +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 ` 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 ` convolution. diff --git a/doc/source/users/user_manual/models_implementation.rst b/doc/source/users/user_manual/models_implementation.rst new file mode 100644 index 0000000..fd2512b --- /dev/null +++ b/doc/source/users/user_manual/models_implementation.rst @@ -0,0 +1,95 @@ +.. _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 ````), +a spatial (tags ```` and ````), and +a temporal component (tag ````). + +.. 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) ` +and that is multiplied by a spectral power law function. + +.. code-block:: xml + + + + + + + + + + + + + + + + + + + + + + + + + + +Model parameters are specified by a ```` 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 ```` tag. + +.. note:: + For compatibility reasons with the Fermi/LAT XML format the ```` + 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. diff --git a/doc/source/users/user_manual/models_spatial_bgd.rst b/doc/source/users/user_manual/models_spatial_bgd.rst new file mode 100644 index 0000000..9f61aae --- /dev/null +++ b/doc/source/users/user_manual/models_spatial_bgd.rst @@ -0,0 +1,181 @@ +.. _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 ```` tag as the spatial component. + + The first radial model is of type ``Gaussian`` + + .. code-block:: xml + + + + + + + + + + + + + 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 + + + + + + + + + + + + + + + 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 + + + + + + + + + + + + + + + + + + + + + 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) ` 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 + + + + + + + + + + If the observation is an On/Off observation, do not forget to switch the instrument + to ``CTAOnOff``: + + .. code-block:: xml + + + + + + + + + + +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 + + + + + + + + + + +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 + + + + + + + + diff --git a/doc/source/users/user_manual/models_spatial_src.rst b/doc/source/users/user_manual/models_spatial_src.rst new file mode 100644 index 0000000..8fb2e8e --- /dev/null +++ b/doc/source/users/user_manual/models_spatial_src.rst @@ -0,0 +1,354 @@ +.. _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 + + + + + + + + ... + + + + 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 + + + + + + + + ... + + + + 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 + + + + + + + + + ... + + + + 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 + + + + + + + + + ... + + + + 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 + + + + + + + + + + ... + + + + 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 + + + + + + + + + + + ... + + + + 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 + + + + + + + + + + + ... + + + + 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 + + + + + + + ... + + + + 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 + + + + + + + ... + + + + 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 + + + + + + + ... + + + + The ``DiffuseMapCube`` model describes an arbitrary energy-dependent intensity + distribution in form of a map cube + + .. code-block:: xml + + + + + + + ... + + + + 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 + + + + + + + + + + + + + + + ... + + + + 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. diff --git a/doc/source/users/user_manual/models_spectral.rst b/doc/source/users/user_manual/models_spectral.rst new file mode 100644 index 0000000..f4a9604 --- /dev/null +++ b/doc/source/users/user_manual/models_spectral.rst @@ -0,0 +1,503 @@ +.. _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 + + + + + + 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 + + + + + + + + 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 + + + + + + + + + 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 + + + + + + + + + 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 + + + + + + + + + + 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 + + + + + + + + + + +Broken power law +^^^^^^^^^^^^^^^^ + + .. code-block:: xml + + + + + + + + + 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 + + + + + + + + + + + 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 + + + + + + + + + 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 + + + + + + + + + where + + * ``alpha`` = -``Index`` + * ``beta`` = -``Curvature`` + + +Gaussian +^^^^^^^^ + + .. code-block:: xml + + + + + + + + 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 + + + + + + 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 + + + + + + + + + + + + + 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 + + + + + + + + + + + + + + + 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 + + + + + + + + + + + + + + + + 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. diff --git a/doc/source/users/user_manual/models_temporal.rst b/doc/source/users/user_manual/models_temporal.rst new file mode 100644 index 0000000..77ab3d0 --- /dev/null +++ b/doc/source/users/user_manual/models_temporal.rst @@ -0,0 +1,130 @@ +.. _sec_models_temporal: + +Temporal model components +------------------------- + +The following sections present the temporal model components that are available +in ctools. + +Constant +^^^^^^^^ + + .. code-block:: xml + + + + + + 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 + + + + + + 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 + + + + + + + + + + + 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. diff --git a/doc/source/users/user_manual/plotting.rst b/doc/source/users/user_manual/plotting.rst index 13b629a..b905b75 100644 --- a/doc/source/users/user_manual/plotting.rst +++ b/doc/source/users/user_manual/plotting.rst @@ -1,70 +1,120 @@ 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 `_ +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 - + *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* diff --git a/doc/source/users/user_manual/plotting_butterfly.png b/doc/source/users/user_manual/plotting_butterfly.png new file mode 100644 index 0000000..ed36522 Binary files /dev/null and b/doc/source/users/user_manual/plotting_butterfly.png differ diff --git a/doc/source/users/user_manual/plotting_lightcurve.png b/doc/source/users/user_manual/plotting_lightcurve.png new file mode 100644 index 0000000..4766a35 Binary files /dev/null and b/doc/source/users/user_manual/plotting_lightcurve.png differ diff --git a/doc/source/users/user_manual/plotting_obs.png b/doc/source/users/user_manual/plotting_obs.png new file mode 100644 index 0000000..bec5639 Binary files /dev/null and b/doc/source/users/user_manual/plotting_obs.png differ diff --git a/doc/source/users/user_manual/plotting_phases.png b/doc/source/users/user_manual/plotting_phases.png new file mode 100644 index 0000000..d893ad0 Binary files /dev/null and b/doc/source/users/user_manual/plotting_phases.png differ diff --git a/doc/source/users/user_manual/plotting_pointings.png b/doc/source/users/user_manual/plotting_pointings.png new file mode 100644 index 0000000..938c51c Binary files /dev/null and b/doc/source/users/user_manual/plotting_pointings.png differ diff --git a/doc/source/users/user_manual/plotting_spectrum.png b/doc/source/users/user_manual/plotting_spectrum.png new file mode 100644 index 0000000..6621fd2 Binary files /dev/null and b/doc/source/users/user_manual/plotting_spectrum.png differ diff --git a/doc/source/users/user_manual/response.rst b/doc/source/users/user_manual/response.rst index 599cd09..089a510 100644 --- a/doc/source/users/user_manual/response.rst +++ b/doc/source/users/user_manual/response.rst @@ -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 `_. 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 -- libgit2 0.21.4