From 90577754fd4e6e8ab52e58d21f6b10b297033622 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Mon, 12 Dec 2016 15:26:19 +0100 Subject: [PATCH] Enhance 1DC tutorial --- doc/source/users/tutorials/1dc/first_detect.rst | 12 ++++++++---- doc/source/users/tutorials/1dc/first_fitting.rst | 102 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------- doc/source/users/tutorials/1dc/first_residuals.rst | 45 ++++++++++++++++++++++++++++++++++----------- doc/source/users/tutorials/1dc/first_skymap_fitted.png | Bin 462581 -> 0 bytes doc/source/users/tutorials/1dc/first_skymap_residual.png | Bin 0 -> 436870 bytes 5 files changed, 134 insertions(+), 25 deletions(-) create mode 100644 doc/source/users/tutorials/1dc/first_skymap_residual.png diff --git a/doc/source/users/tutorials/1dc/first_detect.rst b/doc/source/users/tutorials/1dc/first_detect.rst index ac61cdb..dbf6f38 100644 --- a/doc/source/users/tutorials/1dc/first_detect.rst +++ b/doc/source/users/tutorials/1dc/first_detect.rst @@ -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 `_ 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 ` 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 ` for a model fitting with :ref:`ctlike`. diff --git a/doc/source/users/tutorials/1dc/first_fitting.rst b/doc/source/users/tutorials/1dc/first_fitting.rst index 7f3b762..440439d 100644 --- a/doc/source/users/tutorials/1dc/first_fitting.rst +++ b/doc/source/users/tutorials/1dc/first_fitting.rst @@ -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 ` +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 ` +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 `_ +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. diff --git a/doc/source/users/tutorials/1dc/first_residuals.rst b/doc/source/users/tutorials/1dc/first_residuals.rst index 59550af..faa4763 100644 --- a/doc/source/users/tutorials/1dc/first_residuals.rst +++ b/doc/source/users/tutorials/1dc/first_residuals.rst @@ -1,21 +1,44 @@ .. _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 ` +``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 `_. +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* diff --git a/doc/source/users/tutorials/1dc/first_skymap_fitted.png b/doc/source/users/tutorials/1dc/first_skymap_fitted.png index 94fd532..795e71f 100644 Binary files a/doc/source/users/tutorials/1dc/first_skymap_fitted.png and b/doc/source/users/tutorials/1dc/first_skymap_fitted.png differ diff --git a/doc/source/users/tutorials/1dc/first_skymap_residual.png b/doc/source/users/tutorials/1dc/first_skymap_residual.png new file mode 100644 index 0000000..9d612c6 Binary files /dev/null and b/doc/source/users/tutorials/1dc/first_skymap_residual.png differ -- libgit2 0.21.4