Enhance 1DC tutorial

Jürgen Knödlseder
1 parent 721eb80b
Showing 5 changed files with 134 additions and 25 deletions
doc/source/users/tutorials/1dc/first_detect.rst
doc/source/users/tutorials/1dc/first_fitting.rst
doc/source/users/tutorials/1dc/first_residuals.rst
doc/source/users/tutorials/1dc/first_skymap_fitted.png
doc/source/users/tutorials/1dc/first_skymap_residual.png
@@ -30,6 +30,7 @@ iteratively finds all peaks in the sky map above a given threshold.
 The detection threshold was set here to ``10.0`` which means that only peaks
 are considered that are ten standard deviations above the mean value of the
 sky map pixels.
+
 The script generated also a `ds9 <http://ds9.si.edu>`_ region file that can
 be overlaid on the sky map.
 Due to the relatively high threshold of ten standard deviations, only two
@@ -44,10 +45,13 @@ sources were detected.
 Below is the content of the
 :ref:`model definition XML file <glossary_moddef>`
 that was generated by :ref:`cssrcdetect`.
-In addition to the two point sources the file also contains a background model
-component of type ``CTAIrfBackground``.
-This component was added since the background model type was set to ``IRF``
-in the :ref:`cssrcdetect` run.
+There are three model components in the file and each model component is
+identified using a unique ``name`` attribute.
+The first two components, ``"Src001"`` and ``"Src002"`` are the two point sources
+detected by :ref:`cssrcdetect`, the third component ``"Background"`` is a
+background model component of type ``"CTAIrfBackground"``.
+The latter component was added because the background model type was set to
+``IRF`` in the :ref:`cssrcdetect` run.
 Having a background model component added allows to directly use the
 :ref:`model definition XML file <glossary_moddef>`
 for a model fitting with :ref:`ctlike`.
@@ -3,13 +3,75 @@
 Fitting the model components to the data
 ----------------------------------------
+Now you are ready to fit the model to the data and to determine its maximum
+likelihood parameters.
+You do this with the :ref:`ctlike` tool that adjusts all parameters in the
+:ref:`model definition XML file <glossary_moddef>`
+that have the attribute ``free`` set to ``"1"``.
+In the current example, the free model parameters are the positions and spectral
+parameters of the two point sources and the spectral normalisation of the
+background component.
+You run the :ref:`ctlike` tool as follows:
+
 .. code-block:: bash
    $ ctlike
    Input event list, counts cube or observation definition XML file [events.fits] obs.xml
-   Input model XML file [$CTOOLS/share/models/crab.xml] models.xml
-   Output model XML file [crab_results.xml] results.xml
+   Input model definition XML file [$CTOOLS/share/models/crab.xml] models.xml
+   Output model definition XML file [crab_results.xml] results.xml
+
+The tool will take a few minutes to perform the model fitting, and will write
+the results into an updated
+:ref:`model definition XML file <glossary_moddef>`
+containing the fitted model parameters and their statistical uncertainties.
+You may inspect the log file ``ctlike.log`` to verify that the model fit
+converged properly, as illustrated in the example below:
+
+.. code-block:: bash
+
+   2016-12-02T12:06:03: +=================================+
+   2016-12-02T12:06:03: | Maximum likelihood optimisation |
+   2016-12-02T12:06:03: +=================================+
+   2016-12-02T12:06:19:  >Iteration   0: -logL=13460292.282, Lambda=1.0e-03
+   2016-12-02T12:06:34:  >Iteration   1: -logL=13447430.338, Lambda=1.0e-03, delta=12861.944, max(|grad|)=38643.139247 [Index:13]
+   2016-12-02T12:06:49:  >Iteration   2: -logL=13444219.743, Lambda=1.0e-04, delta=3210.595, max(|grad|)=18223.630095 [RA:0]
+   2016-12-02T12:07:05:  >Iteration   3: -logL=13442740.130, Lambda=1.0e-05, delta=1479.613, max(|grad|)=12560.528406 [RA:0]
+   2016-12-02T12:07:20:  >Iteration   4: -logL=13442178.499, Lambda=1.0e-06, delta=561.631, max(|grad|)=5836.127526 [RA:0]
+   2016-12-02T12:07:35:  >Iteration   5: -logL=13442096.088, Lambda=1.0e-07, delta=82.411, max(|grad|)=1183.057138 [RA:0]
+   2016-12-02T12:07:50:  >Iteration   6: -logL=13442094.149, Lambda=1.0e-08, delta=1.939, max(|grad|)=55.233391 [RA:0]
+   2016-12-02T12:08:06:  >Iteration   7: -logL=13442094.147, Lambda=1.0e-09, delta=0.002, max(|grad|)=-3.893441 [RA:6]
+   2016-12-02T12:08:21:
+   2016-12-02T12:08:21: +=========================================+
+   2016-12-02T12:08:21: | Maximum likelihood optimisation results |
+   2016-12-02T12:08:21: +=========================================+
+   2016-12-02T12:08:21: === GOptimizerLM ===
+   2016-12-02T12:08:21:  Optimized function value ..: 13442094.147
+   2016-12-02T12:08:21:  Absolute precision ........: 0.005
+   2016-12-02T12:08:21:  Acceptable value decrease .: 2
+   2016-12-02T12:08:21:  Optimization status .......: converged
+   2016-12-02T12:08:21:  Number of parameters ......: 16
+   2016-12-02T12:08:21:  Number of free parameters .: 10
+   2016-12-02T12:08:21:  Number of iterations ......: 7
+   2016-12-02T12:08:21:  Lambda ....................: 1e-10
+   2016-12-02T12:08:21:  Maximum log likelihood ....: -13442094.147
+   2016-12-02T12:08:21:  Observed events  (Nobs) ...: 1899265.000
+   2016-12-02T12:08:21:  Predicted events (Npred) ..: 1899225.997 (Nobs - Npred = 39.0025871018879)
+
+You may also convert the fitted model positions into a `ds9 <http://ds9.si.edu>`_
+region file using the :ref:`csmodelinfo` script so that you can overlay the
+fit results over a sky map:
+.. code-block:: bash
+
+   $ csmodelinfo
+   Input model definition XML file [model.xml] results.xml
+   Output DS9 region file [ds9.reg] positions.reg
+
+The following image shows a zoom of the sky map that comprises both point
+sources, with the initial source positions determined by :ref:`cssrcdetect`
+as green crosses and the positions fitted by :ref:`ctlike` as blue circles.
+Obviously, the initial positions were already nearby the fitted position,
+which is required to assure the proper convergence of the fit.
 .. figure:: first_skymap_fitted.png
    :width: 600px
@@ -17,21 +79,41 @@ Fitting the model components to the data
    *Background subtracted sky map of the events recorded during the Galactic Plane Survey around the Galactic Centre with the fitted positions of the sources shown as blue circles*
+You can also convert the spectral parameters of the point sources into a
+butterfly diagram for each source using the :ref:`ctbutterfly` tool.
+The butterfly diagram shows the envelope of all power laws that are
+statistically compatible with the data.
+You create the butterfly diagram for the first source using:
+
 .. code-block:: bash
    $ ctbutterfly
    Input event list, counts cube or observation definition XML file [events.fits] obs.xml
    Source of interest [Crab] Src001
-   Input model XML file [$CTOOLS/share/models/crab.xml] ctlike_unbinned_bkg_3ptsrc.xml
+   Input model definition XML file [$CTOOLS/share/models/crab.xml] results.xml
    Start value for first energy bin in TeV [0.1]
    Stop value for last energy bin in TeV [100.0]
    Output ASCII file [butterfly.txt] butterfly_src001.txt
+You create the butterfly diagram for the second source by selecting ``Src002``
+as ``"Source of interest"``.
+The butterfly diagrams for both sources are displayed in the figure below:
+
+.. figure:: first_spectrum.png
+   :width: 600px
+   :align: center
+
+   *Butterfly diagrams determined with ctbutterfly and spectral points determined with csspec for Src001 (top) and Src002 (bottom)*
+
+The figure also shows spectral points for each source that were determined
+using the :ref:`csspec` script.
+You create the spectrum for the first source using:
+
 .. code-block:: bash
    $ csspec
    Input event list, counts cube, or observation definition XML file [events.fits] obs.xml
-   Input model XML file [results.xml] ctlike_unbinned_bkg_3ptsrc.xml
+   Input model definition XML file [results.xml] results.xml
    Source name [Crab] Src001
    Binning algorithm (LIN|LOG|FILE) [LOG]
    Lower energy limit (TeV) [0.1]
@@ -39,9 +121,9 @@ Fitting the model components to the data
    Number of energy bins (0=unbinned) [10]
    Output spectrum file [spectrum.fits] spectrum_src001.fits
-.. figure:: first_spectrum.png
-   :width: 600px
-   :align: center
-
-   *Butterfly diagrams determined with ctbutterfly and spectral points determined with csspec for Src001 (top) and Src002 (bottom)*
-
+Alike the :ref:`ctbutterfly` tool, you have to run the :ref:`csspec` script
+for each source for which you want to obtain spectral points.
+The script will divide the data into a number of logarithmically spaced energy
+bins and determine the flux in each of the bins using a maximum likelihood
+model fit.
+Ten energy bins were used in this example.
 .. _1dc_first_residuals:
-Inspecting the residuals
-------------------------
+Inspecting the fit residuals
+----------------------------
+
+After doing the model fit you should investigate the residuals to verify that
+the model components properly describe the observed event distribution.
+You do this with the :ref:`csresmap` script by providing the output
+:ref:`model definition XML file <glossary_moddef>`
+``results.xml`` produced by :ref:`ctlike` as input:
 .. code-block:: bash
    $ csresmap
-   Input event list, counts cube, or observation definition XML file [obs.xml]
-   Input model XML file [results.xml]
-   First coordinate of image center in degrees (RA or galactic l) (0-360) [0.0]
-   Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [0.0]
+   Input event list, counts cube, or observation definition XML file [events.fits] obs.xml
+   Input model definition XML file [$CTOOLS/share/models/crab.xml] results.xml
+   First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63] 0.0
+   Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.01] 0.0
    Lower energy limit (TeV) [0.1]
    Upper energy limit (TeV) [100.0]
-   Coordinate System (CEL|GAL) [GAL]
+   Coordinate System (CEL|GAL) [CEL] GAL
    Projection method (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR]
-   Size of the X axis in pixels [200]
-   Size of the Y axis in pixels [200] 100
-   Pixel size (deg/pixel) [0.1] 0.05
-   Residual map computation algorithm (SUB|SUBDIV|SUBDIVSQRT) [SUB]
+   Size of the X axis in pixels [200] 500
+   Size of the Y axis in pixels [200] 250
+   Pixel size (deg/pixel) [0.02] 0.02
+   Residual map computation algorithm (SUB|SUBDIV|SUBDIVSQRT) [SUBDIV] SUB
    Output residual map file [resmap.fits]
+
+This produces a residual map that you can display for example with
+`ds9 <http://ds9.si.edu>`_.
+The figure below shows the result, with the fitted source positions indicated
+by white circles.
+The residuals near the positions of both point sources are relatively flat,
+demonstrating that both excesses are well described by point sources.
+Residuals along the Galactic plane indicate that not all sources that are
+present in the data are described by the model, which is explained by the
+relatively high detection threshold that you used in the source detection
+step.
+
+.. figure:: first_skymap_residual.png
+   :width: 600px
+   :align: center
+
+   *Residual sky map after subtraction of the fitted model*