From 58fd52f72d65dc7e3548964d23729ad45c089237 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Thu, 20 Jul 2017 21:25:23 +0200 Subject: [PATCH] Update the User documentation --- doc/source/users/reference_manual/csroot2caldb.rst | 108 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/users/reference_manual/index.rst | 1 + doc/source/users/user_manual/getting_started/index.rst | 12 ------------ doc/source/users/user_manual/getting_started/irf-aeff.png | Bin 3319 -> 0 bytes doc/source/users/user_manual/getting_started/irf-file.png | Bin 7266 -> 0 bytes doc/source/users/user_manual/getting_started/irfs.jpg | Bin 28217 -> 0 bytes doc/source/users/user_manual/getting_started/model_lightcurve.png | Bin 66994 -> 0 bytes doc/source/users/user_manual/getting_started/model_phasecurve.png | Bin 27527 -> 0 bytes doc/source/users/user_manual/getting_started/models.rst | 1068 ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ doc/source/users/user_manual/getting_started/python.rst | 448 ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- doc/source/users/user_manual/getting_started/response.rst | 192 ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ doc/source/users/user_manual/index.rst | 8 +++++--- doc/source/users/user_manual/introduction.rst | 3 +-- doc/source/users/user_manual/irf-aeff.png | Bin 0 -> 3319 bytes doc/source/users/user_manual/irf-file.png | Bin 0 -> 7266 bytes doc/source/users/user_manual/irfs.jpg | Bin 0 -> 28217 bytes doc/source/users/user_manual/model_lightcurve.png | Bin 0 -> 66994 bytes doc/source/users/user_manual/model_phasecurve.png | Bin 0 -> 27527 bytes doc/source/users/user_manual/models.rst | 1068 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/users/user_manual/python.rst | 454 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/users/user_manual/response.rst | 192 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 21 files changed, 1829 insertions(+), 1725 deletions(-) create mode 100644 doc/source/users/reference_manual/csroot2caldb.rst delete mode 100644 doc/source/users/user_manual/getting_started/index.rst delete mode 100644 doc/source/users/user_manual/getting_started/irf-aeff.png delete mode 100644 doc/source/users/user_manual/getting_started/irf-file.png delete mode 100644 doc/source/users/user_manual/getting_started/irfs.jpg delete mode 100644 doc/source/users/user_manual/getting_started/model_lightcurve.png delete mode 100644 doc/source/users/user_manual/getting_started/model_phasecurve.png delete mode 100644 doc/source/users/user_manual/getting_started/models.rst delete mode 100644 doc/source/users/user_manual/getting_started/python.rst delete mode 100644 doc/source/users/user_manual/getting_started/response.rst create mode 100644 doc/source/users/user_manual/irf-aeff.png create mode 100644 doc/source/users/user_manual/irf-file.png create mode 100644 doc/source/users/user_manual/irfs.jpg create mode 100644 doc/source/users/user_manual/model_lightcurve.png create mode 100644 doc/source/users/user_manual/model_phasecurve.png create mode 100644 doc/source/users/user_manual/models.rst create mode 100644 doc/source/users/user_manual/python.rst create mode 100644 doc/source/users/user_manual/response.rst diff --git a/doc/source/users/reference_manual/csroot2caldb.rst b/doc/source/users/reference_manual/csroot2caldb.rst new file mode 100644 index 0000000..ccf4e94 --- /dev/null +++ b/doc/source/users/reference_manual/csroot2caldb.rst @@ -0,0 +1,108 @@ +.. _csroot2caldb: + +csroot2caldb +============ + +Creates a caldb entry from a ROOT performance file. + + +Synopsis +-------- + +This script creates an entry in the local calibration database from a ROOT +performance file. + + +General parameters +------------------ + +``infile [file]`` + Input ROOT performance file file. + +``(outdir = ) [string]`` + Output CALDB directory. + +``inst [string]`` + Instrument name. + +``id [string]`` + Instrument Response Function identifier. + +``version [string]`` + Instrument Response Function version. + +``analysis [string]`` + Instrument Response Function analysis. + +``zenith [real]`` + Zenith angle (in degrees). + +``azimuth [real]`` + Azimuth angle (in degrees). + +``(emin = UNDEF) [real]`` + Minimum energy threshold (TeV). + +``(emax = UNDEF) [real]`` + Maximum energy threshold (TeV). + +``(psftype = Gauss) [string]`` + Point spread function type. + +``(split = no) [boolean]`` + Split IRF components into different files. + +``(norm1d = no) [boolean]`` + Normalize on 1D histograms. + +``(rebin = no) [boolean]`` + Rebin. + +``(eascale = 1.0) [real]`` + Scaling factor for effective areas. + +``(bgdscale = 1.0) [real]`` + Scaling factor for background rates. + +``(bgdoversample = 1) [integer]`` + Spatial oversampling factor for background template. + +``(bgdethres = 1000.0) [real]`` + Energy above which to replace background by power law (TeV). + +``(bgdinfill = no) [boolean]`` + Infill of background template?. + + +Standard parameters +------------------- + +``(chatter = 2) [integer]`` + Verbosity of the executable: + ``chatter = 0``: no information will be logged + + ``chatter = 1``: only errors will be logged + + ``chatter = 2``: errors and actions will be logged + + ``chatter = 3``: report about the task execution + + ``chatter = 4``: detailed report about the task execution + +``(clobber = yes) [boolean]`` + Specifies whether an existing output file should be overwritten. + +``(debug = no) [boolean]`` + Enables debug mode. In debug mode the executable will dump any log file output to the console. + +``(mode = ql) [string]`` + Mode of automatic parameters (default is "ql", i.e. "query and learn"). + +``(logfile = csroot2caldb.log) [filename]`` + Log filename. + + +Related tools or scripts +------------------------ + +:doc:`cscaldb` diff --git a/doc/source/users/reference_manual/index.rst b/doc/source/users/reference_manual/index.rst index f81202a..9e6c17b 100644 --- a/doc/source/users/reference_manual/index.rst +++ b/doc/source/users/reference_manual/index.rst @@ -53,6 +53,7 @@ cscripts csphasecrv --- Computes phase curve cspull --- Generates pull distribution csresmap --- Generates residual map + csroot2caldb --- Creates a caldb entry from a ROOT file cssens --- Computes CTA sensitivity csspec --- Computes spectral points cssrcdetect --- Detects sources in sky map diff --git a/doc/source/users/user_manual/getting_started/index.rst b/doc/source/users/user_manual/getting_started/index.rst deleted file mode 100644 index 2ce9cfd..0000000 --- a/doc/source/users/user_manual/getting_started/index.rst +++ /dev/null @@ -1,12 +0,0 @@ -.. _sec_using: - -Using ctools -============ - -.. toctree:: - :maxdepth: 1 - - python - models - response - diff --git a/doc/source/users/user_manual/getting_started/irf-aeff.png b/doc/source/users/user_manual/getting_started/irf-aeff.png deleted file mode 100644 index eabbe69..0000000 Binary files a/doc/source/users/user_manual/getting_started/irf-aeff.png and /dev/null differ diff --git a/doc/source/users/user_manual/getting_started/irf-file.png b/doc/source/users/user_manual/getting_started/irf-file.png deleted file mode 100644 index 31f1599..0000000 Binary files a/doc/source/users/user_manual/getting_started/irf-file.png and /dev/null differ diff --git a/doc/source/users/user_manual/getting_started/irfs.jpg b/doc/source/users/user_manual/getting_started/irfs.jpg deleted file mode 100644 index c60bb63..0000000 Binary files a/doc/source/users/user_manual/getting_started/irfs.jpg and /dev/null differ diff --git a/doc/source/users/user_manual/getting_started/model_lightcurve.png b/doc/source/users/user_manual/getting_started/model_lightcurve.png deleted file mode 100644 index 27119e0..0000000 Binary files a/doc/source/users/user_manual/getting_started/model_lightcurve.png and /dev/null differ diff --git a/doc/source/users/user_manual/getting_started/model_phasecurve.png b/doc/source/users/user_manual/getting_started/model_phasecurve.png deleted file mode 100644 index 9ddaa6f..0000000 Binary files a/doc/source/users/user_manual/getting_started/model_phasecurve.png and /dev/null differ diff --git a/doc/source/users/user_manual/getting_started/models.rst b/doc/source/users/user_manual/getting_started/models.rst deleted file mode 100644 index f9ebdd2..0000000 --- a/doc/source/users/user_manual/getting_started/models.rst +++ /dev/null @@ -1,1068 +0,0 @@ -.. _sec_models: - -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. diff --git a/doc/source/users/user_manual/getting_started/python.rst b/doc/source/users/user_manual/getting_started/python.rst deleted file mode 100644 index 8263211..0000000 --- a/doc/source/users/user_manual/getting_started/python.rst +++ /dev/null @@ -1,448 +0,0 @@ -.. _sec_python: - -Using ctools from Python ------------------------- - -Getting started -~~~~~~~~~~~~~~~ - -ctools provides a Python module that allows using all tools and scripts as -Python classes. -To use ctools from Python all you have to do is to import the ctools -module into Python. -You should also import the GammaLib module, as ctools without GammaLib is -generally not very useful. - -.. code-block:: python - - >>> import gammalib - >>> import ctools - -To illustrate how to use ctools from Python, below is a working example -of an event list simulation using the :ref:`ctobssim` class. -An instance ``sim`` of the :ref:`ctobssim` class is generated and user -parameters are set for this instance using the ``[ ]`` operator. -See the :ref:`reference` for a list of parameters and their types). -The ``execute()`` method executes the :ref:`ctobssim` class in the same -way as if it were executed from the command line. - -.. code-block:: python - - >>> import ctools - >>> sim = ctools.ctobssim() - >>> sim["inmodel"] = "${CTOOLS}/share/models/crab.xml" - >>> sim["outevents"] = "events.fits" - >>> sim["caldb"] = "prod2" - >>> sim["irf"] = "South_0.5h" - >>> sim["ra"] = 83.63 - >>> sim["dec"] = 22.01 - >>> sim["rad"] = 5.0 - >>> sim["tmin"] = 0.0 - >>> sim["tmax"] = 1800.0 - >>> sim["emin"] = 0.1 - >>> sim["emax"] = 100.0 - >>> sim.execute() - -Alternatively, you may "run" the :ref:`ctobssim` tool using - -.. code-block:: python - - >>> sim.run() - -The main difference to the ``execute()`` method is that the ``run()`` -will not write the simulated event file to disk. -Why is this useful? -Well, after having typed ``sim.run()`` the :ref:`ctobssim` class still -exists as an object in memory, including all the simulated events. -The :ref:`ctobssim` class has an ``obs()`` method that returns an -observation container that holds the simulated CTA observation with its -associated events. -To visualise this container, type: - -.. code-block:: python - - >>> print(sim.obs()) - === GObservations === - Number of observations ....: 1 - Number of models ..........: 2 - Number of observed events .: 23099 - Number of predicted events : 0 - -There is one CTA observation in the container and to visualise the events -in that observation you may type: - -.. code-block:: python - - >>> print(sim.obs()[0].events()) - === GCTAEventList === - Number of events ..........: 23099 (disposed in "events.fits") - Time interval .............: 51544.5 - 51544.5 days - === GEbounds === - Number of intervals .......: 1 - Energy range ..............: 100 GeV - 100 TeV - === GCTARoi === - ROI centre ................: RA=83.63, DEC=22.01 [0,0] - ROI radius ................: 5 deg - -The ``obs()[0]`` operator returns the first observation in the observation -container, the ``events()`` operator returns the event list in that -observation. -To see what kind of object you actually got, use: - -.. code-block:: python - - >>> type(sim.obs()[0].events()) - - -The CTA event list is implement by the ``GCTAEventList`` class in the -``cta`` module of GammaLib. -To visualise the individual events you can iterate over the events using a -for loop. -This will show the simulated celestial coordinates (RA, DEC), the -coordinate in the camera system [DETX, DETY], energies and -terrestrial times (TT) of all events. - -.. code-block:: python - - >>> events = sim.obs()[0].events() - >>> for event in events: - ... print(event) - ... - Dir=RA=83.6308, DEC=21.8881 [-0.00212759,1.33661e-05] Energy=106.465 GeV Time=-3.15576e+08 s (TT) - Dir=RA=83.7518, DEC=21.8064 [-0.0035525,0.00197398] Energy=117.706 GeV Time=-3.15576e+08 s (TT) - Dir=RA=83.5545, DEC=22.0933 [0.00145377,-0.00122121] Energy=138.624 GeV Time=-3.15576e+08 s (TT) - ... - - -Performing a maximum likelihood analysis in Python -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -We can now benefit from the fact that we have some simulated events in -memory to fit a model to these events using the :ref:`ctlike` class. -We will do this in unbinned mode. -Here is what you have to do: - -.. code-block:: python - - >>> like = ctools.ctlike(sim.obs()) - >>> like.run() - -This is pretty compact. -Where are the user parameters? -:ref:`ctlike` doesn't in fact need any as all the relevant information is -already contained in the observation container produced by the -:ref:`ctobssim` class. -And you may have recognised that we constructed the :ref:`ctlike` -instance by using the :ref:`ctobssim` observation container as -constructor argument. - -To check how the fit went you can inspect the optimiser used by -:ref:`ctlike`: - -.. code-block:: python - - >>> print(like.opt()) - === GOptimizerLM === - Optimized function value ..: 154553.422 - Absolute precision ........: 0.005 - Acceptable value decrease .: 2 - Optimization status .......: converged - Number of parameters ......: 10 - Number of free parameters .: 4 - Number of iterations ......: 2 - Lambda ....................: 1e-05 - -Apparently, the fit converged after 2 iterations. -Out of 10 parameters in the model 4 have been fitted (the others were kept -fixed). -To inspect the fit results you can print the model container that is a -member of the observation container: - -.. code-block:: python - - >>> print(like.obs().models()) - === GModels === - Number of models ..........: 2 - Number of parameters ......: 10 - === GModelSky === - Name ......................: Crab - Instruments ...............: all - Instrument scale factors ..: unity - Observation identifiers ...: all - Model type ................: PointSource - Model components ..........: "PointSource" * "PowerLaw" * "Constant" - Number of parameters ......: 6 - Number of spatial par's ...: 2 - RA .......................: 83.6331 [-360,360] deg (fixed,scale=1) - DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1) - Number of spectral par's ..: 3 - Prefactor ................: 5.82698e-16 +/- 1.02186e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) - Index ....................: -2.47534 +/- 0.0154764 [-0,-5] (free,scale=-1,gradient) - PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) - Number of temporal par's ..: 1 - Normalization ............: 1 (relative value) (fixed,scale=1,gradient) - === GCTAModelIrfBackground === - Name ......................: CTABackgroundModel - Instruments ...............: CTA - Instrument scale factors ..: unity - Observation identifiers ...: all - Model type ................: "PowerLaw" * "Constant" - Number of parameters ......: 4 - Number of spectral par's ..: 3 - Prefactor ................: 1.01266 +/- 0.0119676 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient) - Index ....................: 0.00474762 +/- 0.00731725 [-5,5] (free,scale=1,gradient) - PivotEnergy ..............: 1e+06 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) - Number of temporal par's ..: 1 - Normalization ............: 1 (relative value) (fixed,scale=1,gradient) - -Suppose you want to repeat the fit by optimising also the position of the -point source. -This is easy from Python: - -.. code-block:: python - - >>> like.obs().models()["Crab"]["RA"].free() - >>> like.obs().models()["Crab"]["DEC"].free() - >>> like.run() - >>> print(like.obs().models()) - ... - RA .......................: 83.6327 +/- 0.000916983 [-360,360] deg (free,scale=1) - DEC ......................: 22.0141 +/- 0.00086378 [-90,90] deg (free,scale=1) - -The ``like.obs().models()`` method provides the model container, using the -``["Crab"]`` operator we access the Crab model in that container and using -the ``["RA"]`` and ``["DEC"]`` methods we access the relevant model -parameters. -The ``free()`` method frees a parameter, the opposite would be a call to -the ``fix()`` method. - -.. warning:: - - Passing observation containers between ctools classes is a very - convenient and powerful way of building in-memory analysis pipelines. - However, this implies that you need some computing ressources when - dealing with large observation containers (for example if you want to - analyse a few 100 hours of data at once). - -.. warning:: - - You have to be aware about the scope of the objects you're - dealing with. - In the above example, the ``sim.obs()`` container is allocated by the - :ref:`ctobssim` class, hence it disappears (a.k.a. goes out of scope) - once the :ref:`ctobssim` class is deleted, as illustrated by the - following example: - - .. code-block:: python - - >>> obs = sim.obs() - >>> del sim - >>> print(obs) - Segmentation fault - - Note that - - .. code-block:: python - - >>> obs = sim.obs() - >>> del sim - >>> print(like.obs()) - === GObservations === - Number of observations ....: 1 - Number of models ..........: 2 - Number of observed events .: 23099 - Number of predicted events : 23099 - - is okay as the ``ctools.ctlike(sim.obs())`` constructor will create - a copy of the observation container that lives within the :ref:`ctlike` - instance. - To preserve an observation container after a ctools object goes out - of scope you have to create a local copy of the container using the - ``copy()`` method: - - .. code-block:: python - - >>> obs = sim.obs().copy() - >>> del sim - >>> print(obs) - === GObservations === - Number of observations ....: 1 - Number of models ..........: 2 - Number of observed events .: 23099 - Number of predicted events : 0 - - -Using obsutils -~~~~~~~~~~~~~~ - -ctools provides the Python module ``obsutils`` that may further simplify -your analysis efforts. -``obsutils`` is a Python script that makes use of the GammaLib and ctools -modules to create standard analysis steps. -As all Python scripts, ``obsutils`` is part of the cscripts module that -is imported using - -.. code-block:: python - - >>> import cscripts - -Here an example of how to use ``obsutils``: - -.. code-block:: python - - >>> import gammalib - >>> import ctools - >>> from cscripts import obsutils - >>> pattern = obsutils.set_obs_patterns("four", ra=83.63, dec=22.01, offset=1.0) - >>> obs = obsutils.set_obs_list(pattern, duration=1800, emin=0.1, emax=100.0, rad=5.0, caldb="prod2", irf="South_0.5h") - >>> print(obs) - === GObservations === - Number of observations ....: 4 - Number of models ..........: 0 - Number of observed events .: 0 - Number of predicted events : 0 - >>> obs.models(gammalib.GModels("${CTOOLS}/share/models/crab.xml")) - >>> obs = obsutils.sim(obs) - >>> like = obsutils.fit(obs) - >>> print(like.obs().models()) - === GModels === - Number of models ..........: 2 - Number of parameters ......: 10 - ... - -The module is imported using the ``from cscripts import obsutils`` directive. -The ``obsutils.set_obs_patterns()`` function will create a pointing -pattern of four observations located at offset angles of 1 degree from the -nominal location of the Crab nebula. -The ``obsutils.set_obs_list()`` will build an observation container from -that pattern where each pointing will have a duration of 1800 seconds, -cover the 0.1-100 TeV energy range and a field of view of 5°. -The ``South_0.5h`` IRF from the Prod2 calibration database will be used. -A model is then appended to the observation container using the -``obs.models()`` method. -The ``obsutils.sim()`` function then simulates the event data, the -``obsutils.fit()`` function performs a maximum likelihood fit. - -.. note:: - - The ``obsutils`` module is not yet fully developed and more convenience - functions will be added in the future. - - -Access analysis results -~~~~~~~~~~~~~~~~~~~~~~~ - -Here are some examples that show how to access your analysis results in -python. -In the following it is assumed that you have a ctlike object called -``like`` which is setup and runs: - -.. code-block:: python - - import gammalib - import ctools - like = ctools.ctlike(sim.obs()) - like.run() - -* Best-fit parameters: - - The following command shows how to access the fit parameters and errors: - - .. code-block:: python - - obs = like.obs() # Get observations object - obs.models()["Crab"]["Prefactor"].value() # This returns the actual fit value - obs.models()["Crab"]["Prefactor"].error() # This returns the actual fit error - -* Open XML model file from python and print models on screen: - - .. code-block:: python - - models = gammalib.GModels("$CTOOLS/share/models/crab.xml") - print(models) - -* Likelihood value - - .. code-block:: python - - like.opt().value() # Returns likelihood value - -* Curvature Matrix (aka Hessian) - - .. code-block:: python - - curvature = like.opt().curvature().invert() # Return GMatrix object - print(curvature) - - To get the covariance matrix, the curvature matrix needs to be - inverted: - - .. code-block:: python - - covariance = curvature.invert() - print(covariance) - - -Modify and work on XML models -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -* Fix a certain parameter - - .. code-block:: python - - models["Crab"]["Prefactor"].fix() - -* Release a certain parameter - - .. code-block:: python - - models["Crab"]["Prefactor"].free() - -* Check if parameter is free - - .. code-block:: python - - if models["Crab"]["Prefactor"].is_free(): - ... - -* Set parameter value - - .. code-block:: python - - models["Crab"]["Index"].value(-2.5) - -* Set parameter boundaries - - .. code-block:: python - - models["Crab"]["Index"].min(-5.0) - models["Crab"]["Index"].max(-1.0) - - or quicker - - .. code-block:: python - - models["Crab"]["Index"].range(-5.0,-1.0) - -* Loop over models and parameters (e.g. fix all parameters) - - .. code-block:: python - - for model in models: # Loop over models - for par in model: # Loop over parameters - par.fix() # fix parameter - -* Save modified XML model - - .. code-block:: python - - models.save("my_model.xml") - - -Beyond the first steps -~~~~~~~~~~~~~~~~~~~~~~ - -You now have learned the basics of using ctools and GammaLib within Python. -To go beyond these initial steps you may check the Python scripts in the -``examples`` folder that provide useful analysis examples. -Check the ``README.md`` file in that folder for an explanation of the scripts. diff --git a/doc/source/users/user_manual/getting_started/response.rst b/doc/source/users/user_manual/getting_started/response.rst deleted file mode 100644 index 599cd09..0000000 --- a/doc/source/users/user_manual/getting_started/response.rst +++ /dev/null @@ -1,192 +0,0 @@ -.. _sec_response: - -CTA Instrument Response Functions ---------------------------------- - -What are instrument response functions? -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -The instrument response functions provide a mathematical description that -links the measured quantities :math:`\vec{d}` of an event to the physical -quantities :math:`\vec{p}` of the incident photon. The following figure -illustrates this relationship: - -.. figure:: irfs.jpg - :width: 70% - -:math:`I(\vec{p})` is the gamma-ray intensity arriving at Earth as a -function of photon properties :math:`\vec{p}` -(which usually are true photon energy, true photon incident direction, -and true photon arrival time), -while :math:`e(\vec{d})` is the expected event rate as function of event -properties :math:`\vec{d}` (which usually are the measured photon energy, -measured or reconstructued photon incident direction, and measured photon -arrival time). The expected event rate is obtained by integrating the -product of -the instrumental response function :math:`R(\vec{d}|\vec{p},\vec{a}`) -and the emitted intensity :math:`I(\vec{p})` over the photon properties -:math:`\vec{p}`. -The argument :math:`\vec{a}` in the response function comprises any -auxiliary parameter on which the response function may depend on (e.g. -pointing direction, triggered telescopes, optical efficiencies, -atmospheric conditions, etc.). All these quantities and hence the -instrument response function may depend on time. - - -CTA response functions -~~~~~~~~~~~~~~~~~~~~~~ - -The instrument response functions for CTA are factorised into -the effective area :math:`A_{\rm eff}(p, E, t)` (units :math:`cm^2`), -the point spread function :math:`PSF(p' | p, E, t)`, -and the energy dispersion :math:`E_{\rm disp}(E' | p, E, t)` -following: - -.. math:: - R(p', E', t' | p, E, t) = - A_{\rm eff}(p, E, t) \times - 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. -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``. - -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. -In addition, the response files contain an additional table that -describes the background rate as function of energy and position in -the field of view. -An example of a CTA response file is shown below: - -.. figure:: irf-file.png - :width: 100% - -Each table in the response file is in a standardised format that is -the one that is also used for the Fermi/LAT telescope. -As an example, the effective area component of the response file -is shown below. -Response information is stored in a n-dimensional cube, and each axis -of this cube is described by the lower and upper edges of the axis bins. -In this example the effective area is stored as a 2D matrix with the -first axis being energy and the second axis being offaxis angle. -Effective area information is stored for true (``EFFAREA``) and -reconstructed (``EFFAREA_RECO``) energy. -Vector columns are used to store all information. - -.. figure:: irf-aeff.png - :width: 100% - - -Specifying CTA response functions -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -The specification of the CTA Instrument Response Functions depends on the -way how ctools are used. Common to all methods is that the IRFs are -defined by a response name and a calibration database name. -ctools makes use of HEASARC's CALDB format to index and store -IRFs, and specification of the database and response names is -sufficient to access the response. - - -Specifying the response function as input parameters -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -ctools that require instrument response functions have two parameters -to specify the calibration database name and the response function name. -The following example shows a ``ctobssim`` run using the ``prod2`` -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] - -Running the other tools is equivalent. - - -Specifying the response function in an observation definition file -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -In the above example, only a single global response function can be -used for all CTA observations. If you need to specify response functions -per observation you can add the information directly in the XML observation -definition file: - -.. code-block:: xml - - - - - - - - -The ``Calibration`` parameter specifies the calibration database and -response name. You can then pass this file directly to, e.g., ``ctlike``: - -.. code-block:: bash - - $ ctlike - Input event list, counts cube or observation definition XML file [events.fits] obs_irf.xml - Input model XML file [$CTOOLS/share/models/crab.xml] - Output model XML file [crab_results.xml] - -Note that ``ctlike`` does not ask for the calibration database and -response name as it found the relevant information in the XML file. - -.. _sec_cta_rsp_abspath: - - -Specifying individual instrument response files -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -If you need even more control over individual response files, you can -specify them individually in the XML observation file as follows: - -.. code-block:: xml - - - - - - - - - - - - -From within a Python script -^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -The following example illustrates how to set the calibration database -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) - -The calibration database is set by creating a ``GCaldb`` object. The -constructor takes as argument the mission (always ``cta``) and the -database name, in our case ``prod2``. The response function is then set -by passing the response name (here ``South_0.5h``) and the calibration -database object to the ``response`` method. diff --git a/doc/source/users/user_manual/index.rst b/doc/source/users/user_manual/index.rst index 243ca69..993b830 100644 --- a/doc/source/users/user_manual/index.rst +++ b/doc/source/users/user_manual/index.rst @@ -1,4 +1,5 @@ -=========== +.. _sec_user_manual: + User Manual =========== @@ -6,6 +7,7 @@ User Manual :maxdepth: 1 introduction - getting_started/index + response + models plotting - + python diff --git a/doc/source/users/user_manual/introduction.rst b/doc/source/users/user_manual/introduction.rst index 52515ab..03f1d1a 100644 --- a/doc/source/users/user_manual/introduction.rst +++ b/doc/source/users/user_manual/introduction.rst @@ -53,8 +53,7 @@ command line or used as Python classes, as illustrated below: >>> import cscripts >>> pull = cscripts.cspull() -In this User Manual you will now learn how to install and use the ctools -package. +In this User Manual you will now learn how to use the ctools package. If you encounter any problems, please read the :ref:`sec_help` section. And you should read the :ref:`sec_develop` section if you would like to contribute to the ctools developments. diff --git a/doc/source/users/user_manual/irf-aeff.png b/doc/source/users/user_manual/irf-aeff.png new file mode 100644 index 0000000..eabbe69 Binary files /dev/null and b/doc/source/users/user_manual/irf-aeff.png differ diff --git a/doc/source/users/user_manual/irf-file.png b/doc/source/users/user_manual/irf-file.png new file mode 100644 index 0000000..31f1599 Binary files /dev/null and b/doc/source/users/user_manual/irf-file.png differ diff --git a/doc/source/users/user_manual/irfs.jpg b/doc/source/users/user_manual/irfs.jpg new file mode 100644 index 0000000..c60bb63 Binary files /dev/null and b/doc/source/users/user_manual/irfs.jpg differ diff --git a/doc/source/users/user_manual/model_lightcurve.png b/doc/source/users/user_manual/model_lightcurve.png new file mode 100644 index 0000000..27119e0 Binary files /dev/null and b/doc/source/users/user_manual/model_lightcurve.png differ diff --git a/doc/source/users/user_manual/model_phasecurve.png b/doc/source/users/user_manual/model_phasecurve.png new file mode 100644 index 0000000..9ddaa6f Binary files /dev/null and b/doc/source/users/user_manual/model_phasecurve.png differ diff --git a/doc/source/users/user_manual/models.rst b/doc/source/users/user_manual/models.rst new file mode 100644 index 0000000..f9ebdd2 --- /dev/null +++ b/doc/source/users/user_manual/models.rst @@ -0,0 +1,1068 @@ +.. _sec_models: + +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. diff --git a/doc/source/users/user_manual/python.rst b/doc/source/users/user_manual/python.rst new file mode 100644 index 0000000..c530fc9 --- /dev/null +++ b/doc/source/users/user_manual/python.rst @@ -0,0 +1,454 @@ +.. _sec_python: + +Using ctools from Python +------------------------ + +Getting started +~~~~~~~~~~~~~~~ + +ctools provides a Python module that allows using all tools and scripts as +Python classes. +To use ctools from Python all you have to do is to import the ctools +module into Python. +You should also import the GammaLib module, as ctools without GammaLib is +generally not very useful. + +.. code-block:: python + + >>> import gammalib + >>> import ctools + +To illustrate how to use ctools from Python, below is a working example +of an event list simulation using the :ref:`ctobssim` class. +An instance ``sim`` of the :ref:`ctobssim` class is generated and user +parameters are set for this instance using the ``[ ]`` operator. +See the :ref:`reference` for a list of parameters and their types). +The ``execute()`` method executes the :ref:`ctobssim` class in the same +way as if it were executed from the command line. + +.. code-block:: python + + >>> import ctools + >>> sim = ctools.ctobssim() + >>> sim["inmodel"] = "${CTOOLS}/share/models/crab.xml" + >>> sim["outevents"] = "events.fits" + >>> sim["caldb"] = "prod2" + >>> sim["irf"] = "South_0.5h" + >>> sim["ra"] = 83.63 + >>> sim["dec"] = 22.01 + >>> sim["rad"] = 5.0 + >>> sim["tmin"] = 0.0 + >>> sim["tmax"] = 1800.0 + >>> sim["emin"] = 0.1 + >>> sim["emax"] = 100.0 + >>> sim.execute() + +Alternatively, you may "run" the :ref:`ctobssim` tool using + +.. code-block:: python + + >>> sim.run() + +The main difference to the ``execute()`` method is that the ``run()`` +will not write the simulated event file to disk. +Why is this useful? +Well, after having typed ``sim.run()`` the :ref:`ctobssim` class still +exists as an object in memory, including all the simulated events. +The :ref:`ctobssim` class has an ``obs()`` method that returns an +observation container that holds the simulated CTA observation with its +associated events. +To visualise this container, type: + +.. code-block:: python + + >>> print(sim.obs()) + === GObservations === + Number of observations ....: 1 + Number of models ..........: 2 + Number of observed events .: 23099 + Number of predicted events : 0 + +There is one CTA observation in the container and to visualise the events +in that observation you may type: + +.. code-block:: python + + >>> print(sim.obs()[0].events()) + === GCTAEventList === + Number of events ..........: 23099 (disposed in "events.fits") + Time interval .............: 51544.5 - 51544.5 days + === GEbounds === + Number of intervals .......: 1 + Energy range ..............: 100 GeV - 100 TeV + === GCTARoi === + ROI centre ................: RA=83.63, DEC=22.01 [0,0] + ROI radius ................: 5 deg + +The ``obs()[0]`` operator returns the first observation in the observation +container, the ``events()`` operator returns the event list in that +observation. +To see what kind of object you actually got, use: + +.. code-block:: python + + >>> type(sim.obs()[0].events()) + + +The CTA event list is implement by the ``GCTAEventList`` class in the +``cta`` module of GammaLib. +To visualise the individual events you can iterate over the events using a +for loop. +This will show the simulated celestial coordinates (RA, DEC), the +coordinate in the camera system [DETX, DETY], energies and +terrestrial times (TT) of all events. + +.. code-block:: python + + >>> events = sim.obs()[0].events() + >>> for event in events: + ... print(event) + ... + Dir=RA=83.6308, DEC=21.8881 [-0.00212759,1.33661e-05] Energy=106.465 GeV Time=-3.15576e+08 s (TT) + Dir=RA=83.7518, DEC=21.8064 [-0.0035525,0.00197398] Energy=117.706 GeV Time=-3.15576e+08 s (TT) + Dir=RA=83.5545, DEC=22.0933 [0.00145377,-0.00122121] Energy=138.624 GeV Time=-3.15576e+08 s (TT) + ... + + +Performing a maximum likelihood analysis in Python +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +We can now benefit from the fact that we have some simulated events in +memory to fit a model to these events using the :ref:`ctlike` class. +We will do this in unbinned mode. +Here is what you have to do: + +.. code-block:: python + + >>> like = ctools.ctlike(sim.obs()) + >>> like.run() + +This is pretty compact. +Where are the user parameters? +:ref:`ctlike` doesn't in fact need any as all the relevant information is +already contained in the observation container produced by the +:ref:`ctobssim` class. +And you may have recognised that we constructed the :ref:`ctlike` +instance by using the :ref:`ctobssim` observation container as +constructor argument. + +To check how the fit went you can inspect the optimiser used by +:ref:`ctlike`: + +.. code-block:: python + + >>> print(like.opt()) + === GOptimizerLM === + Optimized function value ..: 154553.422 + Absolute precision ........: 0.005 + Acceptable value decrease .: 2 + Optimization status .......: converged + Number of parameters ......: 10 + Number of free parameters .: 4 + Number of iterations ......: 2 + Lambda ....................: 1e-05 + +Apparently, the fit converged after 2 iterations. +Out of 10 parameters in the model 4 have been fitted (the others were kept +fixed). +To inspect the fit results you can print the model container that is a +member of the observation container: + +.. code-block:: python + + >>> print(like.obs().models()) + === GModels === + Number of models ..........: 2 + Number of parameters ......: 10 + === GModelSky === + Name ......................: Crab + Instruments ...............: all + Instrument scale factors ..: unity + Observation identifiers ...: all + Model type ................: PointSource + Model components ..........: "PointSource" * "PowerLaw" * "Constant" + Number of parameters ......: 6 + Number of spatial par's ...: 2 + RA .......................: 83.6331 [-360,360] deg (fixed,scale=1) + DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1) + Number of spectral par's ..: 3 + Prefactor ................: 5.82698e-16 +/- 1.02186e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) + Index ....................: -2.47534 +/- 0.0154764 [-0,-5] (free,scale=-1,gradient) + PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) + Number of temporal par's ..: 1 + Normalization ............: 1 (relative value) (fixed,scale=1,gradient) + === GCTAModelIrfBackground === + Name ......................: CTABackgroundModel + Instruments ...............: CTA + Instrument scale factors ..: unity + Observation identifiers ...: all + Model type ................: "PowerLaw" * "Constant" + Number of parameters ......: 4 + Number of spectral par's ..: 3 + Prefactor ................: 1.01266 +/- 0.0119676 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient) + Index ....................: 0.00474762 +/- 0.00731725 [-5,5] (free,scale=1,gradient) + PivotEnergy ..............: 1e+06 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) + Number of temporal par's ..: 1 + Normalization ............: 1 (relative value) (fixed,scale=1,gradient) + +Suppose you want to repeat the fit by optimising also the position of the +point source. +This is easy from Python: + +.. code-block:: python + + >>> like.obs().models()["Crab"]["RA"].free() + >>> like.obs().models()["Crab"]["DEC"].free() + >>> like.run() + >>> print(like.obs().models()) + ... + RA .......................: 83.6327 +/- 0.000916983 [-360,360] deg (free,scale=1) + DEC ......................: 22.0141 +/- 0.00086378 [-90,90] deg (free,scale=1) + +The ``like.obs().models()`` method provides the model container, using the +``["Crab"]`` operator we access the Crab model in that container and using +the ``["RA"]`` and ``["DEC"]`` methods we access the relevant model +parameters. +The ``free()`` method frees a parameter, the opposite would be a call to +the ``fix()`` method. + +.. warning:: + + Passing observation containers between ctools classes is a very + convenient and powerful way of building in-memory analysis pipelines. + However, this implies that you need some computing ressources when + dealing with large observation containers (for example if you want to + analyse a few 100 hours of data at once). + +.. warning:: + + You have to be aware about the scope of the objects you're + dealing with. + In the above example, the ``sim.obs()`` container is allocated by the + :ref:`ctobssim` class, hence it disappears (a.k.a. goes out of scope) + once the :ref:`ctobssim` class is deleted, as illustrated by the + following example: + + .. code-block:: python + + >>> obs = sim.obs() + >>> del sim + >>> print(obs) + Segmentation fault + + Note that + + .. code-block:: python + + >>> obs = sim.obs() + >>> del sim + >>> print(like.obs()) + === GObservations === + Number of observations ....: 1 + Number of models ..........: 2 + Number of observed events .: 23099 + Number of predicted events : 23099 + + is okay as the ``ctools.ctlike(sim.obs())`` constructor will create + a copy of the observation container that lives within the :ref:`ctlike` + instance. + To preserve an observation container after a ctools object goes out + of scope you have to create a local copy of the container using the + ``copy()`` method: + + .. code-block:: python + + >>> obs = sim.obs().copy() + >>> del sim + >>> print(obs) + === GObservations === + Number of observations ....: 1 + Number of models ..........: 2 + Number of observed events .: 23099 + Number of predicted events : 0 + + +Using obsutils +~~~~~~~~~~~~~~ + +ctools provides the Python module ``obsutils`` that may further simplify +your analysis efforts. +``obsutils`` is a Python script that makes use of the GammaLib and ctools +modules to create standard analysis steps. +As all Python scripts, ``obsutils`` is part of the cscripts module that +is imported using + +.. code-block:: python + + >>> import cscripts + +Here an example of how to use ``obsutils``: + +.. code-block:: python + + >>> import gammalib + >>> import ctools + >>> from cscripts import obsutils + >>> pattern = obsutils.set_obs_patterns("four", ra=83.63, dec=22.01, offset=1.0) + >>> obs = obsutils.set_obs_list(pattern, duration=1800, emin=0.1, emax=100.0, rad=5.0, caldb="prod2", irf="South_0.5h") + >>> print(obs) + === GObservations === + Number of observations ....: 4 + Number of models ..........: 0 + Number of observed events .: 0 + Number of predicted events : 0 + >>> obs.models(gammalib.GModels("${CTOOLS}/share/models/crab.xml")) + >>> obs = obsutils.sim(obs) + >>> like = ctools.ctlike(obs) + >>> like.run() + >>> print(like.obs().models()) + === GModels === + Number of models ..........: 2 + Number of parameters ......: 10 + ... + +The module is imported using the ``from cscripts import obsutils`` directive. +The ``obsutils.set_obs_patterns()`` function will create a pointing +pattern of four observations located at offset angles of 1 degree from the +nominal location of the Crab nebula. +The ``obsutils.set_obs_list()`` will build an observation container from +that pattern where each pointing will have a duration of 1800 seconds, +cover the 0.1-100 TeV energy range and a field of view of 5°. +The ``South_0.5h`` IRF from the Prod2 calibration database will be used. +A model is then appended to the observation container using the +``obs.models()`` method. +The ``obsutils.sim()`` function then simulates the event data and returns +an observation container with the simulated observations. +The observation container is then passed to ``ctlike`` for maximum likelihood +fitting. + + +Access analysis results +~~~~~~~~~~~~~~~~~~~~~~~ + +Here are some examples that show how to access your analysis results in +python. +In the following it is assumed that you have a ctlike object called +``like`` which is setup and runs: + +.. code-block:: python + + import gammalib + import ctools + like = ctools.ctlike(sim.obs()) + like.run() + +* Best-fit parameters: + + The following command shows how to access the fit parameters and errors: + + .. code-block:: python + + obs = like.obs() # Get observations object + obs.models()["Crab"]["Prefactor"].value() # This returns the actual fit value + obs.models()["Crab"]["Prefactor"].error() # This returns the actual fit error + +* Open XML model file from python and print models on screen: + + .. code-block:: python + + models = gammalib.GModels("$CTOOLS/share/models/crab.xml") + print(models) + +* Likelihood value + + .. code-block:: python + + logL = like.opt().value() # Returns log-likelihood value + print(logL) + + or alternatively + + .. code-block:: python + + logL = like.obs().function().value() # Returns log-likelihood value + print(logL) + +* Curvature Matrix (aka Hessian) + + .. code-block:: python + + curvature = like.obs().function().curvature() # Return GMatrixSparse object + print(curvature) + + To get the covariance matrix, the curvature matrix needs to be + inverted: + + .. code-block:: python + + covariance = curvature.invert() + print(covariance) + + +Modify and work on XML models +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +* Fix a certain parameter + + .. code-block:: python + + models["Crab"]["Prefactor"].fix() + +* Release a certain parameter + + .. code-block:: python + + models["Crab"]["Prefactor"].free() + +* Check if parameter is free + + .. code-block:: python + + if models["Crab"]["Prefactor"].is_free(): + ... + +* Set parameter value + + .. code-block:: python + + models["Crab"]["Index"].value(-2.5) + +* Set parameter boundaries + + .. code-block:: python + + models["Crab"]["Index"].min(-5.0) + models["Crab"]["Index"].max(-1.0) + + or quicker + + .. code-block:: python + + models["Crab"]["Index"].range(-5.0,-1.0) + +* Loop over models and parameters (e.g. fix all parameters) + + .. code-block:: python + + for model in models: # Loop over models + for par in model: # Loop over parameters + par.fix() # fix parameter + +* Save modified XML model + + .. code-block:: python + + models.save("my_model.xml") + + +Beyond the first steps +~~~~~~~~~~~~~~~~~~~~~~ + +You now have learned the basics of using ctools and GammaLib within Python. +To go beyond these initial steps you may check the Python scripts in the +``examples`` folder that provide useful analysis examples. +Check the ``README.md`` file in that folder for an explanation of the scripts. diff --git a/doc/source/users/user_manual/response.rst b/doc/source/users/user_manual/response.rst new file mode 100644 index 0000000..599cd09 --- /dev/null +++ b/doc/source/users/user_manual/response.rst @@ -0,0 +1,192 @@ +.. _sec_response: + +CTA Instrument Response Functions +--------------------------------- + +What are instrument response functions? +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The instrument response functions provide a mathematical description that +links the measured quantities :math:`\vec{d}` of an event to the physical +quantities :math:`\vec{p}` of the incident photon. The following figure +illustrates this relationship: + +.. figure:: irfs.jpg + :width: 70% + +:math:`I(\vec{p})` is the gamma-ray intensity arriving at Earth as a +function of photon properties :math:`\vec{p}` +(which usually are true photon energy, true photon incident direction, +and true photon arrival time), +while :math:`e(\vec{d})` is the expected event rate as function of event +properties :math:`\vec{d}` (which usually are the measured photon energy, +measured or reconstructued photon incident direction, and measured photon +arrival time). The expected event rate is obtained by integrating the +product of +the instrumental response function :math:`R(\vec{d}|\vec{p},\vec{a}`) +and the emitted intensity :math:`I(\vec{p})` over the photon properties +:math:`\vec{p}`. +The argument :math:`\vec{a}` in the response function comprises any +auxiliary parameter on which the response function may depend on (e.g. +pointing direction, triggered telescopes, optical efficiencies, +atmospheric conditions, etc.). All these quantities and hence the +instrument response function may depend on time. + + +CTA response functions +~~~~~~~~~~~~~~~~~~~~~~ + +The instrument response functions for CTA are factorised into +the effective area :math:`A_{\rm eff}(p, E, t)` (units :math:`cm^2`), +the point spread function :math:`PSF(p' | p, E, t)`, +and the energy dispersion :math:`E_{\rm disp}(E' | p, E, t)` +following: + +.. math:: + R(p', E', t' | p, E, t) = + A_{\rm eff}(p, E, t) \times + 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. +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``. + +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. +In addition, the response files contain an additional table that +describes the background rate as function of energy and position in +the field of view. +An example of a CTA response file is shown below: + +.. figure:: irf-file.png + :width: 100% + +Each table in the response file is in a standardised format that is +the one that is also used for the Fermi/LAT telescope. +As an example, the effective area component of the response file +is shown below. +Response information is stored in a n-dimensional cube, and each axis +of this cube is described by the lower and upper edges of the axis bins. +In this example the effective area is stored as a 2D matrix with the +first axis being energy and the second axis being offaxis angle. +Effective area information is stored for true (``EFFAREA``) and +reconstructed (``EFFAREA_RECO``) energy. +Vector columns are used to store all information. + +.. figure:: irf-aeff.png + :width: 100% + + +Specifying CTA response functions +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The specification of the CTA Instrument Response Functions depends on the +way how ctools are used. Common to all methods is that the IRFs are +defined by a response name and a calibration database name. +ctools makes use of HEASARC's CALDB format to index and store +IRFs, and specification of the database and response names is +sufficient to access the response. + + +Specifying the response function as input parameters +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +ctools that require instrument response functions have two parameters +to specify the calibration database name and the response function name. +The following example shows a ``ctobssim`` run using the ``prod2`` +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] + +Running the other tools is equivalent. + + +Specifying the response function in an observation definition file +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +In the above example, only a single global response function can be +used for all CTA observations. If you need to specify response functions +per observation you can add the information directly in the XML observation +definition file: + +.. code-block:: xml + + + + + + + + +The ``Calibration`` parameter specifies the calibration database and +response name. You can then pass this file directly to, e.g., ``ctlike``: + +.. code-block:: bash + + $ ctlike + Input event list, counts cube or observation definition XML file [events.fits] obs_irf.xml + Input model XML file [$CTOOLS/share/models/crab.xml] + Output model XML file [crab_results.xml] + +Note that ``ctlike`` does not ask for the calibration database and +response name as it found the relevant information in the XML file. + +.. _sec_cta_rsp_abspath: + + +Specifying individual instrument response files +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +If you need even more control over individual response files, you can +specify them individually in the XML observation file as follows: + +.. code-block:: xml + + + + + + + + + + + + +From within a Python script +^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The following example illustrates how to set the calibration database +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) + +The calibration database is set by creating a ``GCaldb`` object. The +constructor takes as argument the mission (always ``cta``) and the +database name, in our case ``prod2``. The response function is then set +by passing the response name (here ``South_0.5h``) and the calibration +database object to the ``response`` method. -- libgit2 0.21.4