diff --git a/doc/source/users/tutorials/1dc/first_fitting.rst b/doc/source/users/tutorials/1dc/first_fitting.rst index 7f6a340..81e925e 100644 --- a/doc/source/users/tutorials/1dc/first_fitting.rst +++ b/doc/source/users/tutorials/1dc/first_fitting.rst @@ -1,10 +1,25 @@ .. _1dc_first_fitting: -Fitting the model components to the data ----------------------------------------- +Fitting the model components to the counts cube +----------------------------------------------- -Now you are ready to fit the model to the data and to determine its maximum -likelihood parameters. + .. admonition:: What you will learn + + You will learn how to **fit a parametric model** to the counts cube using + a maximum likelihood algorithm. + + You will also learn how to display the fit results in form of + `ds9 `_ + region files, butterfly diagrams and spectral points. + + Please note that display tools are **not** part of ctools, yet some + scripts for result display that use the ``matplotlib`` Python module can be + found in the + ``$CTOOLS/share/examples/python`` + folder (see :ref:`how to display results<1dc_howto_display>`). + +Now you are ready to fit the model to the counts cube and to determine its +maximum likelihood parameters. You do this with the :ref:`ctlike` tool that adjusts all parameters in the :ref:`model definition file ` @@ -17,9 +32,12 @@ 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_selected.xml - Input model definition XML file [$CTOOLS/share/models/crab.xml] models.xml - Output model definition XML file [crab_results.xml] results.xml + Input event list, counts cube or observation definition XML file [events.fits] cntcube.fits + Input exposure cube file (only needed for stacked analysis) [NONE] expcube.fits + Input PSF cube file (only needed for stacked analysis) [NONE] psfcube.fits + Input background cube file (only needed for stacked analysis) [NONE] bkgcube.fits + Input model definition XML file [$CTOOLS/share/models/crab.xml] stacked_models.xml + Output model definition XML file [crab_results.xml] stacked_results.xml The tool will take a few minutes (on Mac OS X) to perform the model fitting, and will write the results into an updated @@ -30,45 +48,44 @@ converged properly, as illustrated in the example below: .. code-block:: bash - 2017-05-30T14:04:49: +=================================+ - 2017-05-30T14:04:49: | Maximum likelihood optimisation | - 2017-05-30T14:04:49: +=================================+ - 2017-05-30T14:05:14: >Iteration 0: -logL=22930903.429, Lambda=1.0e-03 - 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] - 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] - 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] - 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] - 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] - 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] - 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] - 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] - 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] - 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] - 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] - 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] - 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] - 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] - 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] - 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] - 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] - 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] - 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] - 2017-05-30T14:12:58: - 2017-05-30T14:12:58: +=========================================+ - 2017-05-30T14:12:58: | Maximum likelihood optimisation results | - 2017-05-30T14:12:58: +=========================================+ - 2017-05-30T14:12:58: === GOptimizerLM === - 2017-05-30T14:12:58: Optimized function value ..: 22902055.599 - 2017-05-30T14:12:58: Absolute precision ........: 0.005 - 2017-05-30T14:12:58: Acceptable value decrease .: 2 - 2017-05-30T14:12:58: Optimization status .......: converged - 2017-05-30T14:12:58: Number of parameters ......: 16 - 2017-05-30T14:12:58: Number of free parameters .: 10 - 2017-05-30T14:12:58: Number of iterations ......: 19 - 2017-05-30T14:12:58: Lambda ....................: 1e-22 - 2017-05-30T14:12:58: Maximum log likelihood ....: -22902055.599 - 2017-05-30T14:12:58: Observed events (Nobs) ...: 3348255.000 - 2017-05-30T14:12:58: Predicted events (Npred) ..: 3348253.996 (Nobs - Npred = 1.00355448853225) + 2017-06-01T10:55:30: +=================================+ + 2017-06-01T10:55:30: | Maximum likelihood optimisation | + 2017-06-01T10:55:30: +=================================+ + 2017-06-01T10:55:44: >Iteration 0: -logL=242348.647, Lambda=1.0e-03 + 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] + 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] + 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] + 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] + 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] + 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] + 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] + 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] + 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] + 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] + 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] + 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] + 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] + 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] + 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] + 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] + 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] + 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] + 2017-06-01T11:00:14: + 2017-06-01T11:00:14: +=========================================+ + 2017-06-01T11:00:14: | Maximum likelihood optimisation results | + 2017-06-01T11:00:14: +=========================================+ + 2017-06-01T11:00:14: === GOptimizerLM === + 2017-06-01T11:00:14: Optimized function value ..: 217423.298 + 2017-06-01T11:00:14: Absolute precision ........: 0.005 + 2017-06-01T11:00:14: Acceptable value decrease .: 2 + 2017-06-01T11:00:14: Optimization status .......: converged + 2017-06-01T11:00:14: Number of parameters ......: 16 + 2017-06-01T11:00:14: Number of free parameters .: 10 + 2017-06-01T11:00:14: Number of iterations ......: 18 + 2017-06-01T11:00:14: Lambda ....................: 1e-21 + 2017-06-01T11:00:14: Maximum log likelihood ....: -217423.298 + 2017-06-01T11:00:14: Observed events (Nobs) ...: 2204717.000 + 2017-06-01T11:00:14: Predicted events (Npred) ..: 2204716.997 (Nobs - Npred = 0.00283915968611836) You may also convert the fitted model positions into a `ds9 `_ region file using the :ref:`csmodelinfo` script so that you can overlay the @@ -77,7 +94,7 @@ fit results over a sky map: .. code-block:: bash $ csmodelinfo pnt_type=circle free_color=black show_labels=no - Input model definition XML file [model.xml] results.xml + Input model definition XML file [model.xml] stacked_results.xml Output DS9 region file [ds9.reg] positions.reg 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. 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 +The butterfly diagram shows the envelope of all spectral models 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_selected.xml - Source of interest [Crab] Src001 - 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 + $ ctbutterfly + Input event list, counts cube or observation definition XML file [events.fits] cntcube.fits + Input exposure cube file (only needed for stacked analysis) [ctexpcube.fits] expcube.fits + Input PSF cube file (only needed for stacked analysis) [psfcube.fits] + Input background cube file (only needed for stacked analysis) [bkgcube.fits] + Source of interest [Crab] Src001 + Input model definition XML file [$CTOOLS/share/models/crab.xml] stacked_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 and for the second source using .. code-block:: bash $ ctbutterfly - Input event list, counts cube or observation definition XML file [events.fits] obs_selected.xml - Source of interest [Crab] Src002 - Input model definition XML file [$CTOOLS/share/models/crab.xml] results.xml + Input event list, counts cube or observation definition XML file [cntcube.fits] + Input exposure cube file (only needed for stacked analysis) [expcube.fits] + Input PSF cube file (only needed for stacked analysis) [psfcube.fits] + Input background cube file (only needed for stacked analysis) [bkgcube.fits] + Source of interest [Src001] Src002 + Input model definition XML file [stacked_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_src002.txt - -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)* + Output ASCII file [butterfly_src001.txt] butterfly_src002.txt +The butterfly diagrams for both sources are displayed in the figure below. 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 @@ -140,8 +156,11 @@ 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_selected.xml - Input model definition XML file [$CTOOLS/share/models/crab.xml] results.xml + Input event list, counts cube, or observation definition XML file [events.fits] cntcube.fits + Input exposure cube file (only needed for stacked analysis) [NONE] expcube.fits + Input PSF cube file (only needed for stacked analysis) [NONE] psfcube.fits + Input background cube file (only needed for stacked analysis) [NONE] bkgcube.fits + Input model definition XML file [$CTOOLS/share/models/crab.xml] stacked_results.xml Source name [Crab] Src001 Binning algorithm (LIN|LOG|FILE) [LOG] Lower energy limit (TeV) [0.1] @@ -154,15 +173,35 @@ and for the second source using .. code-block:: bash $ csspec - Input event list, counts cube, or observation definition XML file [events.fits] obs_selected.xml - Input model definition XML file [$CTOOLS/share/models/crab.xml] results.xml - Source name [Crab] Src002 + Input event list, counts cube, or observation definition XML file [cntcube.fits] + Input exposure cube file (only needed for stacked analysis) [expcube.fits] + Input PSF cube file (only needed for stacked analysis) [psfcube.fits] + Input background cube file (only needed for stacked analysis) [bkgcube.fits] + Input model definition XML file [stacked_results.xml] + Source name [Src001] Src002 Binning algorithm (LIN|LOG|FILE) [LOG] Lower energy limit (TeV) [0.1] Upper energy limit (TeV) [100.0] - Number of energy bins (0=unbinned) [20] 10 - Output spectrum file [spectrum.fits] spectrum_src002.fits + Number of energy bins (0=unbinned) [10] + Output spectrum file [spectrum_src001.fits] spectrum_src002.fits The :ref:`csspec` script divided here the data into ten logarithmically spaced energy bins and determined the source flux in each of the bins using a maximum likelihood model fit. + +.. figure:: first_spectrum_stacked.png + :width: 600px + :align: center + + *Butterfly diagrams determined with ctbutterfly and spectral points determined with csspec for Src001 (red) and Src002 (blue)* + +Obviously, ``Src001`` has a spectral cut-off (red flux points) and hence is not +adequately described by a power law model. +Replacing the spectral model of ``Src001`` by an exponentially cut-off +power law improves the fit to the data, as illustrated by the figure below. + +.. figure:: first_spectrum_cutoff_stacked.png + :width: 600px + :align: center + + *Butterfly diagrams determined with ctbutterfly for an exponentially cut-off power law for Src001 (red)*