From c4cb81cc2c76f94ab82c8e9abca9e6f56d0ca533 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Tue, 8 Aug 2017 23:48:40 +0200 Subject: [PATCH] Update quickstart documentation --- doc/source/users/tutorials/quickstart/binning.rst | 82 ++++++++++++++++++++++++++++++++++++++++++++-------------------------------------- doc/source/users/tutorials/quickstart/butterfly.jpg | Bin 60362 -> 0 bytes doc/source/users/tutorials/quickstart/butterfly.rst | 51 --------------------------------------------------- doc/source/users/tutorials/quickstart/dskeys.jpg | Bin 102007 -> 0 bytes doc/source/users/tutorials/quickstart/energy_dispersion.rst | 225 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------------------------------------------------------------------------------------------- doc/source/users/tutorials/quickstart/fitting.rst | 319 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------------------------------ doc/source/users/tutorials/quickstart/index.rst | 16 +++++----------- doc/source/users/tutorials/quickstart/response.rst | 128 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/users/tutorials/quickstart/selecting.rst | 91 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/users/tutorials/quickstart/simulating.rst | 202 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------------------------------------------------------- doc/source/users/tutorials/quickstart/skymap.jpg | Bin 0 -> 188275 bytes doc/source/users/tutorials/quickstart/skymap.rst | 82 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/users/tutorials/quickstart/skymap_subtracted.jpg | Bin 0 -> 130836 bytes doc/source/users/tutorials/quickstart/unbinned.rst | 248 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------------------------------------------------------------------------------------------------------------------------------------------------------- 14 files changed, 842 insertions(+), 602 deletions(-) delete mode 100644 doc/source/users/tutorials/quickstart/butterfly.jpg delete mode 100644 doc/source/users/tutorials/quickstart/butterfly.rst delete mode 100644 doc/source/users/tutorials/quickstart/dskeys.jpg create mode 100644 doc/source/users/tutorials/quickstart/response.rst create mode 100644 doc/source/users/tutorials/quickstart/selecting.rst create mode 100644 doc/source/users/tutorials/quickstart/skymap.jpg create mode 100644 doc/source/users/tutorials/quickstart/skymap.rst create mode 100644 doc/source/users/tutorials/quickstart/skymap_subtracted.jpg diff --git a/doc/source/users/tutorials/quickstart/binning.rst b/doc/source/users/tutorials/quickstart/binning.rst index b0ac1f5..ade56f5 100644 --- a/doc/source/users/tutorials/quickstart/binning.rst +++ b/doc/source/users/tutorials/quickstart/binning.rst @@ -1,40 +1,45 @@ .. _sec_binning_cta: Binning CTA data -~~~~~~~~~~~~~~~~ +---------------- -As next analysis step you will bin the data into a counts cube using -:ref:`ctbin`. -A counts cube is a 3-dimensional data cube, spanned by -Right Ascension (or Galactic longitude), Declination (or Galactic latitude), -and energy (by default logarithmically spaced, but this is under your control). + .. admonition:: What you will learn -:ref:`ctbin` is executed by typing: + You will learn how to use the :ref:`ctbin` tool to **bin the selected + CTA event data into a counts cube**. + + A counts cube is a 3-dimensional data cube that is spanned by Right + Ascension (or Galactic longitude), Declination (or Galactic latitude), + and energy (by default logarithmically spaced, but this is under your + control). + +Let's now bin the selected event data into a counts cube. You do this +using the :ref:`ctbin` tool as follows: .. code-block:: bash - $ ctbin - Input event list or observation definition XML file [events.fits] - First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63] - Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.01] - Projection method (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR] - Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] - Image scale (in degrees/pixel) [0.02] - Size of the X axis in pixels [200] - Size of the Y axis in pixels [200] - Algorithm for defining energy bins (FILE|LIN|LOG) [LOG] - Start value for first energy bin in TeV [0.1] - Stop value for last energy bin in TeV [100.0] - Number of energy bins [20] - Output counts cube file [cntcube.fits] + $ ctbin + Input event list or observation definition XML file [events.fits] selected_events.fits + First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63] + Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.01] + Projection method (AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN) [CAR] + Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] + Image scale (in degrees/pixel) [0.02] + Size of the X axis in pixels [200] + Size of the Y axis in pixels [200] + Algorithm for defining energy bins (FILE|LIN|LOG) [LOG] + Start value for first energy bin in TeV [0.1] + Stop value for last energy bin in TeV [100.0] + Number of energy bins (1-200) [20] + Output counts cube file [cntcube.fits] In this example the events from file ``events.fits`` will be binned into a counts cube stored into the file ``cntcube.fits``. The counts cube -is centred on the location of the Crab nebula (Right Ascension 83.63°, -Declination 22.01°). A cartesian projection aligned in celestial +is centred on the location of the Crab nebula (Right Ascension 83.63 deg, +Declination 22.01 deg). A cartesian projection aligned in celestial coordinates is used and the counts cube has 200 x 200 spatial pixels of 0.02 x 0.02 degrees in size, covering -an area of 4° x 4°, and 20 logarithmically spaced energy bins, covering +an area of 4 deg x 4 deg, and 20 logarithmically spaced energy bins, covering an energy range from 0.1 TeV to 100 TeV. The ``cntcube.fits`` file produced by :ref:`ctbin` contains four extensions: @@ -58,7 +63,7 @@ simulated data. :width: 600px :align: center - *Extensions of a counts cube FITS file* + *Extensions of the counts cube FITS file* An image of the first bin, covering the energy range 100 - 141 GeV, is @@ -68,23 +73,24 @@ shown below: :height: 400px :align: center - *Image of first energy bin of a counts cube* + *Image of first energy bin of the counts cube* For illustration, the last few lines of the log file ``ctbin.log`` are reproduced below: .. code-block:: none - 2016-06-29T10:33:23: +=================+ - 2016-06-29T10:33:23: | Bin observation | - 2016-06-29T10:33:23: +=================+ - 2016-06-29T10:33:23: === CTA observation === - 2016-06-29T10:33:23: Events in list ............: 23099 - 2016-06-29T10:33:23: Events in cube ............: 18532 - 2016-06-29T10:33:23: Event bins outside RoI ....: 0 - 2016-06-29T10:33:23: Events outside cube area ..: 4567 - 2016-06-29T10:33:23: Events outside energy bins : 0 - -From the 23099 events that have been simulated and stored in the -``events.fits`` file, 18532 lie within the cube boundaries and are thus put + 2017-08-08T20:17:06: +=================+ + 2017-08-08T20:17:06: | Bin observation | + 2017-08-08T20:17:06: +=================+ + 2017-08-08T20:17:06: === CTA observation === + 2017-08-08T20:17:06: Events in list ............: 21991 + 2017-08-08T20:17:06: Events in cube ............: 18827 + 2017-08-08T20:17:06: Events outside RoI ........: 0 + 2017-08-08T20:17:06: Events with invalid WCS ...: 0 + 2017-08-08T20:17:06: Events outside cube area ..: 3164 + 2017-08-08T20:17:06: Events outside energy bins : 0 + +From the 21991 events that have been simulated and stored in the +``events.fits`` file, 18827 lie within the cube boundaries and are thus put into the resulting counts cube. diff --git a/doc/source/users/tutorials/quickstart/butterfly.jpg b/doc/source/users/tutorials/quickstart/butterfly.jpg deleted file mode 100644 index 814edc1..0000000 Binary files a/doc/source/users/tutorials/quickstart/butterfly.jpg and /dev/null differ diff --git a/doc/source/users/tutorials/quickstart/butterfly.rst b/doc/source/users/tutorials/quickstart/butterfly.rst deleted file mode 100644 index 4e1b664..0000000 --- a/doc/source/users/tutorials/quickstart/butterfly.rst +++ /dev/null @@ -1,51 +0,0 @@ -.. _sec_butterfly_cta: - -Calculating and visualising a butterfly diagram -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -To visualise the analysis results you can compute a butterfly diagram -which shows the confidence band of the spectral fit using the -:ref:`ctbutterfly` tool. -The butterfly diagram is the envelope of all power law models -that are within a given confidence limit compatible with the data. -By default a confidence level of 68% is used, but the level can be adjusted -using the hidden ``confidence`` parameter. -:ref:`ctbutterfly` computes this envelope by evaluating for each energy the -minimum and maximum intensity of all power law models that fall within the -error ellipse of the prefactor and index parameters. -The error ellipse is derived from the covariance matrix of a maximum -likelihood fit. - -The following example illustrates how you can use the :ref:`ctbutterfly` -tool: - -.. code-block:: bash - - $ ctbutterfly - Input event list, counts cube or observation definition XML file [events.fits] - Calibration database [prod2] - Instrument response function [South_0.5h] - Input model XML file [$CTOOLS/share/models/crab.xml] crab_results.xml - Source of interest [Crab] - 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] - -Now that you have computed the confidence band of the spectral fit and -that you will have an ASCII file named ``butterfly.txt`` on disk you can -visualise the butterfly using the script ``show_butterfly.py`` that is -in the ctools ``example`` folder. You will need matplotlib on your system -to make this work. To launch the script, type: - -.. code-block:: bash - - $CTOOLS/share/examples/python/show_butterfly.py butterfly.txt - -This will result in a canvas which should look like the following: - -.. figure:: butterfly.jpg - :height: 400px - :align: center - - *Confidence band of the power law fit* - \ No newline at end of file diff --git a/doc/source/users/tutorials/quickstart/dskeys.jpg b/doc/source/users/tutorials/quickstart/dskeys.jpg deleted file mode 100644 index 0185466..0000000 Binary files a/doc/source/users/tutorials/quickstart/dskeys.jpg and /dev/null differ diff --git a/doc/source/users/tutorials/quickstart/energy_dispersion.rst b/doc/source/users/tutorials/quickstart/energy_dispersion.rst index 60da875..c024da8 100644 --- a/doc/source/users/tutorials/quickstart/energy_dispersion.rst +++ b/doc/source/users/tutorials/quickstart/energy_dispersion.rst @@ -1,129 +1,132 @@ .. _sec_edisp_cta: -Taking into account energy dispersion -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Taking the energy dispersion into account +----------------------------------------- -You may not have recognised it, but the examples you have exercised so far -have neglected the impact of energy dispersion on the analysis. In reality, + .. admonition:: What you will learn + + You will learn how to **take the energy dispersion (or energy redistribution) + into account**. + + Although the effect of the energy dispersion is often neglegible there + may be cases where you want to consider energy dispersion in an analysis, + for example if you are analysing the data down to very low energies. + + .. warning:: + Energy dispersion is fully implemented in ctools but is computationally + intensive. So be aware that the **tools and scripts will take a substantial + amount of computing time**. + +You may not have recognised, but the examples you have exercised so far +have neglected the impact of the energy dispersion on the analysis. In reality, however, the reconstructed event energy will differ from the true photon energy, and this effect will become particularily important at low energies. There are -therefore cases where you want to take the effects of energy dispersion into -account, and you can do this in all ctools that make use of the instrument -response functions by setting the hidden ``edisp`` parameter to ``yes``. +therefore cases where you want to take the energy dispersion into account. -.. warning:: +To simulate events taking the energy dispersion into account you run the +:ref:`ctobssim` tool with ``edisp=yes`` parameter: - The maximum likelihood computations including energy dispersion are - relatively time consuming, and in many situations the impact of the - energy dispersion on the analysis results will be very small. For this - reason, energy dispersion is by default disabled in ctools. +.. code-block:: bash -To simulate events taking into account the energy dispersion you should -run the :ref:`ctobssim` tool with ``edisp=yes`` parameter: + $ ctobssim edisp=yes + RA of pointing (degrees) (0-360) [83.63] + Dec of pointing (degrees) (-90-90) [22.01] + Radius of FOV (degrees) (0-180) [5.0] + Start time (UTC string, JD, MJD or MET in seconds) [2020-01-01T00:00:00] + Stop time (UTC string, JD, MJD or MET in seconds) [2020-01-01T00:30:00] 2020-01-01T01:00:00 + Lower energy limit (TeV) [0.1] 0.030 + Upper energy limit (TeV) [100.0] 200.0 + Calibration database [prod2] + Instrument response function [South_0.5h] + Input model definition XML file [$CTOOLS/share/models/crab.xml] + Output event data file or observation definition XML file [events.fits] + +You then select events from the simulated data as before using the +:ref:`ctselect` tool: .. code-block:: bash - $ ctobssim edisp=yes - RA of pointing (degrees) (0-360) [83.63] - Dec of pointing (degrees) (-90-90) [22.01] - 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] - -You may then select events as before using the :ref:`ctselect` tool -(since this tool does not need the instrument response function it has no -hidden ``edisp`` parameter): + $ ctselect + Input event list or observation definition XML file [events.fits] + RA for ROI centre (degrees) (0-360) [83.63] + Dec for ROI centre (degrees) (-90-90) [22.01] + Radius of ROI (degrees) (0-180) [3.0] + Start time (UTC string, JD, MJD or MET in seconds) [NONE] 2020-01-01T00:10:00 + Stop time (UTC string, JD, MJD or MET in seconds) [NONE] 2020-01-01T00:40:00 + Lower energy limit (TeV) [0.1] + Upper energy limit (TeV) [100.0] + Output event list or observation definition XML file [selected_events.fits] + +Then you bin as before the selected events into a counts cube using the +:ref:`ctbin` tool: .. code-block:: bash - $ ctselect - Input event list or observation definition XML file [events.fits] - RA for ROI centre (degrees) (0-360) [83.63] - Dec for ROI centre (degrees) (-90-90) [22.01] - Radius of ROI (degrees) (0-180) [3.0] - Start time (CTA MET in seconds) [0.0] - End time (CTA MET in seconds) [0.0] - Lower energy limit (TeV) [0.1] - Upper energy limit (TeV) [100.0] - Output event list or observation definition XML file [selected_events.fits] - -You then perform a maximum likelihood fit using the :ref:`ctlike` tool -with ``edisp=yes`` so that energy dispersion is correctly taken into account -in the likelihood computation: + $ ctbin + Input event list or observation definition XML file [events.fits] selected_events.fits + First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63] + Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.01] + Projection method (AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN) [CAR] + Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] + Image scale (in degrees/pixel) [0.02] + Size of the X axis in pixels [200] + Size of the Y axis in pixels [200] + Algorithm for defining energy bins (FILE|LIN|LOG) [LOG] + Start value for first energy bin in TeV [0.1] + Stop value for last energy bin in TeV [100.0] + Number of energy bins (1-200) [20] + Output counts cube file [cntcube.fits] + +As next step you need to compute the +:ref:`energy dispersion cube ` +using the :ref:`ctexpcube` tool. You run the tool as follows: .. code-block:: bash - $ ctlike edisp=yes - Input event list, counts cube or observation definition XML file [selected_events.fits] - Calibration database [prod2] - Instrument response function [South_0.5h] - Input model XML file [$CTOOLS/share/models/crab.xml] - Output model XML file [crab_results.xml] - -Below you see the output from the ``ctlike.log`` file. - -.. code-block:: none - - 2016-06-29T20:33:36: +=================================+ - 2016-06-29T20:33:36: | Maximum likelihood optimisation | - 2016-06-29T20:33:36: +=================================+ - 2016-06-29T20:33:39: >Iteration 0: -logL=133777.933, Lambda=1.0e-03 - 2016-06-29T20:33:42: >Iteration 1: -logL=133772.125, Lambda=1.0e-03, delta=5.808, max(|grad|)=23.221525 [Index:7] - 2016-06-29T20:33:45: >Iteration 2: -logL=133772.111, Lambda=1.0e-04, delta=0.013, max(|grad|)=-0.042654 [Prefactor:6] - 2016-06-29T20:33:48: >Iteration 3: -logL=133772.111, Lambda=1.0e-05, delta=0.000, max(|grad|)=-0.000611 [Index:3] - ... - 2016-06-29T20:33:50: +=========================================+ - 2016-06-29T20:33:50: | Maximum likelihood optimisation results | - 2016-06-29T20:33:50: +=========================================+ - 2016-06-29T20:33:50: === GOptimizerLM === - 2016-06-29T20:33:50: Optimized function value ..: 133772.111 - 2016-06-29T20:33:50: Absolute precision ........: 0.005 - 2016-06-29T20:33:50: Acceptable value decrease .: 2 - 2016-06-29T20:33:50: Optimization status .......: converged - 2016-06-29T20:33:50: Number of parameters ......: 10 - 2016-06-29T20:33:50: Number of free parameters .: 4 - 2016-06-29T20:33:50: Number of iterations ......: 3 - 2016-06-29T20:33:50: Lambda ....................: 1e-06 - 2016-06-29T20:33:50: Maximum log likelihood ....: -133772.111 - 2016-06-29T20:33:50: Observed events (Nobs) ...: 20952.000 - 2016-06-29T20:33:50: Predicted events (Npred) ..: 20952.000 (Nobs - Npred = 1.01277e-06) - 2016-06-29T20:33:50: === GModels === - 2016-06-29T20:33:50: Number of models ..........: 2 - 2016-06-29T20:33:50: Number of parameters ......: 10 - 2016-06-29T20:33:50: === GModelSky === - 2016-06-29T20:33:50: Name ......................: Crab - 2016-06-29T20:33:50: Instruments ...............: all - 2016-06-29T20:33:50: Instrument scale factors ..: unity - 2016-06-29T20:33:50: Observation identifiers ...: all - 2016-06-29T20:33:50: Model type ................: PointSource - 2016-06-29T20:33:50: Model components ..........: "PointSource" * "PowerLaw" * "Constant" - 2016-06-29T20:33:50: Number of parameters ......: 6 - 2016-06-29T20:33:50: Number of spatial par's ...: 2 - 2016-06-29T20:33:50: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1) - 2016-06-29T20:33:50: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1) - 2016-06-29T20:33:50: Number of spectral par's ..: 3 - 2016-06-29T20:33:50: Prefactor ................: 5.74029e-16 +/- 1.04025e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) - 2016-06-29T20:33:50: Index ....................: -2.49899 +/- 0.015621 [-0,-5] (free,scale=-1,gradient) - 2016-06-29T20:33:50: PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) - 2016-06-29T20:33:50: Number of temporal par's ..: 1 - 2016-06-29T20:33:50: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) - 2016-06-29T20:33:50: === GCTAModelIrfBackground === - 2016-06-29T20:33:50: Name ......................: CTABackgroundModel - 2016-06-29T20:33:50: Instruments ...............: CTA - 2016-06-29T20:33:50: Instrument scale factors ..: unity - 2016-06-29T20:33:50: Observation identifiers ...: all - 2016-06-29T20:33:50: Model type ................: "PowerLaw" * "Constant" - 2016-06-29T20:33:50: Number of parameters ......: 4 - 2016-06-29T20:33:50: Number of spectral par's ..: 3 - 2016-06-29T20:33:50: Prefactor ................: 0.959074 +/- 0.0134496 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient) - 2016-06-29T20:33:50: Index ....................: -0.0162199 +/- 0.00850751 [-5,5] (free,scale=1,gradient) - 2016-06-29T20:33:50: PivotEnergy ..............: 1e+06 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) - 2016-06-29T20:33:50: Number of temporal par's ..: 1 - 2016-06-29T20:33:50: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) + $ ctedispcube + Input event list or observation definition XML file [NONE] selected_events.fits + Calibration database [prod2] + Instrument response function [South_0.5h] + Input counts cube file to extract energy dispersion cube definition [NONE] + First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63] + Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.01] + Projection method (AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN) [CAR] + Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] + Image scale (in degrees/pixel) [1.0] + Size of the X axis in pixels [10] + Size of the Y axis in pixels [10] + Lower energy limit (TeV) [0.1] + Upper energy limit (TeV) [100.0] + Number of energy bins [20] + Output energy dispersion cube file [edispcube.fits] + +Now you are ready to perform a binned maximum likelihood analysis taking the +energy dispersion into account. You do this by running the :ref:`ctlike` tool +with the ``edisp=yes`` parameter. The :ref:`ctlike` tool will now query for the +energy dispersion cube: + +.. code-block:: bash + + $ ctlike edisp=yes + Input event list, counts cube or observation definition XML file [selected_events.fits] cntcube.fits + Input exposure cube file [expcube.fits] + Input PSF cube file [psfcube.fits] + Input energy dispersion cube file [NONE] edispcube.fits + Input background cube file [bkgcube.fits] + Input model definition XML file [$CTOOLS/share/models/crab.xml] models.xml + Output model definition XML file [crab_results.xml] + +You can also perform an unbinned maximum likelihood analysis taking the energy +dispersion into account. In that case the energy dispersion information will be +directly determined from the +:ref:`instrument response functions ` +and no energy dispersion cube is required: + +.. code-block:: bash + $ ctlike edisp=yes + Input event list, counts cube or observation definition XML file [cntcube.fits] selected_events.fits + Calibration database [prod2] + Instrument response function [South_0.5h] + Input model definition XML file [models.xml] $CTOOLS/share/models/crab.xml + Output model definition XML file [crab_results.xml] diff --git a/doc/source/users/tutorials/quickstart/fitting.rst b/doc/source/users/tutorials/quickstart/fitting.rst index dd0994e..2932544 100644 --- a/doc/source/users/tutorials/quickstart/fitting.rst +++ b/doc/source/users/tutorials/quickstart/fitting.rst @@ -1,194 +1,213 @@ .. _sec_fitting_cta: Fitting CTA data -~~~~~~~~~~~~~~~~ +---------------- -Now you are ready to fit the simulated data with a model. For simplicity -the same model is used in the example that you used to simulate the data -with :ref:`ctobssim`. Model fitting is done using the :ref:`ctlike` tool, -and you do the fit by typing: + .. admonition:: What you will learn + + You will learn how to **adjust a parametrised model to the counts cube + by means of a maximum likelihood fit**. + +Now you are ready to fit the data with a model. For simplicity the same model +is used in the example that you used to simulate the data with :ref:`ctobssim`. +Model fitting is done using the :ref:`ctlike` tool, and you do the fit by +typing: .. code-block:: bash - $ ctlike - Input event list, counts cube or observation definition XML file [events.fits] cntcube.fits - Input exposure cube file (only needed for stacked analysis) [NONE] - Input PSF cube file (only needed for stacked analysis) [NONE] - Input background cube file (only needed for stacked analysis) [NONE] - Calibration database [prod2] - Instrument response function [South_0.5h] - Input model XML file [$CTOOLS/share/models/crab.xml] - Output model XML file [crab_results.xml] + $ ctlike + Input event list, counts cube or observation definition XML file [events.fits] cntcube.fits + Input exposure cube file [NONE] expcube.fits + Input PSF cube file [NONE] psfcube.fits + Input background cube file [NONE] bkgcube.fits + Input model definition XML file [$CTOOLS/share/models/crab.xml] models.xml + Output model definition XML file [crab_results.xml] The data are fitted in *binned* mode, which means that the events have been binned into a counts cube and the fit computes the log-likelihood function by summing over all 200 x 200 x 20 bins of the counts cube. There is an alternative method, the so called *unbinned* mode, where the events are not binned into a counts cube and the log-likelihood is computed directly by -summing over all events. We will explore the *unbinned* mode later. - -One of the parameters of :ref:`ctlike` is a source model output file -(we specified ``crab_results.xml`` in the example), and this file will be -a copy of the source model input XML file where the parameter values have -been replaced by the fit results. In addition, the statistical uncertainties -are added for each fitted parameter using the attribute ``error``. -Below we show the XML result file that has been produced by the run: +summing over all events. We will explore the *unbinned* mode later (see +:ref:`here `). + +The :ref:`ctlike` tool produced a +:ref:`model definition XML output file ` +``crab_results.xml`` that contains the fit results. +The file is a copy of the +:ref:`model definition XML input file ` +``models.xml`` where the parameter values were replaced by the fit results. +In addition, the statistical uncertainties were added for each fitted parameter +using the attribute ``error``. +Below the +:ref:`model definition XML output file ` +that was produced by the run: .. code-block:: xml - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + In this example, the ``Prefactor`` and ``Index`` of the spectral model for the Crab nebula as well as the ``Prefactor`` and ``Index`` of the background spectral -model have been fitted (all parameters having the attribute ``free="1"`` are -fitted). Now you see also the effect of having multiplied the background -model with a power law. In that way, the amplitude of the background as -well as it's spectral slope is adjusted by the fit. Obviously, in this -example the adjustment compensates only for the statistical fluctuations -of the background, but with real data, the adjustment may account also for -some of the systematic uncertainties. - -.. warning:: - - As good practice, the amplitude of the background model should always be - left as a free parameter of the fit. Otherwise, any uncertainty in the - background rate will immediately propagate into the flux estimate of the - source. +model were fitted (these parameters had the attribute ``free="1"``). Now you +see also the effect of having multiplied the background model with a power law. +In that way, the amplitude of the background as well as it's spectral slope is +adjusted by the fit. Obviously, in this example the adjustment compensates only +for the statistical fluctuations of the background, but with real data, the +adjustment may account also for some of the systematic uncertainties. + + .. warning:: + + As good practice, the amplitude of the background model should always be + left as a free parameter of the fit. Otherwise, any uncertainty in the + background rate will immediately propagate into the flux estimate of the + source. + + .. warning:: + + You may have recognized the ``scale`` and ``value`` attributes in the + :ref:`model definition XML file `. The value of each + parameter is obtained by multiplying ``value`` with ``scale``. This allows + for a pre-scaling of the parameters, and **you should make use of this + capability to have the value attributes of all parameters that are fitted + of about the same order, typically 1**. This is necessary to assure a + proper convergence of the fitting algorithm. To get more details about the model fitting you can inspect the log file. -Below the last lines of the ctlike.log log file that has been produced by -this run: +Below the last lines of the log file that was produced by this run: .. code-block:: none - 2016-06-29T16:11:05: +=================================+ - 2016-06-29T16:11:05: | Maximum likelihood optimisation | - 2016-06-29T16:11:05: +=================================+ - 2016-06-29T16:11:06: >Iteration 0: -logL=55574.937, Lambda=1.0e-03 - 2016-06-29T16:11:06: >Iteration 1: -logL=55573.383, Lambda=1.0e-03, delta=1.555, max(|grad|)=3.170837 [Index:7] - 2016-06-29T16:11:07: >Iteration 2: -logL=55573.382, Lambda=1.0e-04, delta=0.001, max(|grad|)=0.056280 [Index:3] - ... - 2016-06-29T16:11:08: +=========================================+ - 2016-06-29T16:11:08: | Maximum likelihood optimisation results | - 2016-06-29T16:11:08: +=========================================+ - 2016-06-29T16:11:08: === GOptimizerLM === - 2016-06-29T16:11:08: Optimized function value ..: 55573.382 - 2016-06-29T16:11:08: Absolute precision ........: 0.005 - 2016-06-29T16:11:08: Acceptable value decrease .: 2 - 2016-06-29T16:11:08: Optimization status .......: converged - 2016-06-29T16:11:08: Number of parameters ......: 10 - 2016-06-29T16:11:08: Number of free parameters .: 4 - 2016-06-29T16:11:08: Number of iterations ......: 2 - 2016-06-29T16:11:08: Lambda ....................: 1e-05 - 2016-06-29T16:11:08: Maximum log likelihood ....: -55573.382 - 2016-06-29T16:11:08: Observed events (Nobs) ...: 18532.000 - 2016-06-29T16:11:08: Predicted events (Npred) ..: 18531.998 (Nobs - Npred = 0.00243959) - 2016-06-29T16:11:08: === GModels === - 2016-06-29T16:11:08: Number of models ..........: 2 - 2016-06-29T16:11:08: Number of parameters ......: 10 - 2016-06-29T16:11:08: === GModelSky === - 2016-06-29T16:11:08: Name ......................: Crab - 2016-06-29T16:11:08: Instruments ...............: all - 2016-06-29T16:11:08: Instrument scale factors ..: unity - 2016-06-29T16:11:08: Observation identifiers ...: all - 2016-06-29T16:11:08: Model type ................: PointSource - 2016-06-29T16:11:08: Model components ..........: "PointSource" * "PowerLaw" * "Constant" - 2016-06-29T16:11:08: Number of parameters ......: 6 - 2016-06-29T16:11:08: Number of spatial par's ...: 2 - 2016-06-29T16:11:08: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1) - 2016-06-29T16:11:08: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1) - 2016-06-29T16:11:08: Number of spectral par's ..: 3 - 2016-06-29T16:11:08: Prefactor ................: 5.79064e-16 +/- 1.01828e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) - 2016-06-29T16:11:08: Index ....................: -2.47177 +/- 0.0154935 [-0,-5] (free,scale=-1,gradient) - 2016-06-29T16:11:08: PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) - 2016-06-29T16:11:08: Number of temporal par's ..: 1 - 2016-06-29T16:11:08: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) - 2016-06-29T16:11:08: === GCTAModelIrfBackground === - 2016-06-29T16:11:08: Name ......................: CTABackgroundModel - 2016-06-29T16:11:08: Instruments ...............: CTA - 2016-06-29T16:11:08: Instrument scale factors ..: unity - 2016-06-29T16:11:08: Observation identifiers ...: all - 2016-06-29T16:11:08: Model type ................: "PowerLaw" * "Constant" - 2016-06-29T16:11:08: Number of parameters ......: 4 - 2016-06-29T16:11:08: Number of spectral par's ..: 3 - 2016-06-29T16:11:08: Prefactor ................: 1.01642 +/- 0.0162226 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient) - 2016-06-29T16:11:08: Index ....................: 0.00462945 +/- 0.00964737 [-5,5] (free,scale=1,gradient) - 2016-06-29T16:11:08: PivotEnergy ..............: 1e+06 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) - 2016-06-29T16:11:08: Number of temporal par's ..: 1 - 2016-06-29T16:11:08: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) - -The maximum likelihood optimizer required 2 iterations to converge. This + 2017-08-08T20:24:06: +=================================+ + 2017-08-08T20:24:06: | Maximum likelihood optimisation | + 2017-08-08T20:24:06: +=================================+ + 2017-08-08T20:24:07: >Iteration 0: -logL=56488.484, Lambda=1.0e-03 + 2017-08-08T20:24:07: >Iteration 1: -logL=56431.025, Lambda=1.0e-03, delta=57.459, step=1.0e+00, max(|grad|)=128.772103 [Index:7] + 2017-08-08T20:24:08: >Iteration 2: -logL=56430.794, Lambda=1.0e-04, delta=0.232, step=1.0e+00, max(|grad|)=0.505741 [Index:7] + 2017-08-08T20:24:09: >Iteration 3: -logL=56430.794, Lambda=1.0e-05, delta=0.000, step=1.0e+00, max(|grad|)=0.017423 [Index:3] + 2017-08-08T20:24:09: + 2017-08-08T20:24:09: +=========================================+ + 2017-08-08T20:24:09: | Maximum likelihood optimisation results | + 2017-08-08T20:24:09: +=========================================+ + 2017-08-08T20:24:09: === GOptimizerLM === + 2017-08-08T20:24:09: Optimized function value ..: 56430.794 + 2017-08-08T20:24:09: Absolute precision ........: 0.005 + 2017-08-08T20:24:09: Acceptable value decrease .: 2 + 2017-08-08T20:24:09: Optimization status .......: converged + 2017-08-08T20:24:09: Number of parameters ......: 10 + 2017-08-08T20:24:09: Number of free parameters .: 4 + 2017-08-08T20:24:09: Number of iterations ......: 3 + 2017-08-08T20:24:09: Lambda ....................: 1e-06 + 2017-08-08T20:24:09: Maximum log likelihood ....: -56430.794 + 2017-08-08T20:24:09: Observed events (Nobs) ...: 18827.000 + 2017-08-08T20:24:09: Predicted events (Npred) ..: 18827.000 (Nobs - Npred = 1.92705083463807e-05) + 2017-08-08T20:24:09: === GModels === + 2017-08-08T20:24:09: Number of models ..........: 2 + 2017-08-08T20:24:09: Number of parameters ......: 10 + 2017-08-08T20:24:09: === GModelSky === + 2017-08-08T20:24:09: Name ......................: Crab + 2017-08-08T20:24:09: Instruments ...............: all + 2017-08-08T20:24:09: Instrument scale factors ..: unity + 2017-08-08T20:24:09: Observation identifiers ...: all + 2017-08-08T20:24:09: Model type ................: PointSource + 2017-08-08T20:24:09: Model components ..........: "PointSource" * "PowerLaw" * "Constant" + 2017-08-08T20:24:09: Number of parameters ......: 6 + 2017-08-08T20:24:09: Number of spatial par's ...: 2 + 2017-08-08T20:24:09: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1) + 2017-08-08T20:24:09: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1) + 2017-08-08T20:24:09: Number of spectral par's ..: 3 + 2017-08-08T20:24:09: Prefactor ................: 5.61368377731832e-16 +/- 9.82975197293023e-18 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) + 2017-08-08T20:24:09: Index ....................: -2.45569249143075 +/- 0.0145471260625612 [-0,-5] (free,scale=-1,gradient) + 2017-08-08T20:24:09: PivotEnergy ..............: 300000 [10000,1000000000] MeV (fixed,scale=1000000,gradient) + 2017-08-08T20:24:09: Number of temporal par's ..: 1 + 2017-08-08T20:24:09: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) + 2017-08-08T20:24:09: === GCTAModelCubeBackground === + 2017-08-08T20:24:09: Name ......................: BackgroundModel + 2017-08-08T20:24:09: Instruments ...............: CTA, HESS, MAGIC, VERITAS + 2017-08-08T20:24:09: Instrument scale factors ..: unity + 2017-08-08T20:24:09: Observation identifiers ...: all + 2017-08-08T20:24:09: Model type ................: "PowerLaw" * "Constant" + 2017-08-08T20:24:09: Number of parameters ......: 4 + 2017-08-08T20:24:09: Number of spectral par's ..: 3 + 2017-08-08T20:24:09: Prefactor ................: 0.945994938701748 +/- 0.0148673466828289 [0.01,100] ph/cm2/s/MeV (free,scale=1,gradient) + 2017-08-08T20:24:09: Index ....................: 0.0198073832868029 +/- 0.00948925246912034 [-5,5] (free,scale=1,gradient) + 2017-08-08T20:24:09: PivotEnergy ..............: 1000000 MeV (fixed,scale=1000000,gradient) + 2017-08-08T20:24:09: Number of temporal par's ..: 1 + 2017-08-08T20:24:09: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) + +The maximum likelihood optimizer required 3 iterations to converge. This is pretty fast, but recall that you used the same model file for the simulation and for fitting, hence the initial parameter values were already very close to the best fitting values. To see the impact of the initial parameters on -the fit result, you may re-run :ref:`ctlike` using another copy of the model -XML file where you change the value attributes of the parameters that should be +the fit result, you may re-run :ref:`ctlike` using another copy of the +:ref:`model definition XML input file ` +where you change the value attributes of the parameters that should be fitted. You will see that the optimizer requires a couple of more iterations, but it should converge to the same solution (provided that the initial values are not too far from the best fitting values). -.. note:: + .. note:: - As sanity check you should verify that the predicted number of events - (Npred) is equal to the observed number of events (Nobs). To facilitate - this comparison, :ref:`ctlike` provides the difference Nobs - Npred in - the log file. In real life situations, this difference may not always be - small, in particular if the source model is too constrained. You may - then free some of the model parameters so that the fit can correctly - describe the data. + As sanity check you should verify that the predicted number of events + (Npred) is equal to the observed number of events (Nobs). To facilitate + this comparison, :ref:`ctlike` provides the difference Nobs - Npred in + the log file. In real life situations, this difference may not always be + small, in particular if the source model is too constrained. You may + then free some of the model parameters so that the fit can correctly + describe the data. -.. note:: + .. note:: - The :ref:`ctlike` tool has the ability to estimate the detection - significance for sources in the XML model. This is done by computing - the Test Statistic value which is defined as twice the log-likelihood - difference between fitting a source at a given position on top of a - (background) model or fitting no source. Roughly speaken, the square - root of the Test Statistic value gives the source detection significance - in Gaussian sigmas, although the exact relation depends somewhat on - the formulation of the statistical problem. + The :ref:`ctlike` tool has the ability to estimate the detection + significance for sources in the XML model. This is done by computing + the Test Statistic value which is defined as twice the log-likelihood + difference between fitting a source at a given position on top of a + (background) model or fitting no source. Roughly speaken, the square + root of the Test Statistic value gives the source detection significance + in Gaussian sigmas, although the exact relation depends somewhat on + the formulation of the statistical problem. - To instruct :ref:`ctlike` to compute the Test Statistic value for a - given source you need to add the attribute ``tscalc="1"`` to the XML - file: + To instruct :ref:`ctlike` to compute the Test Statistic value for a + given source you need to add the attribute ``tscalc="1"`` to the XML + file: - .. code-block:: xml + .. code-block:: xml - + - :ref:`ctlike` will then compute the Test Statistic value for that - source and dump the result in the log file: + :ref:`ctlike` will then compute the Test Statistic value for that + source and dump the result in the log file: - .. code-block:: none + .. code-block:: none - 2016-06-29T16:17:38: === GModelSky === - 2016-06-29T16:17:38: Name ......................: Crab - 2016-06-29T16:17:38: Instruments ...............: all - 2016-06-29T16:17:38: Test Statistic ............: 21778.8 + 2017-08-08T20:34:29: === GModelSky === + 2017-08-08T20:34:29: Name ......................: Crab + 2017-08-08T20:34:29: Instruments ...............: all + 2017-08-08T20:34:29: Test Statistic ............: 21873.9558591403 - The Test Statistic value will also be added as new attribute - ``ts`` to the XML result file: + The Test Statistic value will also be added as new attribute + ``ts`` to the XML result file: - .. code-block:: xml + .. code-block:: xml - + diff --git a/doc/source/users/tutorials/quickstart/index.rst b/doc/source/users/tutorials/quickstart/index.rst index ec01bce..b8bf907 100644 --- a/doc/source/users/tutorials/quickstart/index.rst +++ b/doc/source/users/tutorials/quickstart/index.rst @@ -3,23 +3,17 @@ First steps with ctools ----------------------- -This section provides a step-by-step description of a simple CTA data -analysis example using the ctools. In this example we will simulate a -half an hour long CTA observation of the Crab nebula as the single -point source within the CTA field of view. You will learn how to bin the -event data into a counts cube and how to fit a spatial and spectral model -to the data. You will also learn how to analyse directly event lists -without binning the events into a counts cube. We will show how you can -create a butterfly diagram that shows the confidence band of the -fitted spectral model. And you will learn how to take energy dispersion -into account. +This section provides a step-by-step description of a simple CTA data simulation +and analysis example using the ctools. .. toctree:: :maxdepth: 1 simulating + selecting + skymap binning + response fitting unbinned - butterfly energy_dispersion diff --git a/doc/source/users/tutorials/quickstart/response.rst b/doc/source/users/tutorials/quickstart/response.rst new file mode 100644 index 0000000..f7eda4e --- /dev/null +++ b/doc/source/users/tutorials/quickstart/response.rst @@ -0,0 +1,128 @@ +.. _sec_stacked_reponse: + +Pre-computing the response +-------------------------- + + .. admonition:: What you will learn + + You will learn how to **compute the instrument response for a counts + cube**. Pre-computing the instrument response will speed-up + the model fitting later. + +To speed-up the model fitting it is recommended to pre-compute the instrument +response for the counts cube. Specifically, you should compute an exposure +cube, a point spread function cube and a background cube. + +The +:ref:`exposure cube ` +is computed using the :ref:`ctexpcube` tool and provides the effective area +multiplied by the livetime of the observation. +You run :ref:`ctexpcube` as follows: + +.. code-block:: bash + + $ ctexpcube + Input event list or observation definition XML file [NONE] selected_events.fits + Calibration database [prod2] + Instrument response function [South_0.5h] + Input counts cube file to extract exposure cube definition [NONE] cntcube.fits + Output exposure cube file [expcube.fits] + +This produces the FITS file ``expcube.fits`` that contains the exposure as +function of sky position and energy. In the example above you extracted the +binning of the exposure cube from the counts cube, which considerably reduces +the number of parameters that are queried by the tool. + + .. note:: + + The binning of the exposure cube does not need to correspond to the + binning of the counts cube. In any case, exposure values will be + determined by interpolation from the values stored in the exposure cube + file. The same is true for the point spread function and background cubes + that are described below. + +Next, we use the :ref:`ctpsfcube` tool to compute the +:ref:`point spread function (PSF) cube ` for the counts +cube. +You run :ref:`ctpsfcube` as follows: + +.. code-block:: bash + + $ ctpsfcube + Input event list or observation definition XML file [NONE] selected_events.fits + Calibration database [prod2] + Instrument response function [South_0.5h] + Input counts cube file to extract PSF cube definition [NONE] + First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63] + Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.01] + Projection method (AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN) [CAR] + Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] + Image scale (in degrees/pixel) [1.0] + Size of the X axis in pixels [10] + Size of the Y axis in pixels [10] + Lower energy limit (TeV) [0.1] + Upper energy limit (TeV) [100.0] + Number of energy bins [20] + Output PSF cube file [psfcube.fits] + +This produces the FITS file ``psfcube.fits`` that contains the point spread +function as function of sky position and energy. You may have noted in the +example that the definiton of the +:ref:`point spread function cube ` +has not been extracted from the counts cube, since this would lead to a +large FITS file on output. +The point spread function varies in fact only slowly over the field of view +of the camera, and consequently it is sufficient to sample that variation +at a large spatial scale of typically one degree. + +Finally, we use the :ref:`ctbkgcube` tool to compute the +:ref:`background cube `. +You run :ref:`ctbkgcube` as follows: + +.. code-block:: bash + + $ ctbkgcube + Input event list or observation definition XML file [NONE] selected_events.fits + Calibration database [prod2] + Instrument response function [South_0.5h] + Input counts cube file to extract background cube definition [NONE] cntcube.fits + Input model definition XML file [NONE] $CTOOLS/share/models/crab.xml + Output background cube file [bkgcube.fits] + Output model definition XML file [NONE] models.xml + +This produces the FITS file ``bkgcube.fits`` that contains the predicted +background rate as function of sky position and energy. +The tool also produces the +:ref:`model definition file ` +``models.xml`` +on output that will serve as input for the maximum likelihood analysis that +will follow. +The file is a copy of the input +:ref:`model definition file ` +``$CTOOLS/share/models/crab.xml`` +where the input background model has been replaced by a background model of +type ``CTACubeBackground``. Below is the content of the ``models.xml`` file: + +.. code-block:: xml + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc/source/users/tutorials/quickstart/selecting.rst b/doc/source/users/tutorials/quickstart/selecting.rst new file mode 100644 index 0000000..78de090 --- /dev/null +++ b/doc/source/users/tutorials/quickstart/selecting.rst @@ -0,0 +1,91 @@ +.. _sec_selecting_cta: + +Selecting CTA data +------------------ + + .. admonition:: What you will learn + + You will learn how to use the :ref:`ctselect` tool to **select relevant + events from the simulated CTA event data**. + +As next step you may want to select some events from the simulated CTA data. +You may want to restrict the events to a region of interest (ROI) of 3 degrees +around the Crab nebula, you may want to extract a 30 minutes time slice +from the data, and you may want to limit the energy range of the events +to 0.1 - 100 TeV. +You do this event selection using the :ref:`ctselect` tool as follows: + +.. code-block:: bash + + $ ctselect + Input event list or observation definition XML file [events.fits] + RA for ROI centre (degrees) (0-360) [83.63] + Dec for ROI centre (degrees) (-90-90) [22.01] + Radius of ROI (degrees) (0-180) [3.0] + Start time (UTC string, JD, MJD or MET in seconds) [NONE] 2020-01-01T00:10:00 + Stop time (UTC string, JD, MJD or MET in seconds) [NONE] 2020-01-01T00:40:00 + Lower energy limit (TeV) [0.1] + Upper energy limit (TeV) [100.0] + Output event list or observation definition XML file [selected_events.fits] + +This created a new CTA event list FITS file ``selected_events.fits`` +that contains the selected data. +Below an excerpt of the ``ctselect.log`` log file that illustrates the +event selection: + +.. code-block:: none + + 2017-08-08T19:48:43: +===================+ + 2017-08-08T19:48:43: | Input observation | + 2017-08-08T19:48:43: +===================+ + 2017-08-08T19:48:43: === GObservations === + 2017-08-08T19:48:43: Number of observations ....: 1 + 2017-08-08T19:48:43: Number of models ..........: 0 + 2017-08-08T19:48:43: Number of observed events .: 202766 + 2017-08-08T19:48:43: Number of predicted events : 0 + 2017-08-08T19:48:43: + 2017-08-08T19:48:43: +=================+ + 2017-08-08T19:48:43: | Event selection | + 2017-08-08T19:48:43: +=================+ + 2017-08-08T19:48:43: === CTA observation === + 2017-08-08T19:48:43: Input filename ............: events.fits + 2017-08-08T19:48:43: Event extension name ......: EVENTS + 2017-08-08T19:48:43: GTI extension name ........: GTI + 2017-08-08T19:48:43: === Event selection === + 2017-08-08T19:48:43: Time range (MJD) ..........: 58849.0077451852 - 58849.0285785185 days + 2017-08-08T19:48:43: Time range (UTC) ..........: 2020-01-01T00:10:00 - 2020-01-01T00:40:00 + 2017-08-08T19:48:43: Time range (MET) ..........: 631109469.184 - 631111269.184 seconds + 2017-08-08T19:48:43: Selected energy range .....: 0.1 - 100 TeV + 2017-08-08T19:48:43: Requested RoI .............: Centre(RA,DEC)=(83.63, 22.01) deg, Radius=3 deg + 2017-08-08T19:48:43: RoI of data ...............: Centre(RA,DEC)=(83.63, 22.01) deg, Radius=5 deg + 2017-08-08T19:48:43: Selected RoI ..............: Centre(RA,DEC)=(83.63, 22.01) deg, Radius=3 deg + 2017-08-08T19:48:43: cfitsio selection .........: TIME >= 631109469.18399978 && TIME <= 631111269.18400002 && ENERGY >= 0.10000000 && ENERGY <= 100.00000000 && ANGSEP(83.630000,22.010000,RA,DEC) <= 3.000000 + 2017-08-08T19:48:43: FITS filename .............: /var/tmp/tmp.0.FLddl1[EVENTS][TIME >= 631109469.18399978 && TIME <= 631111269.18400002 && ENERGY >= 0.10000000 && ENERGY <= 100.00000000 && ANGSEP(83.630000,22.010000,RA,DEC) <= 3.000000] + 2017-08-08T19:48:43: + 2017-08-08T19:48:43: +====================+ + 2017-08-08T19:48:43: | Output observation | + 2017-08-08T19:48:43: +====================+ + 2017-08-08T19:48:43: === GObservations === + 2017-08-08T19:48:43: Number of observations ....: 1 + 2017-08-08T19:48:43: Number of models ..........: 0 + 2017-08-08T19:48:43: Number of observed events .: 21991 + 2017-08-08T19:48:43: Number of predicted events : 0 + +From the 202766 events in the input event list 21991 events were selected. + + .. note:: + + Event selection by :ref:`ctselect` in a given dimension can be avoided + by specifying ``NONE`` (or ``INDEF``, ``UNDEF`` or ``UNDEFINED``) for + one of the selection parameters. The following run for example only + performs a selection by time but no selection by ROI or energy: + + .. code-block:: bash + + $ ctselect + Input event list or observation definition XML file [events.fits] + RA for ROI centre (degrees) (0-360) [83.63] NONE + Start time (UTC string, JD, MJD or time in seconds) [2020-01-01T00:10:00] + Stop time (UTC string, JD, MJD or time in seconds) [2020-01-01T00:40:00] + Lower energy limit (TeV) [0.1] NONE + Output event list or observation definition XML file [selected_events.fits] diff --git a/doc/source/users/tutorials/quickstart/simulating.rst b/doc/source/users/tutorials/quickstart/simulating.rst index 42076a0..26e9805 100644 --- a/doc/source/users/tutorials/quickstart/simulating.rst +++ b/doc/source/users/tutorials/quickstart/simulating.rst @@ -1,13 +1,26 @@ .. _sec_simulating_cta: Simulating CTA data -~~~~~~~~~~~~~~~~~~~ +------------------- -You can simulate a CTA observation with the :ref:`ctobssim` tool. + .. admonition:: What you will learn + + You will learn how to use the :ref:`ctobssim` tool from the command + line to **simulate CTA event data**. + + :ref:`ctobssim` draws random events using a model of the celestial + gamma-ray intensity distribution that was convolved with the + :ref:`instrument response function `. + The tool adds random events from the expected distribution of the + residual background to the data. + +Let's assume that you want to generate a simulated CTA event list. You +do this by using the :ref:`ctobssim` tool from the command line. To invoke the tool, type :ref:`ctobssim` at the console prompt (which is denoted by ``$``). -:ref:`ctobssim` will query for a number of parameters: +:ref:`ctobssim` will then query for several parameters that define the +characteristics of the simulation: .. code-block:: bash @@ -15,10 +28,10 @@ To invoke the tool, type :ref:`ctobssim` at the console prompt RA of pointing (degrees) (0-360) [83.63] Dec of pointing (degrees) (-90-90) [22.01] 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] + Start time (UTC string, JD, MJD or MET in seconds) [2020-01-01T00:00:00] + Stop time (UTC string, JD, MJD or MET in seconds) [2020-01-01T00:30:00] 2020-01-01T01:00:00 + Lower energy limit (TeV) [0.03] + Upper energy limit (TeV) [200.0] Calibration database [prod2] Instrument response function [South_0.5h] Input model XML file [$CTOOLS/share/models/crab.xml] @@ -28,38 +41,61 @@ Each line represents a query for one parameter value. The line starts with a short description of the parameter, followed by the default value proposed by ctobssim in squared brackets ``[ ]``. -**If no parameter is entered the default value will be used** -(which is the case for all parameters shown here). +**If no parameter value is entered the default value will be used**. Otherwise, the specified value will overwrite the default value. The round brackets ``( )`` indicate the range of possible parameter values (if applicable). + .. note:: -You may have recognised that the environment variable ``$CTOOLS`` has -been used in the path name of the model. ctools will automatically expand -the environment variables in parameter inputs. - -The CTA instrument properties (effective area, PSF width) are taken from -the prod2 database. The response for the southern array using the cuts -optimised for 50 hours of observing time are used. + Times can be entered in various formats in ctools. Times can be provided + as -.. note:: + * UTC date strings (e.g. ``"2020-01-01T00:00:00"``) + * Modified Julian days (e.g. ``"MJD 58849.0"``) + * Julian days (e.g. ``"JD 2458849.5"``) + * Mission elapsed time in seconds, counted from ``2000-01-01T12:00:00`` + (e.g. ``"631108869.18"``) - ctools comes bundled with CTA response functions for the northern and - the southern array. The response functions are based on a Prod2 - analysis with cuts optimised for 0.5 hours, 5 hours and 50 hours of - observing time. The following Instrument Response Functions - (IRFs) are available: ``North_0.5h``, ``North_5h``, ``North_50h``, - ``South_0.5h``, ``South_5h``, and ``South_50h``. + Usually times are given in Terrestial Time (``TT``) but may also be + specified as ``TAI`` or ``UTC`` (e.g. ``"MJD 58849.0 (TAI)"`` or + ``"MJD 58849.0 (UTC)"``). -Events are simulated based on the instrument properties and based on a -source and background model. Only events that fall within the specified -region of interest (ROI), defined as a circle around a sky position in +You may have recognised that the environment variable ``$CTOOLS`` has +been used in the path name of the model. ctools will automatically expand +the environment variables in parameter values. + +The CTA +:ref:`instrument response function ` (IRF) +is taken from the ``prod2`` database. The response for the southern array +using the cuts optimised for 0.5 hours of observing time is used. + + .. note:: + + ctools comes bundled with CTA + :ref:`instrument response functions ` for the northern and + the southern array. The IRFs are based on a ``prod2`` + analysis with cuts optimised for 0.5 hours, 5 hours and 50 hours of + observing time. **Be aware that these times do not need to correspond + to the actual observing time that is simulated.** IRFs optimised for + short observation times correspond to ``loose`` cuts that keep more + photons but also more background events, while IRFs optimised for + long observation times correspond to ``hard`` cuts that limit the + number of background events at the expense of loosing some photons. + The following IRFs are available: + ``North_0.5h``, ``North_5h``, ``North_50h``, + ``South_0.5h``, ``South_5h``, and ``South_50h``. + +Events are simulated based on the +:ref:`instrument response function ` +and based on a source and background model. Only events that fall within the +specified region of interest (ROI), defined as a circle around a sky position in Right Ascension and Declination (in degrees), will be stored in the output -event data file. The duration of the simulation is taken here to 30 minutes -(or 1800 seconds). Events are simulated for energies between 0.1 and 100 TeV. +event data file. The duration of the simulation is taken here to one hour. +Events are simulated for energies between 30 GeV and 200 TeV. -The source and background model is defined by the XML file +The source and background model is defined by the +:ref:`model definition XML file ` ``$CTOOLS/share/models/crab.xml``: .. code-block:: xml @@ -72,7 +108,7 @@ The source and background model is defined by the XML file - + @@ -103,8 +139,10 @@ maximum (``max``) value, and a parameter may be free (``free="1"``) or fixed relevant here for the simulations, but they will be important for the model fitting later. -The spectral intensity I(E) (in units of photons/cm2/s/MeV) of the power law is -given by +The spectral intensity I(E) (in units of +:math:`{\rm photons} \, {\rm cm}^{-2} \, {\rm s}^{-1} \, {\rm MeV}^{-1}`) +of the power law is given by + .. math:: \frac{dN}{dE} = N_0 \left( \frac{E}{E_0} \right)^{\gamma} @@ -115,16 +153,19 @@ where the parameters in the XML definition have the following mappings: * :math:`\gamma` = ``Index`` * :math:`E_0` = ``PivotEnergy`` -.. warning:: +.. - Energies are given in the XML file in MeV units. This is a GammaLib - convention that can not be modified. **So make sure you always use - MeV as energy unit in an XML file.** + .. warning:: + **Energies are given in the XML file in MeV units.** This is a GammaLib + convention that can not be modified. So make sure you always use + MeV as energy unit in an XML file. The instrumental background of CTA is modelled using the background -information provided in the IRF (``CTAIrfBackground``) multipled -by a power law. As it is defined here, the power law represents a -constant of 1, hence the background IRF will be used without any +information provided in the +:ref:`instrument response function ` +(``CTAIrfBackground``), where the energy dependence of the background +model is multipled by a power law. As it is defined here, the power law +represents a constant of 1, hence the background IRF will be used without any modification. The power law will become active when fitting the data later and allows a spectral adjustment of the background model that may account for uncertainties in the background information provided @@ -155,46 +196,45 @@ from this file are shown here: .. code-block:: none - 2016-06-29T10:21:57: +======================+ - 2016-06-29T10:21:57: | Simulate observation | - 2016-06-29T10:21:57: +======================+ - 2016-06-29T10:21:57: === CTA observation === - 2016-06-29T10:21:57: Simulation cone ...........: RA=83.63 deg, Dec=22.01 deg, radius=5.5 deg - 2016-06-29T10:21:57: Time interval .............: 0 - 1800 s - 2016-06-29T10:21:57: Photon energy range .......: 100 GeV - 199.526 GeV - 2016-06-29T10:21:57: Event energy range ........: 100 GeV - 199.526 GeV - 2016-06-29T10:21:57: Simulation area ..........: 5.75561e+09 cm2 - 2016-06-29T10:21:57: Use model ................: Crab - 2016-06-29T10:21:57: Normalization ............: 1 [Crab] - 2016-06-29T10:21:57: Flux .....................: 3.76031e-10 [Crab] photons/cm2/s - 2016-06-29T10:21:57: Normalized flux ..........: 3.76031e-10 [Crab] photons/cm2/s - 2016-06-29T10:21:57: Photon rate ..............: 2.16429 photons/s [Crab] - 2016-06-29T10:21:57: MC source photons ........: 3889 [Crab] - 2016-06-29T10:21:57: MC source events .........: 1226 [Crab] - 2016-06-29T10:21:57: MC source events .........: 1226 (all source models) - 2016-06-29T10:21:57: Photon energy range .......: 199.526 GeV - 398.107 GeV - ... - 2016-06-29T10:21:57: MC source photons .........: 10759 [Crab] - 2016-06-29T10:21:57: MC source events ..........: 3686 [Crab] - 2016-06-29T10:21:57: MC events outside ROI .....: 0 - 2016-06-29T10:21:57: MC background events ......: 19413 - 2016-06-29T10:21:57: MC events .................: 23099 (all models) + 2017-08-08T19:28:54: === CTA observation === + 2017-08-08T19:28:54: Simulation cone ...........: RA=83.63 deg, Dec=22.01 deg, radius=5.5 deg + 2017-08-08T19:28:54: Time interval .............: 6.31109e+08 - 6.31112e+08 s + 2017-08-08T19:28:54: Photon energy range .......: 30 GeV - 72.3622611060088 GeV + 2017-08-08T19:28:54: Event energy range ........: 30 GeV - 72.3622611060088 GeV + 2017-08-08T19:28:54: Simulation area ..........: 1.97769e+09 cm2 + 2017-08-08T19:28:54: Use model ................: Crab + 2017-08-08T19:28:54: Normalization ............: 1 [Crab] + 2017-08-08T19:28:54: Flux .....................: 2.5413e-09 [Crab] photons/cm2/s + 2017-08-08T19:28:54: Normalized flux ..........: 2.5413e-09 [Crab] photons/cm2/s + 2017-08-08T19:28:54: Photon rate ..............: 5.0259 photons/s [Crab] + 2017-08-08T19:28:54: MC source photons ........: 18186 [Crab] + 2017-08-08T19:28:54: MC source events .........: 4175 [Crab] + 2017-08-08T19:28:54: MC source events .........: 4175 (all source models) + 2017-08-08T19:28:54: Photon energy range .......: 72.3622611060088 GeV - 174.543227745807 GeV + ... + 2017-08-08T19:28:54: MC source photons .........: 47371 [Crab] + 2017-08-08T19:28:54: MC source events ..........: 13227 [Crab] + 2017-08-08T19:29:02: MC events outside ROI .....: 0 + 2017-08-08T19:29:02: MC background events ......: 189539 + 2017-08-08T19:29:02: MC identifier 1 ...........: Crab + 2017-08-08T19:29:02: MC identifier 2 ...........: CTABackgroundModel + 2017-08-08T19:29:02: MC events .................: 202766 (all models) Each line starts with the UTC time at which the line has been written. In -this run, 10759 Crab photons have been thrown during a time interval of 1800 -seconds. 3686 of these photons have been registered by CTA as events. In the -same time interval, 19413 background events have been registred by CTA. +this run, 47371 Crab photons have been thrown. 13227 of these photons have been +registered by CTA as events. In the same time interval, 189539 background +events have been registred by CTA. -.. note:: + .. note:: - :ref:`ctobssim` will split the simulated energy range into a number of - slices, controlled via the hidden ``eslices`` parameter (ten energy slices - are used by default). For each energy slice, the simulation area - will be adapted to the effective area of the array in that energy slice, - which helps to keep the computing time low. The log file will provide - information about the simulation in each slice. In the example above, the - simulation results for the first energy slice are shown, followed by a - summary of the results for all slices. + :ref:`ctobssim` will split the simulated energy range into a number of + slices, controlled via the hidden ``eslices`` parameter (ten energy slices + are used by default). For each energy slice, the simulation area + will be adapted to the effective area of the array in that energy slice, + which helps to keep the computing time low. The log file will provide + information about the simulation in each slice. In the example above, the + simulation results for the first energy slice are shown, followed by a + summary of the results for all slices. You may change the name of the log file using the hidden parameter ``logfile``: @@ -222,3 +262,15 @@ the hidden ``debug`` parameter to yes: .. code-block:: bash $ ctobssim debug=yes + +.. + + .. note:: + + All tools have the hidden parameters ``logfile``, ``chatter``, and + ``debug`` and you can use these parameters to control the log file + output. In addition, all tools have the hidden parameter ``clobber`` + that allows to overwrite existing files (set to ``yes`` by default) + and ``mode`` that defines the mode of automatic tool parameters (set to + ``ql`` by default). + diff --git a/doc/source/users/tutorials/quickstart/skymap.jpg b/doc/source/users/tutorials/quickstart/skymap.jpg new file mode 100644 index 0000000..d65a80c Binary files /dev/null and b/doc/source/users/tutorials/quickstart/skymap.jpg differ diff --git a/doc/source/users/tutorials/quickstart/skymap.rst b/doc/source/users/tutorials/quickstart/skymap.rst new file mode 100644 index 0000000..1262c0b --- /dev/null +++ b/doc/source/users/tutorials/quickstart/skymap.rst @@ -0,0 +1,82 @@ +.. _sec_skymap: + +Generating a sky map +-------------------- + + .. admonition:: What you will learn + + You will learn how to use the :ref:`ctskymap` tool to **generate a sky + map from the selected CTA event data**. + +Now let's have a look into the selected data by generating a sky map. You do +this with the :ref:`ctskymap` tool as follows: + +.. code-block:: bash + + $ ctskymap + Input event list or observation definition XML file [events.fits] selected_events.fits + First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63] + Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.01] + Projection method (AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN) [CAR] + Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] + Image scale (in degrees/pixel) [0.02] + Size of the X axis in pixels [200] + Size of the Y axis in pixels [200] + Lower energy limit (TeV) [0.1] + Upper energy limit (TeV) [100.0] + Background subtraction method (NONE|IRF) [NONE] + Output skymap file [skymap.fits] + +The tool produces the file ``skymap.fits`` which contains a sky map of the +events in FITS format. The sky map is centred on the location of the Crab +nebula (Right Ascension 83.63 deg, Declination 22.01 deg) and consists of +200 x 200 spatial pixels of 0.02 x 0.02 degrees in size, covering +an area of 4 deg x 4 deg. + +Below an image of the sky map, displayed using +`ds9 `_ in logarithmic color scale with a small Gaussian +kernel smoothing applied: + +.. figure:: skymap.jpg + :width: 400px + :align: center + + *Sky map of the selected events* + +The sky map shows the Crab nebula on top of a wide-spread distribution of +events that originates from the instrumental background. +The +:ref:`instrument response functions ` +contain templates that describe the spatial and spectral distribution of +the background, and :ref:`ctskymap` can make use of these templates +to subtract the background contribution from the sky map. +You enable the background subtraction by running the :ref:`ctskymap` tool with +the background subtraction method set to ``IRF``, as shown in the following +example: + +.. code-block:: bash + + $ ctskymap + Input event list or observation definition XML file [events.fits] selected_events.fits + First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63] + Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.01] + Projection method (AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN) [CAR] + Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] + Image scale (in degrees/pixel) [0.02] + Size of the X axis in pixels [200] + Size of the Y axis in pixels [200] + Lower energy limit (TeV) [0.1] + Upper energy limit (TeV) [100.0] + Background subtraction method (NONE|IRF) [NONE] IRF + Output skymap file [skymap.fits] skymap_subtracted.fits + +The resulting sky map, displayed using +`ds9 `_ in logarithmic color scale with a small Gaussian +kernel smoothing applied is shown below: + +.. figure:: skymap_subtracted.jpg + :width: 400px + :align: center + + *Background subtracted sky map of the selected events* + diff --git a/doc/source/users/tutorials/quickstart/skymap_subtracted.jpg b/doc/source/users/tutorials/quickstart/skymap_subtracted.jpg new file mode 100644 index 0000000..155eab1 Binary files /dev/null and b/doc/source/users/tutorials/quickstart/skymap_subtracted.jpg differ diff --git a/doc/source/users/tutorials/quickstart/unbinned.rst b/doc/source/users/tutorials/quickstart/unbinned.rst index b2a5452..5654047 100644 --- a/doc/source/users/tutorials/quickstart/unbinned.rst +++ b/doc/source/users/tutorials/quickstart/unbinned.rst @@ -1,128 +1,44 @@ .. _sec_unbinned_cta: Doing an unbinned analysis -~~~~~~~~~~~~~~~~~~~~~~~~~~ +-------------------------- -As gamma-ray events are rare, the counts cubes generated by :ref:`ctbin` -will in general be sparsly populated, having many empty pixels, in -particular at high energies. -An alternative analysis technique consists of working directly on the event -list without binning the events in a counts cube. You will see the benefit of -such an analysis later once you re-run :ref:`ctlike` in unbinned mode. - -For unbinned analysis you first have to define the data space region over -which the analysis is done. This is similiar to the :ref:`ctbin` step for -a binned analysis where you defined the size of the counts cube, the energy -range, and the time interval. For unbinned analysis you have no such thing -as a counts cube, but you have to define over which region of the data space -the selected events are spread (because the ctools have to integrate over -this region to compute the total number of predicted events in the data space -that you analyse). Furthermore, you have to define what energy range is -covered, and what time interval is spanned by the data. All this is done -by the :ref:`ctselect` tool, which replaces the :ref:`ctbin` step in an -unbinned analysis. - -:ref:`ctselect` performs an event selection by choosing only events within -a given region-of-interest (ROI), within a given energy band, and within a -given time interval from the input event list. The ROI is a circular region on -the sky, for which you define the centre (in celestial coordinates) and the -radius. Such a circular ROI is sometimes also called an acceptance cone. The -following example shows how to run ctselect: - -.. code-block:: bash - - $ ctselect - Input event list or observation definition XML file [events.fits] - RA for ROI centre (degrees) (0-360) [83.63] - Dec for ROI centre (degrees) (-90-90) [22.01] - Radius of ROI (degrees) (0-180) [3.0] - Start time (CTA MET in seconds) [0.0] - End time (CTA MET in seconds) [0.0] - Lower energy limit (TeV) [0.1] - Upper energy limit (TeV) [100.0] - Output event list or observation definition XML file [selected_events.fits] - -:ref:`ctselect` takes the input event list ``events.fits``, performs an -event selection, and writes the selected event into the file -``selected_events.fits``. The parameters it will query for are the centre -of the ROI, the radius of the ROI, the start and stop time (in seconds), -and the energy lower and upper limits (in TeV). The event selection information -is also written as a set of data selection keywords to the output events -file ``selected_events.fits``, by respecting the same syntax that has been -implemented for Fermi/LAT. The following image is a screen dump of the data -selection keywords that have been written to the ``EVENTS`` header in the -file ``selected_events.fits``: - -.. figure:: dskeys.jpg - :width: 400px - :align: center + .. admonition:: What you will learn - *Data selection keywords* + You will learn how to **adjust a parametrised model to the events without + binning the data**. -**It is mandatory for an unbinned analysis that these data selection keywords -exist in the FITS file.** -If they don't exist, :ref:`ctlike` will not execute in unbinned mode. + Gamma-ray events are rare, hence the counts cubes generated by + :ref:`ctbin` may be sparsly populated, having many empty pixels, in + particular at high energies. In that case it may be worth to + analyse the events directly with an unbinned maximum likelihood + analysis. -Below some lines of the ``ctselect.log`` file that show the data selection -part: + An unbinned analysis is generally preferred over a binned analysis for + short observation times (a few tens of hours) or if you want to assure + that the analysis results are not biased by the selected binning. Also + the fitting of a light or phase curve is only be possible with an + unbinned maximum likelihood analysis + (see :ref:`light curve <1dc_howto_ligthcurve>` and + :ref:`phase curve <1dc_howto_phasecurve>` fitting). -.. code-block:: none - - 2016-06-29T16:20:10: +=================+ - 2016-06-29T16:20:10: | Event selection | - 2016-06-29T16:20:10: +=================+ - 2016-06-29T16:20:10: === CTA observation === - 2016-06-29T16:20:10: Input filename ............: events.fits - 2016-06-29T16:20:10: Event extension name ......: EVENTS - 2016-06-29T16:20:10: GTI extension name ........: GTI - 2016-06-29T16:20:10: Selected energy range .....: 0.1 - 100 TeV - 2016-06-29T16:20:10: Requested ROI .............: Centre(RA,DEC)=(83.63, 22.01) deg, Radius=3 deg - 2016-06-29T16:20:10: ROI of data ...............: Centre(RA,DEC)=(83.63, 22.01) deg, Radius=5 deg - 2016-06-29T16:20:10: Selected ROI ..............: Centre(RA,DEC)=(83.63, 22.01) deg, Radius=3 deg - 2016-06-29T16:20:10: cfitsio selection .........: ENERGY >= 0.10000000 && ENERGY <= 100.00000000 && ANGSEP(83.630000,22.010000,RA,DEC) <= 3.000000 - 2016-06-29T16:20:10: FITS filename .............: /var/tmp/tmp.0.3aUIXn[EVENTS][ENERGY >= 0.10000000 && ENERGY <= 100.00000000 && ANGSEP(83.630000,22.010000,RA,DEC) <= 3.000000] - -.. note:: - - :ref:`ctobssim` will also write data selection keywords in the event - list FITS file, hence you can run :ref:`ctlike` directly on a FITS file - produced by :ref:`ctobssim`. Any selection performed by :ref:`ctselect` - needs to be fully enclosed within any previous selection, e.g. the ROI - needs to be fully enclosed in the acceptance cone used for event - simulation, the energy selection must be fully comprised in the - range of simulated energies, the same applies for the temporal selection. - :ref:`ctselect` will automatically adjust - the selection parameters to guarantee full enclosure. To keep track of - this adjustment, the :ref:`ctselect` log file quotes the requested - selection, any existing selections, and the selection that was finally - applied. - -.. note:: - - :ref:`ctselect` may of course also be used for event selection prior to - binned analysis, for example to select events for a given period in - time. The ``WEIGHT`` extension of the counts cube will track the overlap - between the ROI applied by :ref:`ctselect` and the counts cube, and this - information will be transparently used in any further analysis to restrict - the analysis on counts cube bins that contain data. - -Now that you have selected the events of interest, you can run ctlike in -unbinned mode. To do this you have to specify the selected event list -instead of the counts cube: +An alternative analysis technique consists of working directly on the event +list without binning the events in a counts cube. +You do this with the :ref:`ctlike` tool as follows: .. code-block:: bash - $ ctlike - Input event list, counts cube or observation definition XML file [cntcube.fits] selected_events.fits - Calibration database [prod2] - Instrument response function [South_0.5h] - Input model XML file [$CTOOLS/share/models/crab.xml] - Output model XML file [crab_results.xml] + $ ctlike + Input event list, counts cube or observation definition XML file [selected_events.fits] + Calibration database [prod2] + Instrument response function [South_0.5h] + Input model definition XML file [models.xml] $CTOOLS/share/models/crab.xml + Output model definition XML file [crab_results.xml] You will recognise that :ref:`ctlike` runs much faster in unbinned mode compared to binned mode. This is understandable as the selected event list contains -only 21604 events, while the binned counts cube you used before had +only 21991 events, while the binned counts cube you used before had 200 x 200 x 20 = 800000 bins. As unbinned maximum likelihood fitting loops over the events (while binned maximum likelihood loops over the bins), there are much less operations to perform in unbinned than in binned mode @@ -130,7 +46,7 @@ there are much less operations to perform in unbinned than in binned mode integrating the models over the region of interest, yet this is negligible compared to the operations needed when looping over all pixels). So as long as you work with small event lists, unbinned mode is faster (this -typically holds up to 50 hours of observing time for point sources). +typically holds up to few tens of hours of observing time). Unbinned :ref:`ctlike` should also be more precise as no binning is performed, hence there is no loss of information due to histogramming. @@ -144,58 +60,58 @@ the same events for both analyses. .. code-block:: none - 2016-06-29T16:24:44: +=================================+ - 2016-06-29T16:24:44: | Maximum likelihood optimisation | - 2016-06-29T16:24:44: +=================================+ - 2016-06-29T16:24:45: >Iteration 0: -logL=137562.094, Lambda=1.0e-03 - 2016-06-29T16:24:45: >Iteration 1: -logL=137560.613, Lambda=1.0e-03, delta=1.481, max(|grad|)=2.548241 [Index:7] - 2016-06-29T16:24:45: >Iteration 2: -logL=137560.613, Lambda=1.0e-04, delta=0.001, max(|grad|)=0.027191 [Index:3] - ... - 2016-06-29T16:24:45: +=========================================+ - 2016-06-29T16:24:45: | Maximum likelihood optimisation results | - 2016-06-29T16:24:45: +=========================================+ - 2016-06-29T16:24:45: === GOptimizerLM === - 2016-06-29T16:24:45: Optimized function value ..: 137560.613 - 2016-06-29T16:24:45: Absolute precision ........: 0.005 - 2016-06-29T16:24:45: Acceptable value decrease .: 2 - 2016-06-29T16:24:45: Optimization status .......: converged - 2016-06-29T16:24:45: Number of parameters ......: 10 - 2016-06-29T16:24:45: Number of free parameters .: 4 - 2016-06-29T16:24:45: Number of iterations ......: 2 - 2016-06-29T16:24:45: Lambda ....................: 1e-05 - 2016-06-29T16:24:45: Maximum log likelihood ....: -137560.613 - 2016-06-29T16:24:45: Observed events (Nobs) ...: 21604.000 - 2016-06-29T16:24:45: Predicted events (Npred) ..: 21603.998 (Nobs - Npred = 0.00202556) - 2016-06-29T16:24:45: === GModels === - 2016-06-29T16:24:45: Number of models ..........: 2 - 2016-06-29T16:24:45: Number of parameters ......: 10 - 2016-06-29T16:24:45: === GModelSky === - 2016-06-29T16:24:45: Name ......................: Crab - 2016-06-29T16:24:45: Instruments ...............: all - 2016-06-29T16:24:45: Instrument scale factors ..: unity - 2016-06-29T16:24:45: Observation identifiers ...: all - 2016-06-29T16:24:45: Model type ................: PointSource - 2016-06-29T16:24:45: Model components ..........: "PointSource" * "PowerLaw" * "Constant" - 2016-06-29T16:24:45: Number of parameters ......: 6 - 2016-06-29T16:24:45: Number of spatial par's ...: 2 - 2016-06-29T16:24:45: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1) - 2016-06-29T16:24:45: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1) - 2016-06-29T16:24:45: Number of spectral par's ..: 3 - 2016-06-29T16:24:45: Prefactor ................: 5.828e-16 +/- 1.0219e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) - 2016-06-29T16:24:45: Index ....................: -2.47546 +/- 0.0154751 [-0,-5] (free,scale=-1,gradient) - 2016-06-29T16:24:45: PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) - 2016-06-29T16:24:45: Number of temporal par's ..: 1 - 2016-06-29T16:24:45: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) - 2016-06-29T16:24:45: === GCTAModelIrfBackground === - 2016-06-29T16:24:45: Name ......................: CTABackgroundModel - 2016-06-29T16:24:45: Instruments ...............: CTA - 2016-06-29T16:24:45: Instrument scale factors ..: unity - 2016-06-29T16:24:45: Observation identifiers ...: all - 2016-06-29T16:24:45: Model type ................: "PowerLaw" * "Constant" - 2016-06-29T16:24:45: Number of parameters ......: 4 - 2016-06-29T16:24:45: Number of spectral par's ..: 3 - 2016-06-29T16:24:45: Prefactor ................: 1.0132 +/- 0.0138017 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient) - 2016-06-29T16:24:45: Index ....................: 0.00653002 +/- 0.00831374 [-5,5] (free,scale=1,gradient) - 2016-06-29T16:24:45: PivotEnergy ..............: 1e+06 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) - 2016-06-29T16:24:45: Number of temporal par's ..: 1 - 2016-06-29T16:24:45: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) + 2017-08-08T20:56:27: +=================================+ + 2017-08-08T20:56:27: | Maximum likelihood optimisation | + 2017-08-08T20:56:27: +=================================+ + 2017-08-08T20:56:27: >Iteration 0: -logL=140091.731, Lambda=1.0e-03 + 2017-08-08T20:56:28: >Iteration 1: -logL=140089.754, Lambda=1.0e-03, delta=1.977, step=1.0e+00, max(|grad|)=5.681267 [Index:3] + 2017-08-08T20:56:28: >Iteration 2: -logL=140089.751, Lambda=1.0e-04, delta=0.003, step=1.0e+00, max(|grad|)=0.255879 [Index:3] + 2017-08-08T20:56:28: + 2017-08-08T20:56:28: +=========================================+ + 2017-08-08T20:56:28: | Maximum likelihood optimisation results | + 2017-08-08T20:56:28: +=========================================+ + 2017-08-08T20:56:28: === GOptimizerLM === + 2017-08-08T20:56:28: Optimized function value ..: 140089.751 + 2017-08-08T20:56:28: Absolute precision ........: 0.005 + 2017-08-08T20:56:28: Acceptable value decrease .: 2 + 2017-08-08T20:56:28: Optimization status .......: converged + 2017-08-08T20:56:28: Number of parameters ......: 10 + 2017-08-08T20:56:28: Number of free parameters .: 4 + 2017-08-08T20:56:28: Number of iterations ......: 2 + 2017-08-08T20:56:28: Lambda ....................: 1e-05 + 2017-08-08T20:56:28: Maximum log likelihood ....: -140089.751 + 2017-08-08T20:56:28: Observed events (Nobs) ...: 21991.000 + 2017-08-08T20:56:28: Predicted events (Npred) ..: 21990.996 (Nobs - Npred = 0.00358908515045187) + 2017-08-08T20:56:28: === GModels === + 2017-08-08T20:56:28: Number of models ..........: 2 + 2017-08-08T20:56:28: Number of parameters ......: 10 + 2017-08-08T20:56:28: === GModelSky === + 2017-08-08T20:56:28: Name ......................: Crab + 2017-08-08T20:56:28: Instruments ...............: all + 2017-08-08T20:56:28: Instrument scale factors ..: unity + 2017-08-08T20:56:28: Observation identifiers ...: all + 2017-08-08T20:56:28: Model type ................: PointSource + 2017-08-08T20:56:28: Model components ..........: "PointSource" * "PowerLaw" * "Constant" + 2017-08-08T20:56:28: Number of parameters ......: 6 + 2017-08-08T20:56:28: Number of spatial par's ...: 2 + 2017-08-08T20:56:28: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1) + 2017-08-08T20:56:28: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1) + 2017-08-08T20:56:28: Number of spectral par's ..: 3 + 2017-08-08T20:56:28: Prefactor ................: 5.69800719463692e-16 +/- 9.95311294944208e-18 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) + 2017-08-08T20:56:28: Index ....................: -2.46101639180337 +/- 0.0145278326558106 [-0,-5] (free,scale=-1,gradient) + 2017-08-08T20:56:28: PivotEnergy ..............: 300000 [10000,1000000000] MeV (fixed,scale=1000000,gradient) + 2017-08-08T20:56:28: Number of temporal par's ..: 1 + 2017-08-08T20:56:28: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) + 2017-08-08T20:56:28: === GCTAModelIrfBackground === + 2017-08-08T20:56:28: Name ......................: CTABackgroundModel + 2017-08-08T20:56:28: Instruments ...............: CTA + 2017-08-08T20:56:28: Instrument scale factors ..: unity + 2017-08-08T20:56:28: Observation identifiers ...: all + 2017-08-08T20:56:28: Model type ................: "PowerLaw" * "Constant" + 2017-08-08T20:56:28: Number of parameters ......: 4 + 2017-08-08T20:56:28: Number of spectral par's ..: 3 + 2017-08-08T20:56:28: Prefactor ................: 1.00310957843794 +/- 0.0134585839142799 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient) + 2017-08-08T20:56:28: Index ....................: 0.00814622034405194 +/- 0.00817635016295375 [-5,5] (free,scale=1,gradient) + 2017-08-08T20:56:28: PivotEnergy ..............: 1000000 [10000,1000000000] MeV (fixed,scale=1000000,gradient) + 2017-08-08T20:56:28: Number of temporal par's ..: 1 + 2017-08-08T20:56:28: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) -- libgit2 0.21.4