Andreas Specovius / ctools

Commit a42feb9f0360adf0059aa093b768953364db0664

Authored by Jürgen Knödlseder
1 parent 224a3bc7

Merge 'update-tutorials'

Inline Side-by-side
Showing 1 changed file with 113 additions and 74 deletions
doc/source/users/tutorials/1dc/first_fitting.rst
  View file @a42feb9
1 1 .. _1dc_first_fitting:
2 2  
3   -Fitting the model components to the data
4   -----------------------------------------
  3 +Fitting the model components to the counts cube
  4 +-----------------------------------------------
5 5  
6   -Now you are ready to fit the model to the data and to determine its maximum
7   -likelihood parameters.
  6 + .. admonition:: What you will learn
  7 +
  8 + You will learn how to **fit a parametric model** to the counts cube using
  9 + a maximum likelihood algorithm.
  10 +
  11 + You will also learn how to display the fit results in form of
  12 + `ds9 <http://ds9.si.edu>`_
  13 + region files, butterfly diagrams and spectral points.
  14 +
  15 + Please note that display tools are **not** part of ctools, yet some
  16 + scripts for result display that use the ``matplotlib`` Python module can be
  17 + found in the
  18 + ``$CTOOLS/share/examples/python``
  19 + folder (see :ref:`how to display results<1dc_howto_display>`).
  20 +
  21 +Now you are ready to fit the model to the counts cube and to determine its
  22 +maximum likelihood parameters.
8 23  
9 24 You do this with the :ref:`ctlike` tool that adjusts all parameters in the
10 25 :ref:`model definition file <glossary_moddef>`
... ... @@ -17,9 +32,12 @@ You run the :ref:`ctlike` tool as follows:
17 32 .. code-block:: bash
18 33  
19 34 $ ctlike
20   - Input event list, counts cube or observation definition XML file [events.fits] obs_selected.xml
21   - Input model definition XML file [$CTOOLS/share/models/crab.xml] models.xml
22   - Output model definition XML file [crab_results.xml] results.xml
  35 + Input event list, counts cube or observation definition XML file [events.fits] cntcube.fits
  36 + Input exposure cube file (only needed for stacked analysis) [NONE] expcube.fits
  37 + Input PSF cube file (only needed for stacked analysis) [NONE] psfcube.fits
  38 + Input background cube file (only needed for stacked analysis) [NONE] bkgcube.fits
  39 + Input model definition XML file [$CTOOLS/share/models/crab.xml] stacked_models.xml
  40 + Output model definition XML file [crab_results.xml] stacked_results.xml
23 41  
24 42 The tool will take a few minutes (on Mac OS X) to perform the model fitting,
25 43 and will write the results into an updated
... ... @@ -30,45 +48,44 @@ converged properly, as illustrated in the example below:
30 48  
31 49 .. code-block:: bash
32 50  
33   - 2017-05-30T14:04:49: +=================================+
34   - 2017-05-30T14:04:49: | Maximum likelihood optimisation |
35   - 2017-05-30T14:04:49: +=================================+
36   - 2017-05-30T14:05:14: >Iteration 0: -logL=22930903.429, Lambda=1.0e-03
37   - 2017-05-30T14:05:38: >Iteration 1: -logL=22906772.079, Lambda=1.0e-03, delta=24131.350, step=1.0e+00, max(|grad|)=53111.535756 [Index:13]
38   - 2017-05-30T14:06:00: >Iteration 2: -logL=22903418.347, Lambda=1.0e-04, delta=3353.732, step=1.0e+00, max(|grad|)=-11638.756860 [RA:0]
39   - 2017-05-30T14:06:23: >Iteration 3: -logL=22902335.786, Lambda=1.0e-05, delta=1082.561, step=1.0e+00, max(|grad|)=-9993.636957 [RA:0]
40   - 2017-05-30T14:06:46: >Iteration 4: -logL=22902144.463, Lambda=1.0e-06, delta=191.323, step=1.0e+00, max(|grad|)=-7053.128713 [RA:0]
41   - 2017-05-30T14:07:09: >Iteration 5: -logL=22902098.908, Lambda=1.0e-07, delta=45.555, step=1.0e+00, max(|grad|)=4767.092710 [RA:6]
42   - 2017-05-30T14:07:33: >Iteration 6: -logL=22902077.092, Lambda=1.0e-08, delta=21.816, step=1.0e+00, max(|grad|)=3512.238312 [RA:6]
43   - 2017-05-30T14:07:56: >Iteration 7: -logL=22902066.319, Lambda=1.0e-09, delta=10.774, step=1.0e+00, max(|grad|)=2564.648439 [RA:6]
44   - 2017-05-30T14:08:21: >Iteration 8: -logL=22902060.969, Lambda=1.0e-10, delta=5.350, step=1.0e+00, max(|grad|)=1861.818754 [RA:6]
45   - 2017-05-30T14:08:45: >Iteration 9: -logL=22902058.299, Lambda=1.0e-11, delta=2.670, step=1.0e+00, max(|grad|)=1346.424825 [RA:6]
46   - 2017-05-30T14:09:08: >Iteration 10: -logL=22902056.961, Lambda=1.0e-12, delta=1.338, step=1.0e+00, max(|grad|)=971.176097 [RA:6]
47   - 2017-05-30T14:09:31: >Iteration 11: -logL=22902056.287, Lambda=1.0e-13, delta=0.674, step=1.0e+00, max(|grad|)=699.274062 [RA:6]
48   - 2017-05-30T14:09:54: >Iteration 12: -logL=22902055.947, Lambda=1.0e-14, delta=0.340, step=1.0e+00, max(|grad|)=502.893171 [RA:6]
49   - 2017-05-30T14:10:18: >Iteration 13: -logL=22902055.774, Lambda=1.0e-15, delta=0.172, step=1.0e+00, max(|grad|)=361.372646 [RA:6]
50   - 2017-05-30T14:10:41: >Iteration 14: -logL=22902055.687, Lambda=1.0e-16, delta=0.088, step=1.0e+00, max(|grad|)=259.519694 [RA:6]
51   - 2017-05-30T14:11:04: >Iteration 15: -logL=22902055.642, Lambda=1.0e-17, delta=0.045, step=1.0e+00, max(|grad|)=186.309625 [RA:6]
52   - 2017-05-30T14:11:27: >Iteration 16: -logL=22902055.620, Lambda=1.0e-18, delta=0.023, step=1.0e+00, max(|grad|)=133.719821 [RA:6]
53   - 2017-05-30T14:11:50: >Iteration 17: -logL=22902055.608, Lambda=1.0e-19, delta=0.012, step=1.0e+00, max(|grad|)=95.962554 [RA:6]
54   - 2017-05-30T14:12:13: >Iteration 18: -logL=22902055.602, Lambda=1.0e-20, delta=0.006, step=1.0e+00, max(|grad|)=68.863372 [RA:6]
55   - 2017-05-30T14:12:35: >Iteration 19: -logL=22902055.599, Lambda=1.0e-21, delta=0.003, step=1.0e+00, max(|grad|)=49.416455 [RA:6]
56   - 2017-05-30T14:12:58:
57   - 2017-05-30T14:12:58: +=========================================+
58   - 2017-05-30T14:12:58: | Maximum likelihood optimisation results |
59   - 2017-05-30T14:12:58: +=========================================+
60   - 2017-05-30T14:12:58: === GOptimizerLM ===
61   - 2017-05-30T14:12:58: Optimized function value ..: 22902055.599
62   - 2017-05-30T14:12:58: Absolute precision ........: 0.005
63   - 2017-05-30T14:12:58: Acceptable value decrease .: 2
64   - 2017-05-30T14:12:58: Optimization status .......: converged
65   - 2017-05-30T14:12:58: Number of parameters ......: 16
66   - 2017-05-30T14:12:58: Number of free parameters .: 10
67   - 2017-05-30T14:12:58: Number of iterations ......: 19
68   - 2017-05-30T14:12:58: Lambda ....................: 1e-22
69   - 2017-05-30T14:12:58: Maximum log likelihood ....: -22902055.599
70   - 2017-05-30T14:12:58: Observed events (Nobs) ...: 3348255.000
71   - 2017-05-30T14:12:58: Predicted events (Npred) ..: 3348253.996 (Nobs - Npred = 1.00355448853225)
  51 + 2017-06-01T10:55:30: +=================================+
  52 + 2017-06-01T10:55:30: | Maximum likelihood optimisation |
  53 + 2017-06-01T10:55:30: +=================================+
  54 + 2017-06-01T10:55:44: >Iteration 0: -logL=242348.647, Lambda=1.0e-03
  55 + 2017-06-01T10:55:58: >Iteration 1: -logL=222180.556, Lambda=1.0e-03, delta=20168.091, step=1.0e+00, max(|grad|)=50061.839794 [Index:13]
  56 + 2017-06-01T10:56:11: >Iteration 2: -logL=218667.504, Lambda=1.0e-04, delta=3513.051, step=1.0e+00, max(|grad|)=-11303.423190 [RA:6]
  57 + 2017-06-01T10:56:24: >Iteration 3: -logL=217668.566, Lambda=1.0e-05, delta=998.938, step=1.0e+00, max(|grad|)=-9596.810651 [RA:6]
  58 + 2017-06-01T10:56:38: >Iteration 4: -logL=217500.656, Lambda=1.0e-06, delta=167.910, step=1.0e+00, max(|grad|)=-6620.600218 [RA:6]
  59 + 2017-06-01T10:56:52: >Iteration 5: -logL=217460.888, Lambda=1.0e-07, delta=39.769, step=1.0e+00, max(|grad|)=-4306.170420 [RA:6]
  60 + 2017-06-01T10:57:05: >Iteration 6: -logL=217441.548, Lambda=1.0e-08, delta=19.339, step=1.0e+00, max(|grad|)=3170.609703 [RA:0]
  61 + 2017-06-01T10:57:20: >Iteration 7: -logL=217432.039, Lambda=1.0e-09, delta=9.510, step=1.0e+00, max(|grad|)=2239.034002 [RA:0]
  62 + 2017-06-01T10:57:34: >Iteration 8: -logL=217427.577, Lambda=1.0e-10, delta=4.461, step=1.0e+00, max(|grad|)=1582.498191 [RA:0]
  63 + 2017-06-01T10:57:48: >Iteration 9: -logL=217425.390, Lambda=1.0e-11, delta=2.188, step=1.0e+00, max(|grad|)=1162.116922 [RA:0]
  64 + 2017-06-01T10:58:03: >Iteration 10: -logL=217424.221, Lambda=1.0e-12, delta=1.169, step=1.0e+00, max(|grad|)=815.765881 [RA:0]
  65 + 2017-06-01T10:58:17: >Iteration 11: -logL=217423.683, Lambda=1.0e-13, delta=0.538, step=1.0e+00, max(|grad|)=542.812704 [RA:0]
  66 + 2017-06-01T10:58:31: >Iteration 12: -logL=217423.482, Lambda=1.0e-14, delta=0.200, step=1.0e+00, max(|grad|)=288.286783 [RA:0]
  67 + 2017-06-01T10:58:45: >Iteration 13: -logL=217423.422, Lambda=1.0e-15, delta=0.060, step=1.0e+00, max(|grad|)=221.413071 [RA:0]
  68 + 2017-06-01T10:59:00: >Iteration 14: -logL=217423.378, Lambda=1.0e-16, delta=0.044, step=1.0e+00, max(|grad|)=226.131427 [RA:0]
  69 + 2017-06-01T10:59:17: >Iteration 15: -logL=217423.330, Lambda=1.0e-17, delta=0.048, step=1.0e+00, max(|grad|)=180.270743 [RA:0]
  70 + 2017-06-01T10:59:32: >Iteration 16: -logL=217423.309, Lambda=1.0e-18, delta=0.020, step=1.0e+00, max(|grad|)=104.148460 [RA:0]
  71 + 2017-06-01T10:59:46: >Iteration 17: -logL=217423.303, Lambda=1.0e-19, delta=0.007, step=1.0e+00, max(|grad|)=63.327107 [RA:0]
  72 + 2017-06-01T11:00:00: >Iteration 18: -logL=217423.298, Lambda=1.0e-20, delta=0.004, step=1.0e+00, max(|grad|)=42.005734 [RA:0]
  73 + 2017-06-01T11:00:14:
  74 + 2017-06-01T11:00:14: +=========================================+
  75 + 2017-06-01T11:00:14: | Maximum likelihood optimisation results |
  76 + 2017-06-01T11:00:14: +=========================================+
  77 + 2017-06-01T11:00:14: === GOptimizerLM ===
  78 + 2017-06-01T11:00:14: Optimized function value ..: 217423.298
  79 + 2017-06-01T11:00:14: Absolute precision ........: 0.005
  80 + 2017-06-01T11:00:14: Acceptable value decrease .: 2
  81 + 2017-06-01T11:00:14: Optimization status .......: converged
  82 + 2017-06-01T11:00:14: Number of parameters ......: 16
  83 + 2017-06-01T11:00:14: Number of free parameters .: 10
  84 + 2017-06-01T11:00:14: Number of iterations ......: 18
  85 + 2017-06-01T11:00:14: Lambda ....................: 1e-21
  86 + 2017-06-01T11:00:14: Maximum log likelihood ....: -217423.298
  87 + 2017-06-01T11:00:14: Observed events (Nobs) ...: 2204717.000
  88 + 2017-06-01T11:00:14: Predicted events (Npred) ..: 2204716.997 (Nobs - Npred = 0.00283915968611836)
72 89  
73 90 You may also convert the fitted model positions into a `ds9 <http://ds9.si.edu>`_
74 91 region file using the :ref:`csmodelinfo` script so that you can overlay the
... ... @@ -77,7 +94,7 @@ fit results over a sky map:
77 94 .. code-block:: bash
78 95  
79 96 $ csmodelinfo pnt_type=circle free_color=black show_labels=no
80   - Input model definition XML file [model.xml] results.xml
  97 + Input model definition XML file [model.xml] stacked_results.xml
81 98 Output DS9 region file [ds9.reg] positions.reg
82 99  
83 100 The command line arguments ``pnt_type``, ``free_color`` and ``show_labels``
... ... @@ -99,40 +116,39 @@ which is required to assure the proper convergence of the fit.
99 116  
100 117 You can also convert the spectral parameters of the point sources into a
101 118 butterfly diagram for each source using the :ref:`ctbutterfly` tool.
102   -The butterfly diagram shows the envelope of all power laws that are
  119 +The butterfly diagram shows the envelope of all spectral models that are
103 120 statistically compatible with the data.
104 121 You create the butterfly diagram for the first source using
105 122  
106 123 .. code-block:: bash
107 124  
108   - $ ctbutterfly
109   - Input event list, counts cube or observation definition XML file [events.fits] obs_selected.xml
110   - Source of interest [Crab] Src001
111   - Input model definition XML file [$CTOOLS/share/models/crab.xml] results.xml
112   - Start value for first energy bin in TeV [0.1]
113   - Stop value for last energy bin in TeV [100.0]
114   - Output ASCII file [butterfly.txt] butterfly_src001.txt
  125 + $ ctbutterfly
  126 + Input event list, counts cube or observation definition XML file [events.fits] cntcube.fits
  127 + Input exposure cube file (only needed for stacked analysis) [ctexpcube.fits] expcube.fits
  128 + Input PSF cube file (only needed for stacked analysis) [psfcube.fits]
  129 + Input background cube file (only needed for stacked analysis) [bkgcube.fits]
  130 + Source of interest [Crab] Src001
  131 + Input model definition XML file [$CTOOLS/share/models/crab.xml] stacked_results.xml
  132 + Start value for first energy bin in TeV [0.1]
  133 + Stop value for last energy bin in TeV [100.0]
  134 + Output ASCII file [butterfly.txt] butterfly_src001.txt
115 135  
116 136 and for the second source using
117 137  
118 138 .. code-block:: bash
119 139  
120 140 $ ctbutterfly
121   - Input event list, counts cube or observation definition XML file [events.fits] obs_selected.xml
122   - Source of interest [Crab] Src002
123   - Input model definition XML file [$CTOOLS/share/models/crab.xml] results.xml
  141 + Input event list, counts cube or observation definition XML file [cntcube.fits]
  142 + Input exposure cube file (only needed for stacked analysis) [expcube.fits]
  143 + Input PSF cube file (only needed for stacked analysis) [psfcube.fits]
  144 + Input background cube file (only needed for stacked analysis) [bkgcube.fits]
  145 + Source of interest [Src001] Src002
  146 + Input model definition XML file [stacked_results.xml]
124 147 Start value for first energy bin in TeV [0.1]
125 148 Stop value for last energy bin in TeV [100.0]
126   - Output ASCII file [butterfly.txt] butterfly_src002.txt
127   -
128   -The butterfly diagrams for both sources are displayed in the figure below:
129   -
130   -.. figure:: first_spectrum.png
131   - :width: 600px
132   - :align: center
133   -
134   - *Butterfly diagrams determined with ctbutterfly and spectral points determined with csspec for Src001 (top) and Src002 (bottom)*
  149 + Output ASCII file [butterfly_src001.txt] butterfly_src002.txt
135 150  
  151 +The butterfly diagrams for both sources are displayed in the figure below.
136 152 The figure also shows spectral points for each source that were determined
137 153 using the :ref:`csspec` script.
138 154 You create the spectrum for the first source using
... ... @@ -140,8 +156,11 @@ You create the spectrum for the first source using
140 156 .. code-block:: bash
141 157  
142 158 $ csspec
143   - Input event list, counts cube, or observation definition XML file [events.fits] obs_selected.xml
144   - Input model definition XML file [$CTOOLS/share/models/crab.xml] results.xml
  159 + Input event list, counts cube, or observation definition XML file [events.fits] cntcube.fits
  160 + Input exposure cube file (only needed for stacked analysis) [NONE] expcube.fits
  161 + Input PSF cube file (only needed for stacked analysis) [NONE] psfcube.fits
  162 + Input background cube file (only needed for stacked analysis) [NONE] bkgcube.fits
  163 + Input model definition XML file [$CTOOLS/share/models/crab.xml] stacked_results.xml
145 164 Source name [Crab] Src001
146 165 Binning algorithm (LIN|LOG|FILE) [LOG]
147 166 Lower energy limit (TeV) [0.1]
... ... @@ -154,15 +173,35 @@ and for the second source using
154 173 .. code-block:: bash
155 174  
156 175 $ csspec
157   - Input event list, counts cube, or observation definition XML file [events.fits] obs_selected.xml
158   - Input model definition XML file [$CTOOLS/share/models/crab.xml] results.xml
159   - Source name [Crab] Src002
  176 + Input event list, counts cube, or observation definition XML file [cntcube.fits]
  177 + Input exposure cube file (only needed for stacked analysis) [expcube.fits]
  178 + Input PSF cube file (only needed for stacked analysis) [psfcube.fits]
  179 + Input background cube file (only needed for stacked analysis) [bkgcube.fits]
  180 + Input model definition XML file [stacked_results.xml]
  181 + Source name [Src001] Src002
160 182 Binning algorithm (LIN|LOG|FILE) [LOG]
161 183 Lower energy limit (TeV) [0.1]
162 184 Upper energy limit (TeV) [100.0]
163   - Number of energy bins (0=unbinned) [20] 10
164   - Output spectrum file [spectrum.fits] spectrum_src002.fits
  185 + Number of energy bins (0=unbinned) [10]
  186 + Output spectrum file [spectrum_src001.fits] spectrum_src002.fits
165 187  
166 188 The :ref:`csspec` script divided here the data into ten logarithmically
167 189 spaced energy bins and determined the source flux in each of the bins using
168 190 a maximum likelihood model fit.
  191 +
  192 +.. figure:: first_spectrum_stacked.png
  193 + :width: 600px
  194 + :align: center
  195 +
  196 + *Butterfly diagrams determined with ctbutterfly and spectral points determined with csspec for Src001 (red) and Src002 (blue)*
  197 +
  198 +Obviously, ``Src001`` has a spectral cut-off (red flux points) and hence is not
  199 +adequately described by a power law model.
  200 +Replacing the spectral model of ``Src001`` by an exponentially cut-off
  201 +power law improves the fit to the data, as illustrated by the figure below.
  202 +
  203 +.. figure:: first_spectrum_cutoff_stacked.png
  204 + :width: 600px
  205 + :align: center
  206 +
  207 + *Butterfly diagrams determined with ctbutterfly for an exponentially cut-off power law for Src001 (red)*
... ...