From c1b6ac1b03ef28ec29f57db5b6366b3b424cec76 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Mon, 7 Dec 2015 23:46:03 +0100 Subject: [PATCH] Update ctools documentation. --- doc/source/download.rst | 5 ++++- doc/source/user_manual/enhancing.rst | 23 ----------------------- doc/source/user_manual/expert.rst | 19 ------------------- doc/source/user_manual/getting.rst | 35 +++++++++++++---------------------- doc/source/user_manual/getting_started/beyond.rst | 18 +++++------------- doc/source/user_manual/getting_started/beyond_combine.rst | 238 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------------------------------------------------------- doc/source/user_manual/getting_started/beyond_irf.rst | 288 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------------------------------------------------------------------------------------------------------------------------- doc/source/user_manual/getting_started/beyond_model.rst | 41 ++++++++++++++++++++++------------------- doc/source/user_manual/getting_started/beyond_residual.rst | 23 ++++++++++++----------- doc/source/user_manual/getting_started/beyond_stacked.rst | 208 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------------------------------------------------------------------------------------- doc/source/user_manual/getting_started/beyond_tsmap.rst | 19 ++++++++++--------- doc/source/user_manual/getting_started/binning.rst | 83 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/user_manual/getting_started/butterfly.rst | 51 +++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/user_manual/getting_started/fitting.rst | 195 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/user_manual/getting_started/quickstart.rst | 771 +++++++++++++++++++-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- doc/source/user_manual/getting_started/resmap-crab.png | Bin 257530 -> 0 bytes doc/source/user_manual/getting_started/response.rst | 19 ++++++++++--------- doc/source/user_manual/getting_started/simulating.rst | 210 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/user_manual/getting_started/tsmap-crab.png | Bin 72360 -> 0 bytes doc/source/user_manual/getting_started/unbinned.rst | 200 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/user_manual/glossary.rst | 2 +- doc/source/user_manual/index.rst | 2 -- doc/source/user_manual/introduction.rst | 11 +++++++++-- doc/source/user_manual/plotting.rst | 56 +++++++++++++++++++++++++++++++++----------------------- 24 files changed, 1259 insertions(+), 1258 deletions(-) delete mode 100644 doc/source/user_manual/enhancing.rst delete mode 100644 doc/source/user_manual/expert.rst create mode 100644 doc/source/user_manual/getting_started/binning.rst create mode 100644 doc/source/user_manual/getting_started/butterfly.rst create mode 100644 doc/source/user_manual/getting_started/fitting.rst create mode 100644 doc/source/user_manual/getting_started/simulating.rst create mode 100644 doc/source/user_manual/getting_started/unbinned.rst diff --git a/doc/source/download.rst b/doc/source/download.rst index 5b2f28d..1a7a805 100644 --- a/doc/source/download.rst +++ b/doc/source/download.rst @@ -13,7 +13,7 @@ bugs. Releases -------- -The latest ctools release is ``ctools-0.10.0`` (1 November 2015). +The latest ctools release is ``ctools-1.0.0`` (6 December 2015). Below a table of ctools releases. Please note that at this stage of the project there is a strict link between the ctools and gammalib versions. @@ -28,6 +28,9 @@ gammalib. * - `ctools `_ - `gammalib `_ - Mac OS X package + * - `1.0.0 `_ + - `1.0.0 `_ + - `ctools-1.0.0-macosx10.3.dmg `_ * - `0.10.0 `_ - `0.11.0 `_ - `ctools-0.10.0-macosx10.3.dmg `_ diff --git a/doc/source/user_manual/enhancing.rst b/doc/source/user_manual/enhancing.rst deleted file mode 100644 index 42169c4..0000000 --- a/doc/source/user_manual/enhancing.rst +++ /dev/null @@ -1,23 +0,0 @@ -Enhancing ctools -================ - -Most ctools users will run the existing ctools or cscripts, -which cover the most common simulation and analysis needs. - -Advanced users and developers can write their own tools and scripts. -The following C++ and Python examples can help get you started: - -* :doxy:`ctlike` -* :doxy:`obsutils` -* :doxy:`cstsdist` - -Building your own cscripts --------------------------- - -To be written ... - - -Building your own ctools ------------------------- - -To be written ... diff --git a/doc/source/user_manual/expert.rst b/doc/source/user_manual/expert.rst deleted file mode 100644 index 61b6ee9..0000000 --- a/doc/source/user_manual/expert.rst +++ /dev/null @@ -1,19 +0,0 @@ -Behind the scenes -================= - -CTA data structure ------------------- - -To be written ... - - -Instrument response functions ------------------------------ - -To be written ... - - -Maximum likelihood fitting --------------------------- - -To be written ... diff --git a/doc/source/user_manual/getting.rst b/doc/source/user_manual/getting.rst index 7597113..fe477bb 100644 --- a/doc/source/user_manual/getting.rst +++ b/doc/source/user_manual/getting.rst @@ -29,15 +29,15 @@ No system administrator privileges are needed to install ctools. .. _sec_installing_gammalib: -Installing gammalib +Installing GammaLib ------------------- -ctools is built on top of gammalib, hence gammalib needs to be built, +ctools is built on top of GammaLib, hence GammaLib needs to be built, installed and configured before ctools can be installed. In case that ctools is already installed on your system you may skip reading this section and continue with the :ref:`sec_building_ctools` section. -You will need the following to build the gammalib: +You will need the following to build the GammaLib: - About 200 MB of free disk space. @@ -58,12 +58,12 @@ You will need the following to build the gammalib: Furthermore it is recommended to have readline and ncurses installed (including the developer packages). -After downloading the gammalib tarball (see :ref:`download`), save it +After downloading the GammaLib tarball (see :ref:`download`), save it in an appropriate location (for example ``$HOME/builds``), and type .. code-block:: bash - $ tar xvfz gammalib-0.11.0.tar.gz + $ tar xvfz gammalib-1.0.0.tar.gz (the ``$`` symbol indicates the console prompt and is not part of the command that you should type in). @@ -73,7 +73,7 @@ typing .. code-block:: bash - $ cd gammalib-0.11.0 + $ cd gammalib-1.0.0 $ ./configure $ make $ make check @@ -82,6 +82,8 @@ typing The last step may require system administator privileges. In this case, use +.. code-block:: bash + $ sudo make install If you do not have system administrator privileges, change the install @@ -133,7 +135,7 @@ an appropriate location (for example ``$HOME/builds``), and type .. code-block:: bash - $ tar xvfz ctools-0.10.0.tar.gz + $ tar xvfz ctools-1.0.0.tar.gz (the ``$`` symbol indicates the console prompt and is not part of the command that you should type in). @@ -142,7 +144,7 @@ Step in the directory and build the ctools by typing .. code-block:: bash - $ cd ctools-0.10.0 + $ cd ctools-1.0.0 $ ./configure $ make @@ -181,7 +183,7 @@ following output at the end of the unit testing: PASS: test_python.py make[4]: Nothing to be done for `all'. ============================================================================ - Testsuite summary for ctools 0.10.0 + Testsuite summary for ctools 1.0.0 ============================================================================ # TOTAL: 3 # PASS: 3 @@ -345,19 +347,8 @@ As a quick check that your setup is okay you can run ``csinfo check``: ===> Your Gammalib / ctools setup is OK. If the setup is not okay, run the ``csinfo info`` command to print -detailed information about your setup. -There's also a ``csinfo list`` command to quickly list the available tools: - -.. code-block:: bash - - $ csinfo - - Print info about Gammalib and ctools to the console. - - Available commands: - csinfo list List available ctools and cscripts - csinfo check Check Gammalib / ctools setup - csinfo info Print Gammalib / ctools setup info +detailed information about your setup. There's also a ``csinfo list`` +command to quickly list the available tools. .. _sec_known_problems: diff --git a/doc/source/user_manual/getting_started/beyond.rst b/doc/source/user_manual/getting_started/beyond.rst index 5019107..92225ab 100644 --- a/doc/source/user_manual/getting_started/beyond.rst +++ b/doc/source/user_manual/getting_started/beyond.rst @@ -3,12 +3,12 @@ Beyond the first steps ---------------------- -Assuming that you have read the :ref:`quickstart ` tutorial, -here now some chapters to read if you intend to do more complex -analyses. -Once you have gone through these chapters you should have a pretty complete +Now that you are familiar with the ctools basics, here are some sections +to read if you intend to do more complex analyses. +Once you have gone through these sections you should have a pretty complete picture about the capabilities of ctools. -And you should be ready for doing your own cutting edge analyses. +And you should be ready for doing your own cutting edge CTA science data +analyses. .. toctree:: :maxdepth: 2 @@ -19,11 +19,3 @@ And you should be ready for doing your own cutting edge analyses. beyond_residual beyond_irf beyond_model - - - - - - - - diff --git a/doc/source/user_manual/getting_started/beyond_combine.rst b/doc/source/user_manual/getting_started/beyond_combine.rst index 487eea8..63cc418 100644 --- a/doc/source/user_manual/getting_started/beyond_combine.rst +++ b/doc/source/user_manual/getting_started/beyond_combine.rst @@ -25,9 +25,9 @@ nebula, each offset by 0.5 deg from the nebula in opposite directions: Lower energy limit (TeV) [0.1] Upper energy limit (TeV) [100.0] Calibration database [prod2] - Instrument response function [South_50h] - Model [$CTOOLS/share/models/crab.xml] - Output event data file or observation definition file [events.fits] events1.fits + 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] events1.fits .. code-block:: bash @@ -40,16 +40,16 @@ nebula, each offset by 0.5 deg from the nebula in opposite directions: Lower energy limit (TeV) [0.1] Upper energy limit (TeV) [100.0] Calibration database [prod2] - Instrument response function [South_50h] - Model [$CTOOLS/share/models/crab.xml] - Output event data file or observation definition file [events1.fits] events2.fits + Instrument response function [South_0.5h] + Input model XML file [$CTOOLS/share/models/crab.xml] + Output event data file or observation definition XML file [events1.fits] events2.fits This will produce the two event files ``events1.fits`` and ``events2.fits`` on disk. -As next step we have to create an **observation definition XML file** that +As next step you have to create an **observation definition XML file** that collects both observations in a list. -Here is how such a file would look like: +Here is how that file looks like: .. code-block:: xml @@ -86,116 +86,140 @@ ctools knows which instrument specific functions need to be called.** Imaging Air Cherenkov Telescope in a joint analysis. **Please recall that instrument strings are case sensitive.** -Now we are ready to do a joint maximum likelihood analysis using -:ref:`ctlike`: +Now you are ready to do a joint maximum likelihood analysis using +:ref:`ctlike`. Type: .. code-block:: bash - $ ctlike - Event list, counts cube or observation definition file [selected_events.fits] obs.xml + $ ctlike chatter=3 + Input event list, counts cube or observation definition XML file [events.fits] obs.xml Calibration database [prod2] - Instrument response function [South_50h] - Source model [$CTOOLS/share/models/crab.xml] - Source model output file [crab_results.xml] + Instrument response function [South_0.5h] + Input model XML file [$CTOOLS/share/models/crab.xml] + Output model XML file [crab_results.xml] -Instead of providing an event list or a counts cube, we now provided the +Instead of providing an event list or a counts cube, you now provided the filename of the observation definition XML file (here ``obs.xml``) as input -parameter. -:ref:`ctlike` recognises this format and automatically performs a -joint analysis, as indicated in the log file: +parameter. :ref:`ctlike` recognises this format and automatically performs a +joint analysis. To see the result of this you needed to specify the +``chatter=3`` parameter so that the chattiness is increased in the log +file: .. code-block:: xml - 2015-05-22T20:17:07: +==============+ - 2015-05-22T20:17:07: | Observations | - 2015-05-22T20:17:07: +==============+ - 2015-05-22T20:17:07: === GObservations === - 2015-05-22T20:17:07: Number of observations ....: 2 - 2015-05-22T20:17:07: Number of predicted events : 0 - 2015-05-22T20:17:07: === GCTAObservation === - 2015-05-22T20:17:07: Name ......................: - 2015-05-22T20:17:07: Identifier ................: 00001 - 2015-05-22T20:17:07: Instrument ................: CTA - 2015-05-22T20:17:07: Event file ................: events1.fits - 2015-05-22T20:17:07: Event type ................: EventList - 2015-05-22T20:17:07: Statistics ................: Poisson - 2015-05-22T20:17:07: Ontime ....................: 1800 s - 2015-05-22T20:17:07: Livetime ..................: 1710 s - 2015-05-22T20:17:07: Deadtime correction .......: 0.95 - 2015-05-22T20:17:07: User energy range .........: undefined - 2015-05-22T20:17:07: === GCTAPointing === - 2015-05-22T20:17:07: Pointing direction ........: (RA,Dec)=(83.63,21.51) + 2015-12-07T21:23:07: +==============+ + 2015-12-07T21:23:07: | Observations | + 2015-12-07T21:23:07: +==============+ + 2015-12-07T21:23:07: === GObservations === + 2015-12-07T21:23:07: Number of observations ....: 2 + 2015-12-07T21:23:07: Number of models ..........: 2 + 2015-12-07T21:23:07: Number of predicted events : 0 + 2015-12-07T21:23:07: === GCTAObservation === + 2015-12-07T21:23:07: Name ......................: Crab + 2015-12-07T21:23:07: Identifier ................: 00001 + 2015-12-07T21:23:07: Instrument ................: CTA + 2015-12-07T21:23:07: Event file ................: events1.fits + 2015-12-07T21:23:07: Event type ................: EventList + 2015-12-07T21:23:07: Statistics ................: Poisson + 2015-12-07T21:23:07: Ontime ....................: 1800 s + 2015-12-07T21:23:07: Livetime ..................: 1710 s + 2015-12-07T21:23:07: Deadtime correction .......: 0.95 + 2015-12-07T21:23:07: User energy range .........: undefined + 2015-12-07T21:23:07: === GCTAPointing === + 2015-12-07T21:23:07: Pointing direction ........: (RA,Dec)=(83.63,21.51) + 2015-12-07T21:23:07: === GCTAResponseIrf === + 2015-12-07T21:23:07: Caldb mission .............: cta + 2015-12-07T21:23:07: Caldb instrument ..........: prod2 + 2015-12-07T21:23:07: Response name .............: South_0.5h + 2015-12-07T21:23:07: Energy dispersion .........: Not used + 2015-12-07T21:23:07: Save energy range .........: undefined + 2015-12-07T21:23:07: === GCTAEventList === + 2015-12-07T21:23:07: Number of events ..........: 22593 + 2015-12-07T21:23:07: Time interval .............: 51544.5 - 51544.5 days + 2015-12-07T21:23:07: Energy interval ...........: 0.1 - 100 TeV + 2015-12-07T21:23:07: Region of interest ........: RA=83.63, DEC=21.51 [0,0] Radius=5 deg + 2015-12-07T21:23:07: === GCTAObservation === + 2015-12-07T21:23:07: Name ......................: Crab + 2015-12-07T21:23:07: Identifier ................: 00002 + 2015-12-07T21:23:07: Instrument ................: CTA + 2015-12-07T21:23:07: Event file ................: events2.fits + 2015-12-07T21:23:07: Event type ................: EventList + 2015-12-07T21:23:07: Statistics ................: Poisson + 2015-12-07T21:23:07: Ontime ....................: 1800 s + 2015-12-07T21:23:07: Livetime ..................: 1710 s + 2015-12-07T21:23:07: Deadtime correction .......: 0.95 + 2015-12-07T21:23:07: User energy range .........: undefined + 2015-12-07T21:23:07: === GCTAPointing === + 2015-12-07T21:23:07: Pointing direction ........: (RA,Dec)=(83.63,22.51) + 2015-12-07T21:23:07: === GCTAResponseIrf === + 2015-12-07T21:23:07: Caldb mission .............: cta + 2015-12-07T21:23:07: Caldb instrument ..........: prod2 + 2015-12-07T21:23:07: Response name .............: South_0.5h + 2015-12-07T21:23:07: Energy dispersion .........: Not used + 2015-12-07T21:23:07: Save energy range .........: undefined + 2015-12-07T21:23:07: === GCTAEventList === + 2015-12-07T21:23:07: Number of events ..........: 22529 + 2015-12-07T21:23:07: Time interval .............: 51544.5 - 51544.5 days + 2015-12-07T21:23:07: Energy interval ...........: 0.1 - 100 TeV + 2015-12-07T21:23:07: Region of interest ........: RA=83.63, DEC=22.51 [0,0] Radius=5 deg ... - 2015-05-22T20:17:07: === GCTAObservation === - 2015-05-22T20:17:07: Name ......................: - 2015-05-22T20:17:07: Identifier ................: 00002 - 2015-05-22T20:17:07: Instrument ................: CTA - 2015-05-22T20:17:07: Event file ................: events2.fits - 2015-05-22T20:17:07: Event type ................: EventList - 2015-05-22T20:17:07: Statistics ................: Poisson - 2015-05-22T20:17:07: Ontime ....................: 1800 s - 2015-05-22T20:17:07: Livetime ..................: 1710 s - 2015-05-22T20:17:07: Deadtime correction .......: 0.95 - 2015-05-22T20:17:07: User energy range .........: undefined - 2015-05-22T20:17:07: === GCTAPointing === - 2015-05-22T20:17:07: Pointing direction ........: (RA,Dec)=(83.63,22.51) + 2015-12-07T21:23:07: +=================================+ + 2015-12-07T21:23:07: | Maximum likelihood optimisation | + 2015-12-07T21:23:07: +=================================+ + 2015-12-07T21:23:07: >Iteration 0: -logL=304176.886, Lambda=1.0e-03 + 2015-12-07T21:23:08: >Iteration 1: -logL=304161.799, Lambda=1.0e-03, delta=15.088, max(|grad|)=12.922804 [Index:3] + 2015-12-07T21:23:08: >Iteration 2: -logL=304161.783, Lambda=1.0e-04, delta=0.015, max(|grad|)=0.098155 [Index:3] + 2015-12-07T21:23:08: >Iteration 3: -logL=304161.783, Lambda=1.0e-05, delta=0.000, max(|grad|)=0.000598 [Index:3] ... - 2015-05-22T20:17:07: +=================================+ - 2015-05-22T20:17:07: | Maximum likelihood optimisation | - 2015-05-22T20:17:07: +=================================+ - 2015-05-22T20:17:07: >Iteration 0: -logL=201836.657, Lambda=1.0e-03 - 2015-05-22T20:17:07: >Iteration 1: -logL=201829.438, Lambda=1.0e-03, delta=7.220, max(|grad|)=22.756507 [Index:7] - 2015-05-22T20:17:07: >Iteration 2: -logL=201829.432, Lambda=1.0e-04, delta=0.005, max(|grad|)=0.086395 [Index:3] - 2015-05-22T20:17:07: >Iteration 3: -logL=201829.432, Lambda=1.0e-05, delta=0.000, max(|grad|)=0.001287 [Index:3] - 2015-05-22T20:17:07: - 2015-05-22T20:17:07: +=========================================+ - 2015-05-22T20:17:07: | Maximum likelihood optimization results | - 2015-05-22T20:17:07: +=========================================+ - 2015-05-22T20:17:07: === GOptimizerLM === - 2015-05-22T20:17:07: Optimized function value ..: 201829.432 - 2015-05-22T20:17:07: Absolute precision ........: 0.005 - 2015-05-22T20:17:07: Acceptable value decrease .: 2 - 2015-05-22T20:17:07: Optimization status .......: converged - 2015-05-22T20:17:07: Number of parameters ......: 10 - 2015-05-22T20:17:07: Number of free parameters .: 4 - 2015-05-22T20:17:07: Number of iterations ......: 3 - 2015-05-22T20:17:07: Lambda ....................: 1e-06 - 2015-05-22T20:17:07: Maximum log likelihood ....: -201829.432 - 2015-05-22T20:17:07: Observed events (Nobs) ...: 31799.000 - 2015-05-22T20:17:07: Predicted events (Npred) ..: 31799.000 (Nobs - Npred = 2.67544e-06) - 2015-05-22T20:17:07: === GModels === - 2015-05-22T20:17:07: Number of models ..........: 2 - 2015-05-22T20:17:07: Number of parameters ......: 10 - 2015-05-22T20:17:07: === GModelSky === - 2015-05-22T20:17:07: Name ......................: Crab - 2015-05-22T20:17:07: Instruments ...............: all - 2015-05-22T20:17:07: Instrument scale factors ..: unity - 2015-05-22T20:17:07: Observation identifiers ...: all - 2015-05-22T20:17:07: Model type ................: PointSource - 2015-05-22T20:17:07: Model components ..........: "SkyDirFunction" * "PowerLaw" * "Constant" - 2015-05-22T20:17:07: Number of parameters ......: 6 - 2015-05-22T20:17:07: Number of spatial par's ...: 2 - 2015-05-22T20:17:07: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1) - 2015-05-22T20:17:07: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1) - 2015-05-22T20:17:07: Number of spectral par's ..: 3 - 2015-05-22T20:17:07: Prefactor ................: 5.8727e-16 +/- 8.04314e-18 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) - 2015-05-22T20:17:07: Index ....................: -2.49825 +/- 0.0118418 [-0,-5] (free,scale=-1,gradient) - 2015-05-22T20:17:07: PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) - 2015-05-22T20:17:07: Number of temporal par's ..: 1 - 2015-05-22T20:17:07: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) - 2015-05-22T20:17:07: === GCTAModelIrfBackground === - 2015-05-22T20:17:07: Name ......................: CTABackgroundModel - 2015-05-22T20:17:07: Instruments ...............: CTA - 2015-05-22T20:17:07: Instrument scale factors ..: unity - 2015-05-22T20:17:07: Observation identifiers ...: all - 2015-05-22T20:17:07: Model type ................: "PowerLaw" * "Constant" - 2015-05-22T20:17:07: Number of parameters ......: 4 - 2015-05-22T20:17:07: Number of spectral par's ..: 3 - 2015-05-22T20:17:07: Prefactor ................: 1.01608 +/- 0.0152751 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient) - 2015-05-22T20:17:07: Index ....................: 0.0226023 +/- 0.0086794 [-5,5] (free,scale=1,gradient) - 2015-05-22T20:17:07: PivotEnergy ..............: 1e+06 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) - 2015-05-22T20:17:07: Number of temporal par's ..: 1 - 2015-05-22T20:17:07: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) + + 2015-12-07T21:23:08: +=========================================+ + 2015-12-07T21:23:08: | Maximum likelihood optimisation results | + 2015-12-07T21:23:08: +=========================================+ + 2015-12-07T21:23:08: === GOptimizerLM === + 2015-12-07T21:23:08: Optimized function value ..: 304161.783 + 2015-12-07T21:23:08: Absolute precision ........: 0.005 + 2015-12-07T21:23:08: Acceptable value decrease .: 2 + 2015-12-07T21:23:08: Optimization status .......: converged + 2015-12-07T21:23:08: Number of parameters ......: 10 + 2015-12-07T21:23:08: Number of free parameters .: 4 + 2015-12-07T21:23:08: Number of iterations ......: 3 + 2015-12-07T21:23:08: Lambda ....................: 1e-06 + 2015-12-07T21:23:08: Maximum log likelihood ....: -304161.783 + 2015-12-07T21:23:08: Observed events (Nobs) ...: 45122.000 + 2015-12-07T21:23:08: Predicted events (Npred) ..: 45122.000 (Nobs - Npred = 1.04131e-06) + 2015-12-07T21:23:08: === GModels === + 2015-12-07T21:23:08: Number of models ..........: 2 + 2015-12-07T21:23:08: Number of parameters ......: 10 + 2015-12-07T21:23:08: === GModelSky === + 2015-12-07T21:23:08: Name ......................: Crab + 2015-12-07T21:23:08: Instruments ...............: all + 2015-12-07T21:23:08: Instrument scale factors ..: unity + 2015-12-07T21:23:08: Observation identifiers ...: all + 2015-12-07T21:23:08: Model type ................: PointSource + 2015-12-07T21:23:08: Model components ..........: "SkyDirFunction" * "PowerLaw" * "Constant" + 2015-12-07T21:23:08: Number of parameters ......: 6 + 2015-12-07T21:23:08: Number of spatial par's ...: 2 + 2015-12-07T21:23:08: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1) + 2015-12-07T21:23:08: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1) + 2015-12-07T21:23:08: Number of spectral par's ..: 3 + 2015-12-07T21:23:08: Prefactor ................: 5.78382e-16 +/- 7.27218e-18 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) + 2015-12-07T21:23:08: Index ....................: -2.5314 +/- 0.0112546 [-0,-5] (free,scale=-1,gradient) + 2015-12-07T21:23:08: PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) + 2015-12-07T21:23:08: Number of temporal par's ..: 1 + 2015-12-07T21:23:08: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) + 2015-12-07T21:23:08: === GCTAModelIrfBackground === + 2015-12-07T21:23:08: Name ......................: CTABackgroundModel + 2015-12-07T21:23:08: Instruments ...............: CTA + 2015-12-07T21:23:08: Instrument scale factors ..: unity + 2015-12-07T21:23:08: Observation identifiers ...: all + 2015-12-07T21:23:08: Model type ................: "PowerLaw" * "Constant" + 2015-12-07T21:23:08: Number of parameters ......: 4 + 2015-12-07T21:23:08: Number of spectral par's ..: 3 + 2015-12-07T21:23:08: Prefactor ................: 1.0008 +/- 0.00839192 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient) + 2015-12-07T21:23:08: Index ....................: 0.0102998 +/- 0.0051965 [-5,5] (free,scale=1,gradient) + 2015-12-07T21:23:08: PivotEnergy ..............: 1e+06 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) + 2015-12-07T21:23:08: Number of temporal par's ..: 1 + 2015-12-07T21:23:08: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) The log file indicates that the fit converged quickly, the spectral parameters of the Crab nebula have now been constrained using the events diff --git a/doc/source/user_manual/getting_started/beyond_irf.rst b/doc/source/user_manual/getting_started/beyond_irf.rst index 2c9d9fa..2275168 100644 --- a/doc/source/user_manual/getting_started/beyond_irf.rst +++ b/doc/source/user_manual/getting_started/beyond_irf.rst @@ -6,9 +6,9 @@ Connecting observations to dedicated instrument response functions So far the same instrument response function has been used for all observations in our analysis examples. You may have recognised that we always specified the calibration database -``prod2`` and the instrument response function ``South_50h`` that +``prod2`` and the instrument response function ``South_0.5h`` that are shipped with ctools. -When combining multiple observations in a joint analysis we may however +When combining multiple observations in a joint analysis you may however need to specify dedicated instrument response functions for the different observations. @@ -26,40 +26,40 @@ the northern array and the other with the southern array: Lower energy limit (TeV) [0.1] Upper energy limit (TeV) [100.0] Calibration database [prod2] - Instrument response function [South_50h] North_50h - Model [$CTOOLS/share/models/crab.xml] - Output event data file or observation definition file [events.fits] events1.fits + Instrument response function [South_0.5h] North_0.5h + Input model XML file [$CTOOLS/share/models/crab.xml] + Output event data file or observation definition XML file [events.fits] north.fits .. code-block:: bash $ ctobssim RA of pointing (degrees) (0-360) [83.63] - Dec of pointing (degrees) (-90-90) [22.01] 22.51 + Dec of pointing (degrees) (-90-90) [21.51] 22.51 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 [North_50h] South_50h - Model [$CTOOLS/share/models/crab.xml] - Output event data file or observation definition file [events1.fits] events2.fits + Instrument response function [North_0.5h] South_0.5h + Input model XML file [$CTOOLS/share/models/crab.xml] + Output event data file or observation definition XML file [north.fits] south.fits -Now we combine the information about these two observations in an +You now need to combine the information about these two observations in an observation definition XML file: .. code-block:: xml - - - - - - - - - - + + + + + + + + + + Each observation now has a ``Calibration`` parameter with the attributes @@ -97,126 +97,136 @@ Within this database there are six response functions: ``North_0.5h``, ``North_5h``, ``North_50h``, ``South_0.5h``, ``South_5h``, and ``South_50h``. -We now jointly fit both observations using :ref:`ctlike`: +You now jointly fit both observations using :ref:`ctlike`: .. code-block:: bash - $ ctlike - Event list, counts cube or observation definition file [obs.xml] obs_irf.xml - Source model [$CTOOLS/share/models/crab.xml] - Source model output file [crab_results.xml] + $ ctlike chatter=3 + Input event list, counts cube or observation definition XML file [events.fits] obs_irf.xml + Input model XML file [$CTOOLS/share/models/crab.xml] + Output model XML file [crab_results.xml] To see the usage of the individual response functions you may inspect the -log file. You will notice that each observation now has a specific response +log file (you need to set ``chatter=3`` to see the details of the +observations in the log file). +You will notice that each observation now has a specific response function and that the filename of the response information differs for both observations. .. code-block:: xml - 2015-05-22T22:06:43: +==============+ - 2015-05-22T22:06:43: | Observations | - 2015-05-22T22:06:43: +==============+ - 2015-05-22T22:06:43: === GObservations === - 2015-05-22T22:06:43: Number of observations ....: 2 - 2015-05-22T22:06:43: Number of predicted events : 0 - 2015-05-22T22:06:43: === GCTAObservation === - 2015-05-22T22:06:43: Name ......................: - 2015-05-22T22:06:43: Identifier ................: 00001 - 2015-05-22T22:06:43: Instrument ................: CTA - 2015-05-22T22:06:43: Event file ................: events1.fits - 2015-05-22T22:06:43: Event type ................: EventList - 2015-05-22T22:06:43: Statistics ................: Poisson - 2015-05-22T22:06:43: Ontime ....................: 1800 s - 2015-05-22T22:06:43: Livetime ..................: 1710 s - 2015-05-22T22:06:43: Deadtime correction .......: 0.95 - 2015-05-22T22:06:43: User energy range .........: undefined - 2015-05-22T22:06:43: === GCTAPointing === - 2015-05-22T22:06:43: Pointing direction ........: (RA,Dec)=(83.63,21.51) - 2015-05-22T22:06:43: === GCTAResponseIrf === - 2015-05-22T22:06:43: Caldb mission .............: cta - 2015-05-22T22:06:43: Caldb instrument ..........: prod2 - 2015-05-22T22:06:43: Response name .............: North_50h - 2015-05-22T22:06:43: Energy dispersion .........: Not used - 2015-05-22T22:06:43: Save energy range .........: undefined + 2015-12-07T22:22:18: +==============+ + 2015-12-07T22:22:18: | Observations | + 2015-12-07T22:22:18: +==============+ + 2015-12-07T22:22:18: === GObservations === + 2015-12-07T22:22:18: Number of observations ....: 2 + 2015-12-07T22:22:18: Number of models ..........: 2 + 2015-12-07T22:22:18: Number of predicted events : 0 + 2015-12-07T22:22:18: === GCTAObservation === + 2015-12-07T22:22:18: Name ......................: Crab + 2015-12-07T22:22:18: Identifier ................: 00001 + 2015-12-07T22:22:18: Instrument ................: CTA + 2015-12-07T22:22:18: Event file ................: north.fits + 2015-12-07T22:22:18: Event type ................: EventList + 2015-12-07T22:22:18: Statistics ................: Poisson + 2015-12-07T22:22:18: Ontime ....................: 1800 s + 2015-12-07T22:22:18: Livetime ..................: 1710 s + 2015-12-07T22:22:18: Deadtime correction .......: 0.95 + 2015-12-07T22:22:18: User energy range .........: undefined + 2015-12-07T22:22:18: === GCTAPointing === + 2015-12-07T22:22:18: Pointing direction ........: (RA,Dec)=(83.63,21.51) + 2015-12-07T22:22:18: === GCTAResponseIrf === + 2015-12-07T22:22:18: Caldb mission .............: cta + 2015-12-07T22:22:18: Caldb instrument ..........: prod2 + 2015-12-07T22:22:18: Response name .............: North_0.5h + 2015-12-07T22:22:18: Energy dispersion .........: Not used + 2015-12-07T22:22:18: Save energy range .........: undefined + 2015-12-07T22:22:18: === GCTAEventList === + 2015-12-07T22:22:18: Number of events ..........: 8469 + 2015-12-07T22:22:18: Time interval .............: 51544.5 - 51544.5 days + 2015-12-07T22:22:18: Energy interval ...........: 0.1 - 100 TeV + 2015-12-07T22:22:18: Region of interest ........: RA=83.63, DEC=21.51 [0,0] Radius=5 deg + 2015-12-07T22:22:18: === GCTAObservation === + 2015-12-07T22:22:18: Name ......................: Crab + 2015-12-07T22:22:18: Identifier ................: 00002 + 2015-12-07T22:22:18: Instrument ................: CTA + 2015-12-07T22:22:18: Event file ................: south.fits + 2015-12-07T22:22:18: Event type ................: EventList + 2015-12-07T22:22:18: Statistics ................: Poisson + 2015-12-07T22:22:18: Ontime ....................: 1800 s + 2015-12-07T22:22:18: Livetime ..................: 1710 s + 2015-12-07T22:22:18: Deadtime correction .......: 0.95 + 2015-12-07T22:22:18: User energy range .........: undefined + 2015-12-07T22:22:18: === GCTAPointing === + 2015-12-07T22:22:18: Pointing direction ........: (RA,Dec)=(83.63,22.51) + 2015-12-07T22:22:18: === GCTAResponseIrf === + 2015-12-07T22:22:18: Caldb mission .............: cta + 2015-12-07T22:22:18: Caldb instrument ..........: prod2 + 2015-12-07T22:22:18: Response name .............: South_0.5h + 2015-12-07T22:22:18: Energy dispersion .........: Not used + 2015-12-07T22:22:18: Save energy range .........: undefined + 2015-12-07T22:22:18: === GCTAEventList === + 2015-12-07T22:22:18: Number of events ..........: 22529 + 2015-12-07T22:22:18: Time interval .............: 51544.5 - 51544.5 days + 2015-12-07T22:22:18: Energy interval ...........: 0.1 - 100 TeV + 2015-12-07T22:22:18: Region of interest ........: RA=83.63, DEC=22.51 [0,0] Radius=5 deg ... - 2015-05-22T22:06:43: === GCTAObservation === - 2015-05-22T22:06:43: Name ......................: - 2015-05-22T22:06:43: Identifier ................: 00002 - 2015-05-22T22:06:43: Instrument ................: CTA - 2015-05-22T22:06:43: Event file ................: events2.fits - 2015-05-22T22:06:43: Event type ................: EventList - 2015-05-22T22:06:43: Statistics ................: Poisson - 2015-05-22T22:06:43: Ontime ....................: 1800 s - 2015-05-22T22:06:43: Livetime ..................: 1710 s - 2015-05-22T22:06:43: Deadtime correction .......: 0.95 - 2015-05-22T22:06:43: User energy range .........: undefined - 2015-05-22T22:06:43: === GCTAPointing === - 2015-05-22T22:06:43: Pointing direction ........: (RA,Dec)=(83.63,22.51) - 2015-05-22T22:06:43: === GCTAResponseIrf === - 2015-05-22T22:06:43: Caldb mission .............: cta - 2015-05-22T22:06:43: Caldb instrument ..........: prod2 - 2015-05-22T22:06:43: Response name .............: South_50h - 2015-05-22T22:06:43: Energy dispersion .........: Not used - 2015-05-22T22:06:43: Save energy range .........: undefined + 2015-12-07T22:22:18: +=================================+ + 2015-12-07T22:22:18: | Maximum likelihood optimisation | + 2015-12-07T22:22:18: +=================================+ + 2015-12-07T22:22:18: >Iteration 0: -logL=214691.468, Lambda=1.0e-03 + 2015-12-07T22:22:18: >Iteration 1: -logL=214679.380, Lambda=1.0e-03, delta=12.088, max(|grad|)=9.713100 [Index:3] + 2015-12-07T22:22:18: >Iteration 2: -logL=214679.365, Lambda=1.0e-04, delta=0.015, max(|grad|)=0.092975 [Index:3] + 2015-12-07T22:22:18: >Iteration 3: -logL=214679.365, Lambda=1.0e-05, delta=0.000, max(|grad|)=0.000749 [Index:3] ... - 2015-05-22T22:06:43: +=================================+ - 2015-05-22T22:06:43: | Maximum likelihood optimisation | - 2015-05-22T22:06:43: +=================================+ - 2015-05-22T22:06:43: >Iteration 0: -logL=138039.472, Lambda=1.0e-03 - 2015-05-22T22:06:43: >Iteration 1: -logL=138035.501, Lambda=1.0e-03, delta=3.971, max(|grad|)=16.203761 [Index:7] - 2015-05-22T22:06:43: >Iteration 2: -logL=138035.496, Lambda=1.0e-04, delta=0.005, max(|grad|)=-0.053284 [Index:7] - 2015-05-22T22:06:43: >Iteration 3: -logL=138035.496, Lambda=1.0e-05, delta=0.000, max(|grad|)=0.001386 [Index:7] - 2015-05-22T22:06:44: - 2015-05-22T22:06:44: +=========================================+ - 2015-05-22T22:06:44: | Maximum likelihood optimization results | - 2015-05-22T22:06:44: +=========================================+ - 2015-05-22T22:06:44: === GOptimizerLM === - 2015-05-22T22:06:44: Optimized function value ..: 138035.496 - 2015-05-22T22:06:44: Absolute precision ........: 0.005 - 2015-05-22T22:06:44: Acceptable value decrease .: 2 - 2015-05-22T22:06:44: Optimization status .......: converged - 2015-05-22T22:06:44: Number of parameters ......: 10 - 2015-05-22T22:06:44: Number of free parameters .: 4 - 2015-05-22T22:06:44: Number of iterations ......: 3 - 2015-05-22T22:06:44: Lambda ....................: 1e-06 - 2015-05-22T22:06:44: Maximum log likelihood ....: -138035.496 - 2015-05-22T22:06:44: Observed events (Nobs) ...: 20814.000 - 2015-05-22T22:06:44: Predicted events (Npred) ..: 20814.000 (Nobs - Npred = 2.67698e-06) - 2015-05-22T22:06:44: === GModels === - 2015-05-22T22:06:44: Number of models ..........: 2 - 2015-05-22T22:06:44: Number of parameters ......: 10 - 2015-05-22T22:06:44: === GModelSky === - 2015-05-22T22:06:44: Name ......................: Crab - 2015-05-22T22:06:44: Instruments ...............: all - 2015-05-22T22:06:44: Instrument scale factors ..: unity - 2015-05-22T22:06:44: Observation identifiers ...: all - 2015-05-22T22:06:44: Model type ................: PointSource - 2015-05-22T22:06:44: Model components ..........: "SkyDirFunction" * "PowerLaw" * "Constant" - 2015-05-22T22:06:44: Number of parameters ......: 6 - 2015-05-22T22:06:44: Number of spatial par's ...: 2 - 2015-05-22T22:06:44: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1) - 2015-05-22T22:06:44: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1) - 2015-05-22T22:06:44: Number of spectral par's ..: 3 - 2015-05-22T22:06:44: Prefactor ................: 5.9184e-16 +/- 9.68924e-18 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) - 2015-05-22T22:06:44: Index ....................: -2.49569 +/- 0.0145376 [-0,-5] (free,scale=-1,gradient) - 2015-05-22T22:06:44: PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) - 2015-05-22T22:06:44: Number of temporal par's ..: 1 - 2015-05-22T22:06:44: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) - 2015-05-22T22:06:44: === GCTAModelIrfBackground === - 2015-05-22T22:06:44: Name ......................: CTABackgroundModel - 2015-05-22T22:06:44: Instruments ...............: CTA - 2015-05-22T22:06:44: Instrument scale factors ..: unity - 2015-05-22T22:06:44: Observation identifiers ...: all - 2015-05-22T22:06:44: Model type ................: "PowerLaw" * "Constant" - 2015-05-22T22:06:44: Number of parameters ......: 4 - 2015-05-22T22:06:44: Number of spectral par's ..: 3 - 2015-05-22T22:06:44: Prefactor ................: 1.02524 +/- 0.0182054 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient) - 2015-05-22T22:06:44: Index ....................: 0.0229728 +/- 0.0104713 [-5,5] (free,scale=1,gradient) - 2015-05-22T22:06:44: PivotEnergy ..............: 1e+06 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) - 2015-05-22T22:06:44: Number of temporal par's ..: 1 - 2015-05-22T22:06:44: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) - - + 2015-12-07T22:22:18: +=========================================+ + 2015-12-07T22:22:18: | Maximum likelihood optimisation results | + 2015-12-07T22:22:18: +=========================================+ + 2015-12-07T22:22:18: === GOptimizerLM === + 2015-12-07T22:22:18: Optimized function value ..: 214679.365 + 2015-12-07T22:22:18: Absolute precision ........: 0.005 + 2015-12-07T22:22:18: Acceptable value decrease .: 2 + 2015-12-07T22:22:18: Optimization status .......: converged + 2015-12-07T22:22:18: Number of parameters ......: 10 + 2015-12-07T22:22:18: Number of free parameters .: 4 + 2015-12-07T22:22:18: Number of iterations ......: 3 + 2015-12-07T22:22:18: Lambda ....................: 1e-06 + 2015-12-07T22:22:18: Maximum log likelihood ....: -214679.365 + 2015-12-07T22:22:18: Observed events (Nobs) ...: 30998.000 + 2015-12-07T22:22:18: Predicted events (Npred) ..: 30998.000 (Nobs - Npred = 1.02755e-06) + 2015-12-07T22:22:18: === GModels === + 2015-12-07T22:22:18: Number of models ..........: 2 + 2015-12-07T22:22:18: Number of parameters ......: 10 + 2015-12-07T22:22:18: === GModelSky === + 2015-12-07T22:22:18: Name ......................: Crab + 2015-12-07T22:22:18: Instruments ...............: all + 2015-12-07T22:22:18: Instrument scale factors ..: unity + 2015-12-07T22:22:18: Observation identifiers ...: all + 2015-12-07T22:22:18: Model type ................: PointSource + 2015-12-07T22:22:18: Model components ..........: "SkyDirFunction" * "PowerLaw" * "Constant" + 2015-12-07T22:22:18: Number of parameters ......: 6 + 2015-12-07T22:22:18: Number of spatial par's ...: 2 + 2015-12-07T22:22:18: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1) + 2015-12-07T22:22:18: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1) + 2015-12-07T22:22:18: Number of spectral par's ..: 3 + 2015-12-07T22:22:18: Prefactor ................: 5.81237e-16 +/- 8.63534e-18 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) + 2015-12-07T22:22:18: Index ....................: -2.53954 +/- 0.0137432 [-0,-5] (free,scale=-1,gradient) + 2015-12-07T22:22:18: PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) + 2015-12-07T22:22:18: Number of temporal par's ..: 1 + 2015-12-07T22:22:18: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) + 2015-12-07T22:22:18: === GCTAModelIrfBackground === + 2015-12-07T22:22:18: Name ......................: CTABackgroundModel + 2015-12-07T22:22:18: Instruments ...............: CTA + 2015-12-07T22:22:18: Instrument scale factors ..: unity + 2015-12-07T22:22:18: Observation identifiers ...: all + 2015-12-07T22:22:18: Model type ................: "PowerLaw" * "Constant" + 2015-12-07T22:22:18: Number of parameters ......: 4 + 2015-12-07T22:22:18: Number of spectral par's ..: 3 + 2015-12-07T22:22:18: Prefactor ................: 1.01067 +/- 0.0103997 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient) + 2015-12-07T22:22:18: Index ....................: 0.0135414 +/- 0.0064452 [-5,5] (free,scale=1,gradient) + 2015-12-07T22:22:18: PivotEnergy ..............: 1e+06 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) + 2015-12-07T22:22:18: Number of temporal par's ..: 1 + 2015-12-07T22:22:18: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) A more fine grained control over the response function can be achieved by specifying individual filenames for the various response components. @@ -228,18 +238,18 @@ This is definitely expert mode, to be used with utmost care. - - - - - + + + + + - - - - - + + + + + @@ -264,5 +274,3 @@ file is given below: - - diff --git a/doc/source/user_manual/getting_started/beyond_model.rst b/doc/source/user_manual/getting_started/beyond_model.rst index 0621daa..2db8335 100644 --- a/doc/source/user_manual/getting_started/beyond_model.rst +++ b/doc/source/user_manual/getting_started/beyond_model.rst @@ -26,20 +26,15 @@ specified by the following observation definition XML file: -Each observation has a **unique identifier string** and we can use this -unique string to connect that observation to a specific model component. +Each observation has a **unique identifier string** and you can use this +unique string to connect an observation to a specific model component. This is illustrated in the model XML file below that now has two background components. -**Model components in a XML file need to have a unique name**, for this -reason we have named the components ``Background_00001`` and -``Background_00002``. -Both components now have an ``id`` attribute that ties the model to -an observation. -In other words, the model is only applicable to a specific observation -(it will be ignored for all other observations). -In our case, model ``Background_00001`` applies to the ``events1.fits`` -event list while model ``Background_00002`` applies to the ``events2.fits`` -event list. +Both components have an ``id`` attribute that ties the model to a specific +observation. +Using ``id="00001"``, the first background model is tied to the +``events1.fits`` file, using ``id="00002"``, the second background model is +tied to the ``events2.fits`` file. .. code-block:: xml @@ -72,13 +67,21 @@ event list. -Note that both background model components distinguish in their event rate; -component ``Background_00002`` has a ten times larger background rate than -component ``Background_00001``. -To illustrate this difference, we simulated two observations offset by -+/- 2 degrees from the Crab nebulae using :ref:`ctobssim`. -We then performed a joint fit of the observations using :ref:`ctlike` -and created a resdiual map using :ref:`csresmap` that is shown below. +.. note:: + + **Model components in a XML file need to have a unique name.** + For this reason the components are named here ``Background_00001`` + and ``Background_00002``. + +Note that both background model components distinguish in their event rate. +``Background_00002`` has a ten times larger background rate than +``Background_00001``. +To illustrate this difference, two observations offset by +/- 2 degrees from +the Crab nebulae have been simulated using :ref:`ctobssim`. +A joint fit of the observations is then performed using :ref:`ctlike` and +a residual map is created using :ref:`csresmap`. + +The result is shown below. Obviously, the upper observation has larger residuals owing to the ten times larger background rate that was used in the simulation of the associated data. diff --git a/doc/source/user_manual/getting_started/beyond_residual.rst b/doc/source/user_manual/getting_started/beyond_residual.rst index 9eb2924..aa8749e 100644 --- a/doc/source/user_manual/getting_started/beyond_residual.rst +++ b/doc/source/user_manual/getting_started/beyond_residual.rst @@ -22,11 +22,11 @@ Suppose that you have done the following unbinned ``ctlike`` analysis: .. code-block:: bash $ ctlike - Event list, counts cube or observation definition file [events.fits] obs.xml + Input event list, counts cube or observation definition XML file [events.fits] obs.xml Calibration database [prod2] - Instrument response function [South_50h] - Source model [$CTOOLS/share/models/crab.xml] - Source model output file [crab_results.xml] + Instrument response function [South_0.5h] + Input model XML file [$CTOOLS/share/models/crab.xml] + Output model XML file [crab_results.xml] The :ref:`csresmap` script is then used as follows: @@ -34,18 +34,19 @@ The :ref:`csresmap` script is then used as follows: $ csresmap Parfile csresmap.par not found. Create default parfile. - Event list, counts cube, or observation definition file [events.fits] obs.xml + Input event list, counts cube, or observation definition XML file [events.fits] obs.xml Calibration database [prod2] - Instrument response function [South_50h] - Source model [$CTOOLS/share/models/crab.xml] crab_results.xml - Output residual map [resmap.fits] - First coordinate of image center in degrees (RA or galactic l) [83.63] - Second coordinate of image center in degrees (DEC or galactic b) [22.01] + Instrument response function [South_0.5h] + Input model XML file [$CTOOLS/share/models/crab.xml] crab_results.xml + 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] Coordinate System (CEL|GAL) [CEL] - Projection method e.g. AIT|AZP|CAR|MER|MOL|STG|TAN (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR] + Projection method (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR] Size of the X axis in pixels [200] Size of the Y axis in pixels [200] Pixel size (deg/pixel) [0.02] + Output residual map file [resmap.fits] + Residual map computation algorithm (SUB|SUBDIV|SUBDIVSQRT) [SUBDIV] In our example, the map will cover a field of 4 x 4 degrees centred on the Crab nebula at a pixel scale of 0.02 degrees per pixel. diff --git a/doc/source/user_manual/getting_started/beyond_stacked.rst b/doc/source/user_manual/getting_started/beyond_stacked.rst index 4e7b9a4..537d443 100644 --- a/doc/source/user_manual/getting_started/beyond_stacked.rst +++ b/doc/source/user_manual/getting_started/beyond_stacked.rst @@ -6,16 +6,15 @@ Performing a stacked analysis A stacked analysis is a binned analysis where all data from multiple observations are stacked into a single counts cube. The event stacking is done using the :ref:`ctbin` tool. -Instead of providing to :ref:`ctbin` an event list we -now specify the observation definition XML file ``obs.xml`` -on input. +Instead of providing to :ref:`ctbin` an event list you should now +specify the observation definition XML file ``obs.xml`` on input. :ref:`ctbin` will then loop over all observations and collect all events into a single counts cube: .. code-block:: bash $ ctbin - Event list or observation definition file [events.fits] obs.xml + Input event list or observation definition XML file [events.fits] obs.xml 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] @@ -27,37 +26,37 @@ into a single counts cube: 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 [cntcube.fits] + Output counts cube file [cntcube.fits] -We now have a stacked counts cube ``cntcube.fits`` on disk. -Before we can use that counts cube in a maximum likelihood -analysis, we have to compute the instrument response and the +You now have a stacked counts cube ``cntcube.fits`` on disk. +Before you can use that counts cube in a maximum likelihood +analysis, you have to compute the instrument response and the background model that are needed to describe the stacked data. -For the former, we have to compute the total exposure for the stacked +For the former, you have to compute the total exposure for the stacked cube (i.e. the sum of the effective areas for each observation multiplied by the corresponding livetimes) and an effective point spread function (i.e. the point spread function of the different observations weighted by the corresponding exposures). -To get both informations we use the :ref:`ctexpcube` and +To get both informations you use the :ref:`ctexpcube` and :ref:`ctpsfcube` tools: .. code-block:: bash $ ctexpcube - Event list or observation definition file [NONE] obs.xml + Input event list or observation definition XML file [NONE] obs.xml Calibration database [prod2] - Instrument response function [South_50h] - Counts cube for exposure cube definition [NONE] cntcube.fits + Instrument response function [South_0.5h] + Input counts cube file to extract exposure cube definition [NONE] cntcube.fits Output exposure cube file [expcube.fits] .. code-block:: bash $ ctpsfcube - Event list or observation definition file [NONE] obs.xml + Input event list or observation definition XML file [NONE] obs.xml Calibration database [prod2] - Instrument response function [South_50h] - Counts cube for psf cube definition [NONE] + 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|MER|MOL|STG|TAN) [CAR] @@ -68,42 +67,41 @@ To get both informations we use the :ref:`ctexpcube` and Lower energy limit (TeV) [0.1] Upper energy limit (TeV) [100.0] Number of energy bins [20] - Output psf cube file [psfcube.fits] - -We provide the ``obs.xml`` file on input to inform both tools which -observations have been combined. -For :ref:`ctexpcube` we further provide the counts cube so that the -tool can copy the exposure cube definition (number of spatial pixels -and pixel size, number of energy bins) from the counts cube. -This minimises the number of further user parameters that need to be -provided and assures an exposure cube that is compatible with the counts -cube. -For :ref:`ctpsfcube` we do not use the counts cube for the PSF cube -definition as this would lead to a hugh file owing to the fine spatial + Output PSF cube file [psfcube.fits] + +The ``obs.xml`` file has been provided on input to specify for both tools +which observations have been combined. +You further provided the counts cube on input so that :ref:`ctexpcube` can +read the cub definition from that file and apply it to the exposure cube. +This is a trick to reduce the number of user parameters that you need to +provide. + +You do not apply this trick when using :ref:`ctpsfcube` as this +would lead to a hugh output file owing to the fine spatial pixelisation of the counts cube. -Since the PSF evolves only slowly over the field of fiew, we provide a -rather coarse spatial binning of 1 degree covering a grid of 10 x 10 -degrees around the centre of the counts cube. -For the energy binning, we use the same logarithmic binning that has -also been used for the counts cube. +Such a fine binning is not needed for the PSF cube, as the PSF evolves +only slowly over the field of view. +It is thus sufficient to compute a PSF cube with a rather coarse spatial +binning; here you used a spatial binning of 1 degree covering a grid of +10 x 10 degrees. -As final step of the analysis preparation, we need to generate a +As final step of the analysis preparation, you need to generate a background cube using the :ref:`ctbkgcube` tool: .. code-block:: bash $ ctbkgcube - Input event list or observation definition file [NONE] obs.xml + Input event list or observation definition XML file [NONE] obs.xml Calibration database [prod2] - Instrument response function [South_50h] - Counts cube for background cube definition [NONE] cntcube.fits + Instrument response function [South_0.5h] + Input counts cube file to extract background cube definition [NONE] cntcube.fits Input model XML file [NONE] $CTOOLS/share/models/crab.xml Output background cube file [bkgcube.fits] Output model XML file [NONE] model.xml The usage of :ref:`ctbkgcube` is very similar to that of :ref:`ctexpcube`, yet it takes the model XML file as an additional input parameter. -We here use the usual ``$CTOOLS/share/models/crab.xml`` model +You used here the usual ``$CTOOLS/share/models/crab.xml`` model file that is shipped with the ctools. :ref:`ctbkgcube` provides on output the background cube file ``bkgcube.fits`` and the model XML file ``model.xml`` that can @@ -115,7 +113,7 @@ modelling works: - + @@ -157,88 +155,88 @@ as an appropriate multiplicator to the background cube. Using the same binning for the exposure cube, the background cube and the counts cube is only a convenience. -Now we have all files at hand to launch a stacked maximum likelihood +Now you have all files at hand to perform a stacked maximum likelihood analysis using the :ref:`ctlike` tool: .. code-block:: bash $ ctlike - Event list, counts cube or observation definition file [obs.xml] cntcube.fits - Exposure cube file (only needed for stacked analysis) [NONE] expcube.fits - PSF cube file (only needed for stacked analysis) [NONE] psfcube.fits - Background cube file (only needed for stacked analysis) [NONE] bkgcube.fits - Source model [$CTOOLS/share/models/crab.xml] model.xml - Source model output file [crab_results.xml] + Input event list, counts cube or observation definition XML file [obs.xml] 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 XML file [$CTOOLS/share/models/crab.xml] model.xml + Output model XML file [crab_results.xml] :ref:`ctlike` recognises that a counts cube should be analysed and queries for the exposure cube, the PSF cube, and the background cube file names. -We specified the names of the files produced by the :ref:`ctexpcube`, +You specified here the names of the files produced by the :ref:`ctexpcube`, the :ref:`ctpsfcube` and the :ref:`ctbkgcube` tools. -Furthermore we provided as model the ``model.xml`` file that has been +Furthermore you provided as model the ``model.xml`` file that has been generated by the :ref:`ctbkgcube` tool. The log file of the :ref:`ctlike` run is shown below. Note that the spectral model that is multiplied with the background -cube has a Prefactor of 1.09 +/- 0.02 and an Index of 0.03 +/- 0.01, +cube has a Prefactor of 0.998 +/- 0.011 and an Index of 0.006 +/- 0.007, indicating a very small correction to the actual spectrum of the background cube. Real life situations may of course require larger correction factors. .. code-block:: xml - 2015-05-22T21:10:09: +=================================+ - 2015-05-22T21:10:09: | Maximum likelihood optimisation | - 2015-05-22T21:10:09: +=================================+ - 2015-05-22T21:10:10: >Iteration 0: -logL=60663.798, Lambda=1.0e-03 - 2015-05-22T21:10:11: >Iteration 1: -logL=60642.413, Lambda=1.0e-03, delta=21.385, max(|grad|)=120.031102 [Index:7] - 2015-05-22T21:10:13: >Iteration 2: -logL=60642.233, Lambda=1.0e-04, delta=0.180, max(|grad|)=1.139187 [Index:7] - 2015-05-22T21:10:14: >Iteration 3: -logL=60642.233, Lambda=1.0e-05, delta=0.000, max(|grad|)=-0.001247 [Index:7] - 2015-05-22T21:10:16: - 2015-05-22T21:10:16: +=========================================+ - 2015-05-22T21:10:16: | Maximum likelihood optimization results | - 2015-05-22T21:10:16: +=========================================+ - 2015-05-22T21:10:16: === GOptimizerLM === - 2015-05-22T21:10:16: Optimized function value ..: 60642.233 - 2015-05-22T21:10:16: Absolute precision ........: 0.005 - 2015-05-22T21:10:16: Acceptable value decrease .: 2 - 2015-05-22T21:10:16: Optimization status .......: converged - 2015-05-22T21:10:16: Number of parameters ......: 10 - 2015-05-22T21:10:16: Number of free parameters .: 4 - 2015-05-22T21:10:16: Number of iterations ......: 3 - 2015-05-22T21:10:16: Lambda ....................: 1e-06 - 2015-05-22T21:10:16: Maximum log likelihood ....: -60642.233 - 2015-05-22T21:10:16: Observed events (Nobs) ...: 25963.000 - 2015-05-22T21:10:16: Predicted events (Npred) ..: 25963.000 (Nobs - Npred = 1.5749e-05) - 2015-05-22T21:10:16: === GModels === - 2015-05-22T21:10:16: Number of models ..........: 2 - 2015-05-22T21:10:16: Number of parameters ......: 10 - 2015-05-22T21:10:16: === GModelSky === - 2015-05-22T21:10:16: Name ......................: Crab - 2015-05-22T21:10:16: Instruments ...............: all - 2015-05-22T21:10:16: Instrument scale factors ..: unity - 2015-05-22T21:10:16: Observation identifiers ...: all - 2015-05-22T21:10:16: Model type ................: PointSource - 2015-05-22T21:10:16: Model components ..........: "SkyDirFunction" * "PowerLaw" * "Constant" - 2015-05-22T21:10:16: Number of parameters ......: 6 - 2015-05-22T21:10:16: Number of spatial par's ...: 2 - 2015-05-22T21:10:16: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1) - 2015-05-22T21:10:16: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1) - 2015-05-22T21:10:16: Number of spectral par's ..: 3 - 2015-05-22T21:10:16: Prefactor ................: 5.8499e-16 +/- 8.02289e-18 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) - 2015-05-22T21:10:16: Index ....................: -2.49909 +/- 0.0118976 [-0,-5] (free,scale=-1,gradient) - 2015-05-22T21:10:16: PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) - 2015-05-22T21:10:16: Number of temporal par's ..: 1 - 2015-05-22T21:10:16: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) - 2015-05-22T21:10:16: === GCTAModelCubeBackground === - 2015-05-22T21:10:16: Name ......................: BackgroundModel - 2015-05-22T21:10:16: Instruments ...............: CTA, HESS, MAGIC, VERITAS - 2015-05-22T21:10:16: Instrument scale factors ..: unity - 2015-05-22T21:10:16: Observation identifiers ...: all - 2015-05-22T21:10:16: Model type ................: "PowerLaw" * "Constant" - 2015-05-22T21:10:16: Number of parameters ......: 4 - 2015-05-22T21:10:16: Number of spectral par's ..: 3 - 2015-05-22T21:10:16: Prefactor ................: 1.08781 +/- 0.0212678 [0,infty[ ph/cm2/s/MeV (free,scale=1,gradient) - 2015-05-22T21:10:16: Index ....................: 0.0269681 +/- 0.0111532 [-10,10] (free,scale=1,gradient) - 2015-05-22T21:10:16: PivotEnergy ..............: 1e+06 MeV (fixed,scale=1e+06,gradient) - 2015-05-22T21:10:16: Number of temporal par's ..: 1 - 2015-05-22T21:10:16: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) + 2015-12-07T21:37:45: +=================================+ + 2015-12-07T21:37:45: | Maximum likelihood optimisation | + 2015-12-07T21:37:45: +=================================+ + 2015-12-07T21:37:46: >Iteration 0: -logL=83748.453, Lambda=1.0e-03 + 2015-12-07T21:37:47: >Iteration 1: -logL=83736.039, Lambda=1.0e-03, delta=12.413, max(|grad|)=15.134616 [Index:3] + 2015-12-07T21:37:48: >Iteration 2: -logL=83736.021, Lambda=1.0e-04, delta=0.018, max(|grad|)=0.108613 [Index:3] + 2015-12-07T21:37:49: >Iteration 3: -logL=83736.021, Lambda=1.0e-05, delta=0.000, max(|grad|)=0.000570 [Index:3] + ... + 2015-12-07T21:37:51: +=========================================+ + 2015-12-07T21:37:51: | Maximum likelihood optimisation results | + 2015-12-07T21:37:51: +=========================================+ + 2015-12-07T21:37:51: === GOptimizerLM === + 2015-12-07T21:37:51: Optimized function value ..: 83736.021 + 2015-12-07T21:37:51: Absolute precision ........: 0.005 + 2015-12-07T21:37:51: Acceptable value decrease .: 2 + 2015-12-07T21:37:51: Optimization status .......: converged + 2015-12-07T21:37:51: Number of parameters ......: 10 + 2015-12-07T21:37:51: Number of free parameters .: 4 + 2015-12-07T21:37:51: Number of iterations ......: 3 + 2015-12-07T21:37:51: Lambda ....................: 1e-06 + 2015-12-07T21:37:51: Maximum log likelihood ....: -83736.021 + 2015-12-07T21:37:51: Observed events (Nobs) ...: 35198.000 + 2015-12-07T21:37:51: Predicted events (Npred) ..: 35198.000 (Nobs - Npred = 7.6773e-07) + 2015-12-07T21:37:51: === GModels === + 2015-12-07T21:37:51: Number of models ..........: 2 + 2015-12-07T21:37:51: Number of parameters ......: 10 + 2015-12-07T21:37:51: === GModelSky === + 2015-12-07T21:37:51: Name ......................: Crab + 2015-12-07T21:37:51: Instruments ...............: all + 2015-12-07T21:37:51: Instrument scale factors ..: unity + 2015-12-07T21:37:51: Observation identifiers ...: all + 2015-12-07T21:37:51: Model type ................: PointSource + 2015-12-07T21:37:51: Model components ..........: "SkyDirFunction" * "PowerLaw" * "Constant" + 2015-12-07T21:37:51: Number of parameters ......: 6 + 2015-12-07T21:37:51: Number of spatial par's ...: 2 + 2015-12-07T21:37:51: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1) + 2015-12-07T21:37:51: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1) + 2015-12-07T21:37:51: Number of spectral par's ..: 3 + 2015-12-07T21:37:51: Prefactor ................: 5.75289e-16 +/- 7.24749e-18 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) + 2015-12-07T21:37:51: Index ....................: -2.53122 +/- 0.0113068 [-0,-5] (free,scale=-1,gradient) + 2015-12-07T21:37:51: PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) + 2015-12-07T21:37:51: Number of temporal par's ..: 1 + 2015-12-07T21:37:51: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) + 2015-12-07T21:37:51: === GCTAModelCubeBackground === + 2015-12-07T21:37:51: Name ......................: BackgroundModel + 2015-12-07T21:37:51: Instruments ...............: CTA, HESS, MAGIC, VERITAS + 2015-12-07T21:37:51: Instrument scale factors ..: unity + 2015-12-07T21:37:51: Observation identifiers ...: all + 2015-12-07T21:37:51: Model type ................: "PowerLaw" * "Constant" + 2015-12-07T21:37:51: Number of parameters ......: 4 + 2015-12-07T21:37:51: Number of spectral par's ..: 3 + 2015-12-07T21:37:51: Prefactor ................: 0.998055 +/- 0.0114979 [0.01,100] ph/cm2/s/MeV (free,scale=1,gradient) + 2015-12-07T21:37:51: Index ....................: 0.00648796 +/- 0.00697365 [-5,5] (free,scale=1,gradient) + 2015-12-07T21:37:51: PivotEnergy ..............: 1e+06 MeV (fixed,scale=1e+06,gradient) + 2015-12-07T21:37:51: Number of temporal par's ..: 1 + 2015-12-07T21:37:51: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) diff --git a/doc/source/user_manual/getting_started/beyond_tsmap.rst b/doc/source/user_manual/getting_started/beyond_tsmap.rst index c34be7e..6d5bd17 100644 --- a/doc/source/user_manual/getting_started/beyond_tsmap.rst +++ b/doc/source/user_manual/getting_started/beyond_tsmap.rst @@ -3,13 +3,14 @@ Generating a Test Statistic map ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -In the above examples we were always considering that we knew exactly +In the above examples you have always considering that you knew exactly where the source is located, but in real life situations this may not necessarily be the case. -We thus need a means to search for significant excess emission in our data. +You thus need a means to search for significant excess emission in your data. One possibility to accomplish this task is to compute the Test Statistic for a specific source at a grid of trial positions. The :ref:`cttsmap` tool will exactly do that job. + The Test Statistic is defined as twice the log-likelihood difference between fitting a source at a given position on top of a model or fitting no source. @@ -18,7 +19,7 @@ the source detection significance in Gaussian sigmas, although the exact relation depends somewhat on the formulation of the statistical problem. -Here an example where we apply the :ref:`cttsmap` tool to our two +Here an example of how to apply the :ref:`cttsmap` tool to the two simulated observations of the Crab in unbinned mode. On input we provide a source model XML file and the name of the source in the XML file that should be moved on a grid of test positions. @@ -31,11 +32,11 @@ nominal position of the Crab nebula with a grid spacing of 0.02 degrees. .. code-block:: bash $ cttsmap - Input event list, counts cube or observation definition file [events.fits] obs.xml + Input event list, counts cube or observation definition XML file [events.fits] obs.xml Calibration database [prod2] - Instrument response function [South_50h] - Source model [$CTOOLS/share/models/crab.xml] - Test source [Crab] + Instrument response function [South_0.5h] + Input model XML file [$CTOOLS/share/models/crab.xml] + Test source name [Crab] 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] @@ -43,14 +44,14 @@ nominal position of the Crab nebula with a grid spacing of 0.02 degrees. Image scale (in degrees/pixel) [0.02] Size of the X axis in pixels [200] 10 Size of the Y axis in pixels [200] 10 - Output Test Statistic map [tsmap.fits] + Output Test Statistic map file [tsmap.fits] :ref:`cttsmap` writes the Test Statistic map in the ``tsmap.fits`` file that contains one extension for the Test Statistic value and further extensions for the spectral parameters that have been fitted for the source at each position of the grid. The figure below show the Test Statistic map which reaches a maximum -value of 37484 near the centre of the map. +value of 40733 near the centre of the map. .. figure:: tsmap-crab.png :height: 400px diff --git a/doc/source/user_manual/getting_started/binning.rst b/doc/source/user_manual/getting_started/binning.rst new file mode 100644 index 0000000..35fac63 --- /dev/null +++ b/doc/source/user_manual/getting_started/binning.rst @@ -0,0 +1,83 @@ +.. _sec_binning_cta: + +Binning CTA data +~~~~~~~~~~~~~~~~ + +As next analysis step you will bin the data in 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 (typically logarithmically spaced, but this is under user +control). + +:ref:`ctbin` is executed by typing: + +.. 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] + +The counts cube will be centred on the location of the Crab (Right Ascension +83.63 degrees, Declination 22.01 degrees) and will be aligned in celestial +coordinates. A cartesian projection has been selected. The counts cube has +200 x 200 spatial pixels of 0.02 x 0.02 degrees in size, hence it covers a +total area of 4 x 4 degrees. + +The counts cube will contain 20 maps, which are logarithmically spaced +in energy, and which cover the energy range from 0.1 TeV to 100 TeV. In this +example, the counts cube will be saved as ``cntcube.fits`` in the working +directory. In addition to the counts cube, that is stored as the primary +image extension, the FITS file also contains an extension named ``EBOUNDS`` +that defines the energy boundaries that were used, and an extension ``GTI`` +that defines the Good Time Intervals that have been used. The following +image shows the resulting FITS file. The ``EBOUNDS`` table has 20 rows, one +for each energy bin, while the ``GTI`` table has just a single row, indicating +the start and stop time of the simulated data. + +.. figure:: cntmap-fits.jpg + :width: 600px + :align: center + + *Counts cube FITS file* + + +An image of the first bin, covering the energy range 100 - 141 GeV, is +shown below: + +.. figure:: cntmap-map.jpg + :height: 400px + :align: center + + *Counts cube for first energy bin* + +For illustration, the last few lines of the log file ``ctbin.log`` are +reproduced below: + +.. code-block:: xml + + 2015-12-07T20:52:41: +=================+ + 2015-12-07T20:52:41: | Bin observation | + 2015-12-07T20:52:41: +=================+ + 2015-12-07T20:52:41: === CTA observation === + 2015-12-07T20:52:41: Events in list ............: 22685 + 2015-12-07T20:52:41: Events in cube ............: 18195 + 2015-12-07T20:52:41: Event bins outside RoI ....: 0 + 2015-12-07T20:52:41: Events outside cube area ..: 4490 + 2015-12-07T20:52:41: Events outside energy bins : 0 + +From the 22685 events that have been simulated and stored in the +``events.fits`` file, 18195 lie within the cube boundaries and are thus put +into the resulting counts cube. The counts cube is stored in a cartesian +projection in a World Coordinate System (WCS) compliant format. diff --git a/doc/source/user_manual/getting_started/butterfly.rst b/doc/source/user_manual/getting_started/butterfly.rst new file mode 100644 index 0000000..26b2b02 --- /dev/null +++ b/doc/source/user_manual/getting_started/butterfly.rst @@ -0,0 +1,51 @@ +.. _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 + + python $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 fit* + \ No newline at end of file diff --git a/doc/source/user_manual/getting_started/fitting.rst b/doc/source/user_manual/getting_started/fitting.rst new file mode 100644 index 0000000..43e05ad --- /dev/null +++ b/doc/source/user_manual/getting_started/fitting.rst @@ -0,0 +1,195 @@ +.. _sec_fitting_cta: + +Fitting CTA data +~~~~~~~~~~~~~~~~ + +Now we are ready to fit the simulated data with a model. For simplicity +we use in this example the same model that we used to simulate the data +with :ref:`ctobssim`. Model fitting is done using the :ref:`ctlike` tool, +and we 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] + +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: + +.. code-block:: xml + + + + + + + + + + + + + + + + + + + + + + + +In this example, the ``Prefactor`` and ``Index`` of the spectral model for the +Crab 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 we 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. + +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: + +.. code-block:: xml + + 2015-12-07T20:55:32: +=================================+ + 2015-12-07T20:55:32: | Maximum likelihood optimisation | + 2015-12-07T20:55:32: +=================================+ + 2015-12-07T20:55:33: >Iteration 0: -logL=54877.412, Lambda=1.0e-03 + 2015-12-07T20:55:35: >Iteration 1: -logL=54875.025, Lambda=1.0e-03, delta=2.387, max(|grad|)=1.856683 [Index:7] + 2015-12-07T20:55:36: >Iteration 2: -logL=54875.024, Lambda=1.0e-04, delta=0.001, max(|grad|)=-0.024661 [Index:3] + 2015-12-07T20:55:36: + ... + 2015-12-07T20:55:38: +=========================================+ + 2015-12-07T20:55:38: | Maximum likelihood optimisation results | + 2015-12-07T20:55:38: +=========================================+ + 2015-12-07T20:55:38: === GOptimizerLM === + 2015-12-07T20:55:38: Optimized function value ..: 54875.024 + 2015-12-07T20:55:38: Absolute precision ........: 0.005 + 2015-12-07T20:55:38: Acceptable value decrease .: 2 + 2015-12-07T20:55:38: Optimization status .......: converged + 2015-12-07T20:55:38: Number of parameters ......: 10 + 2015-12-07T20:55:38: Number of free parameters .: 4 + 2015-12-07T20:55:38: Number of iterations ......: 2 + 2015-12-07T20:55:38: Lambda ....................: 1e-05 + 2015-12-07T20:55:38: Maximum log likelihood ....: -54875.024 + 2015-12-07T20:55:38: Observed events (Nobs) ...: 18195.000 + 2015-12-07T20:55:38: Predicted events (Npred) ..: 18194.999 (Nobs - Npred = 0.00146383) + 2015-12-07T20:55:38: === GModels === + 2015-12-07T20:55:38: Number of models ..........: 2 + 2015-12-07T20:55:38: Number of parameters ......: 10 + 2015-12-07T20:55:38: === GModelSky === + 2015-12-07T20:55:38: Name ......................: Crab + 2015-12-07T20:55:38: Instruments ...............: all + 2015-12-07T20:55:38: Instrument scale factors ..: unity + 2015-12-07T20:55:38: Observation identifiers ...: all + 2015-12-07T20:55:38: Model type ................: PointSource + 2015-12-07T20:55:38: Model components ..........: "SkyDirFunction" * "PowerLaw" * "Constant" + 2015-12-07T20:55:38: Number of parameters ......: 6 + 2015-12-07T20:55:38: Number of spatial par's ...: 2 + 2015-12-07T20:55:38: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1) + 2015-12-07T20:55:38: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1) + 2015-12-07T20:55:38: Number of spectral par's ..: 3 + 2015-12-07T20:55:38: Prefactor ................: 5.78443e-16 +/- 1.01348e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) + 2015-12-07T20:55:38: Index ....................: -2.50664 +/- 0.0153631 [-0,-5] (free,scale=-1,gradient) + 2015-12-07T20:55:38: PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) + 2015-12-07T20:55:38: Number of temporal par's ..: 1 + 2015-12-07T20:55:38: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) + 2015-12-07T20:55:38: === GCTAModelIrfBackground === + 2015-12-07T20:55:38: Name ......................: CTABackgroundModel + 2015-12-07T20:55:38: Instruments ...............: CTA + 2015-12-07T20:55:38: Instrument scale factors ..: unity + 2015-12-07T20:55:38: Observation identifiers ...: all + 2015-12-07T20:55:38: Model type ................: "PowerLaw" * "Constant" + 2015-12-07T20:55:38: Number of parameters ......: 4 + 2015-12-07T20:55:38: Number of spectral par's ..: 3 + 2015-12-07T20:55:38: Prefactor ................: 0.993018 +/- 0.0160113 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient) + 2015-12-07T20:55:38: Index ....................: 0.00269761 +/- 0.00973921 [-5,5] (free,scale=1,gradient) + 2015-12-07T20:55:38: PivotEnergy ..............: 1e+06 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) + 2015-12-07T20:55:38: Number of temporal par's ..: 1 + 2015-12-07T20:55:38: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) + +The maximum likelihood optimizer required 2 iterations to converge. This +is pretty fast, but recall that we 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 +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:: + + 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:: + + 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: + + .. code-block:: xml + + + + :ref:`ctlike` will then compute the Test Statistic value for that + source and dump the result in the log file: + + .. code-block:: xml + + 2015-12-07T20:58:45: === GModelSky === + 2015-12-07T20:58:45: Name ......................: Crab + 2015-12-07T20:58:45: Instruments ...............: all + 2015-12-07T20:58:45: Test Statistic ............: 21164.4 + + The Test Statistic value will also be added as new attribute + ``ts`` to the XML result file: + + .. code-block:: xml + + diff --git a/doc/source/user_manual/getting_started/quickstart.rst b/doc/source/user_manual/getting_started/quickstart.rst index 65fa257..5e8cc69 100644 --- a/doc/source/user_manual/getting_started/quickstart.rst +++ b/doc/source/user_manual/getting_started/quickstart.rst @@ -1,756 +1,23 @@ .. _quickstart: -A quick ctools tutorial ------------------------ - First steps with ctools -~~~~~~~~~~~~~~~~~~~~~~~ - -This section provides a step-by-step description of the CTA data analysis -procedure implemented in ctools. In this example we will simulate data for -the Crab nebula, observed on-axis, as a single source in the CTA field of -view. - -The example consists of the following analysis steps: - -- :ref:`sec_simulating_cta` - -- :ref:`sec_binning_cta` - -- :ref:`sec_fitting_cta` - -- :ref:`sec_unbinned_cta` - -- :ref:`sec_butterfly_cta` - - -.. _sec_simulating_cta: - -Simulating CTA data -~~~~~~~~~~~~~~~~~~~ - -CTA data are simulated using the :ref:`ctobssim` tool. To invoke the tool, -type :ref:`ctobssim` at the console prompt (which is denoted by ``$``). -:ref:`ctobssim` will query for a number of parameters: - -.. code-block:: bash - - $ ctobssim - 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_50h] - Model [$CTOOLS/share/models/crab.xml] - Output event data file or observation definition file [events.fits] - -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). -Otherwise, the specified value will overwrite the default value. -The round brackets ``( )`` indicate the range of possible parameter -values (if applicable). - - -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. - -.. note:: - - 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``. - -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 -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. - -The source and background model is defined by the XML file -``$CTOOLS/share/models/crab.xml``: - -.. code-block:: xml - - - - - - - - - - - - - - - - - - - - - - - -The model consists of a source library that contains two components: -the Crab nebula and an instrumental background model. - -The Crab nebula is modelled by a factorized sky model that has a spectral -and a spatial component (tags ```` and ````, -respectively). The spectrum is modelled by a power law, which is defined by -three parameters: the ``Prefactor``, the ``Index`` and the ``Scale``. -The spatial model has two parameters: Right Ascension in degrees (``RA``), and -Declination in degrees (``DEC``). Each parameter has a value and a scale factor, -the real value of the parameter being the product ``value * scale``. Typically, -``scale`` is chosen so that ``value`` is of the order of 1 (this is relevant for -model fitting). In addition, ``value`` is bound by a minimum (``min``) and -maximum (``max``) value, and a parameter may be free (``free="1"``) or fixed -(``free="0"``). The ``min``, ``max``, and ``free`` attributes are not -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 - -.. math:: - \frac{dN}{dE} = N_0 \left( \frac{E}{E_0} \right)^{\gamma} - -where the parameters in the XML definition have the following mappings: - -* :math:`N_0` = ``Prefactor`` -* :math:`\gamma` = ``Index`` -* :math:`E_0` = ``Scale`` - -.. 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 -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 -in the IRF. - -:ref:`ctobssim` has a couple of hidden parameters, the most important one being -certainly ``seed``. ``seed`` is an integer that specifies the seed value -for the random number generator, and changing this parameter will allow to -generate statistically independent Monte Carlo samples of CTA event data. -To use for example a seed value of 41 you should type: - -.. code-block:: bash - - $ ctobssim seed=41 - -:ref:`ctobssim` will write 2 files in the working directory: ``events.fits`` -and ``ctobssim.log``. The first file contains the simulated events in FITS -format and can be inspected using ``fv`` or ``ds9``. The FITS file will -contain three extensions: an empty primary image, a binary table named -``EVENTS`` that holds the events (one row per event), and a binary table -named ``GTI`` holding the Good Time Intervals (for the moment a single row -with 2 columns providing the start and the stop time of the simulated time -interval). - -The second file produced by :ref:`ctobssim` is a human readable log file that -contains information about the job execution. As example, the last lines -from this file are shown here: - -.. code-block:: xml - - 2015-05-22T19:24:53: +======================+ - 2015-05-22T19:24:53: | Simulate observation | - 2015-05-22T19:24:53: +======================+ - 2015-05-22T19:24:53: === Observation === - 2015-05-22T19:24:53: Simulation area ...........: 1.9635e+11 cm2 - 2015-05-22T19:24:53: Simulation cone ...........: RA=83.63 deg, Dec=22.01 deg, r=5.5 deg - 2015-05-22T19:24:53: Time interval .............: 0 - 1800 s - 2015-05-22T19:24:53: Photon energy range .......: 100 GeV - 100 TeV - 2015-05-22T19:24:53: Event energy range ........: 100 GeV - 100 TeV - 2015-05-22T19:24:53: MC source photons .........: 207547 [Crab] - 2015-05-22T19:24:53: MC source events ..........: 2961 [Crab] - 2015-05-22T19:24:53: MC source events ..........: 2961 (all source models) - 2015-05-22T19:24:55: MC background events ......: 13003 - 2015-05-22T19:24:55: MC events .................: 15964 (all models) - -Each line starts with the UTC time at which the line has been written. In -this run, 207547 Crab photons have been thrown over an area of 19.6 square -kilometres during a time interval of 1800 seconds. 2961 of these photons have -been registered by CTA as events. In the same time interval, 13003 background -events have been registred by CTA. - -You may change the name of the log file using the hidden parameter -``logfile``: - -.. code-block:: bash - - $ ctobssim logfile=my-private-log-file - -Furthermore, you may decide on the amount of information provided in the -log file (the chattiness of the executable) using the hidden parameter -``chatter``: - -.. code-block:: bash - - $ ctobssim chatter=4 - -``chatter`` can vary between 0 and 4, 0 providing no information while 4 -provides the most detailed information. **By default, all ctools have a -chatter level of 2.** - -You may also duplicate the log file information into the console by setting -the hidden ``debug`` parameter to yes: - -.. code-block:: bash - - $ ctobssim debug=yes - - -.. _sec_binning_cta: - -Binning CTA data -~~~~~~~~~~~~~~~~ - -As next analysis step you will bin the data in 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 (typically logarithmically spaced, but this is under user -control). - -:ref:`ctbin` is executed by typing: - -.. code-block:: bash - - $ ctbin - Event list or observation definition 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 [cntcube.fits] - -The counts cube will be centred on the location of the Crab (Right Ascension -83.63 degrees, Declination 22.01 degrees) and will be aligned in celestial -coordinates. A cartesian projection has been selected. The counts cube has -200 x 200 spatial pixels of 0.02 x 0.02 degrees in size, hence it covers a -total area of 4 x 4 degrees. - -The counts cube will contain 20 maps, which are logarithmically spaced -in energy, and which cover the energy range from 0.1 TeV to 100 TeV. In this -example, the counts cube will be saved as ``cntcube.fits`` in the working -directory. In addition to the counts cube, that is stored as the primary -image extension, the FITS file also contains an extension named ``EBOUNDS`` -that defines the energy boundaries that were used, and an extension ``GTI`` -that defines the Good Time Intervals that have been used. The following -image shows the resulting FITS file. The ``EBOUNDS`` table has 20 rows, one -for each energy bin, while the ``GTI`` table has just a single row, indicating -the start and stop time of the simulated data. - -.. figure:: cntmap-fits.jpg - :width: 600px - :align: center - - *Counts cube FITS file* - - -An image of the first bin, covering the energy range 100 - 141 GeV, is -shown below: - -.. figure:: cntmap-map.jpg - :height: 400px - :align: center - - *Counts cube for first energy bin* - - -For illustration, the last few lines of the log file ``ctbin.log`` are -reproduced below: - -.. code-block:: xml - - 2015-05-22T19:40:30: +=================+ - 2015-05-22T19:40:30: | Bin observation | - 2015-05-22T19:40:30: +=================+ - 2015-05-22T19:40:30: === Observation === - 2015-05-22T19:40:30: Events in list ............: 15964 - 2015-05-22T19:40:30: Events in cube ............: 13360 - 2015-05-22T19:40:30: Event bins outside RoI ....: 0 - 2015-05-22T19:40:30: Events outside cube area ..: 2604 - 2015-05-22T19:40:30: Events outside energy bins : 0 - -From the 15964 events that have been simulated and stored in the -``events.fits`` file, 13360 lie within the cube boundaries and are thus put -into the resulting counts cube. The counts cube is stored in a cartesian -projection in a World Coordinate System (WCS) compliant format. - - -.. _sec_fitting_cta: - -Fitting CTA data -~~~~~~~~~~~~~~~~ - -Now we are ready to fit the simulated data with a model. For simplicity -we use in this example the same model that we used to simulate the data -with :ref:`ctobssim`. Model fitting is done using the :ref:`ctlike` tool, -and we do the fit by typing: - -.. code-block:: bash - - $ ctlike - Event list, counts cube or observation definition file [events.fits] cntcube.fits - Exposure cube file (only needed for stacked analysis) [NONE] - PSF cube file (only needed for stacked analysis) [NONE] - Background cube file (only needed for stacked analysis) [NONE] - Calibration database [prod2] - Instrument response function [South_50h] - Source model [$CTOOLS/share/models/crab.xml] - Source model output 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: - -.. code-block:: xml - - - - - - - - - - - - - - - - - - - - - - - -In this example, the ``Prefactor`` and ``Index`` of the spectral model for the -Crab 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 we 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. - -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: - -.. code-block:: xml - - 2015-05-22T19:45:15: +=================================+ - 2015-05-22T19:45:15: | Maximum likelihood optimisation | - 2015-05-22T19:45:15: +=================================+ - 2015-05-22T19:45:16: >Iteration 0: -logL=39910.748, Lambda=1.0e-03 - 2015-05-22T19:45:18: >Iteration 1: -logL=39908.788, Lambda=1.0e-03, delta=1.960, max(|grad|)=6.701856 [Index:7] - 2015-05-22T19:45:19: >Iteration 2: -logL=39908.787, Lambda=1.0e-04, delta=0.001, max(|grad|)=-0.057041 [Index:3] - 2015-05-22T19:45:21: - 2015-05-22T19:45:21: +=========================================+ - 2015-05-22T19:45:21: | Maximum likelihood optimization results | - 2015-05-22T19:45:21: +=========================================+ - 2015-05-22T19:45:21: === GOptimizerLM === - 2015-05-22T19:45:21: Optimized function value ..: 39908.787 - 2015-05-22T19:45:21: Absolute precision ........: 0.005 - 2015-05-22T19:45:21: Acceptable value decrease .: 2 - 2015-05-22T19:45:21: Optimization status .......: converged - 2015-05-22T19:45:21: Number of parameters ......: 10 - 2015-05-22T19:45:21: Number of free parameters .: 4 - 2015-05-22T19:45:21: Number of iterations ......: 2 - 2015-05-22T19:45:21: Lambda ....................: 1e-05 - 2015-05-22T19:45:21: Maximum log likelihood ....: -39908.787 - 2015-05-22T19:45:21: Observed events (Nobs) ...: 13360.000 - 2015-05-22T19:45:21: Predicted events (Npred) ..: 13359.998 (Nobs - Npred = 0.00180799) - 2015-05-22T19:45:21: === GModels === - 2015-05-22T19:45:21: Number of models ..........: 2 - 2015-05-22T19:45:21: Number of parameters ......: 10 - 2015-05-22T19:45:21: === GModelSky === - 2015-05-22T19:45:21: Name ......................: Crab - 2015-05-22T19:45:21: Instruments ...............: all - 2015-05-22T19:45:21: Instrument scale factors ..: unity - 2015-05-22T19:45:21: Observation identifiers ...: all - 2015-05-22T19:45:21: Model type ................: PointSource - 2015-05-22T19:45:21: Model components ..........: "SkyDirFunction" * "PowerLaw" * "Constant" - 2015-05-22T19:45:21: Number of parameters ......: 6 - 2015-05-22T19:45:21: Number of spatial par's ...: 2 - 2015-05-22T19:45:21: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1) - 2015-05-22T19:45:21: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1) - 2015-05-22T19:45:21: Number of spectral par's ..: 3 - 2015-05-22T19:45:21: Prefactor ................: 5.78197e-16 +/- 1.10889e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) - 2015-05-22T19:45:21: Index ....................: -2.49174 +/- 0.0164369 [-0,-5] (free,scale=-1,gradient) - 2015-05-22T19:45:21: PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) - 2015-05-22T19:45:21: Number of temporal par's ..: 1 - 2015-05-22T19:45:21: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) - 2015-05-22T19:45:21: === GCTAModelIrfBackground === - 2015-05-22T19:45:21: Name ......................: CTABackgroundModel - 2015-05-22T19:45:21: Instruments ...............: CTA - 2015-05-22T19:45:21: Instrument scale factors ..: unity - 2015-05-22T19:45:21: Observation identifiers ...: all - 2015-05-22T19:45:21: Model type ................: "PowerLaw" * "Constant" - 2015-05-22T19:45:21: Number of parameters ......: 4 - 2015-05-22T19:45:21: Number of spectral par's ..: 3 - 2015-05-22T19:45:21: Prefactor ................: 1.02048 +/- 0.0280684 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient) - 2015-05-22T19:45:21: Index ....................: 0.0197856 +/- 0.0156686 [-5,5] (free,scale=1,gradient) - 2015-05-22T19:45:21: PivotEnergy ..............: 1e+06 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) - 2015-05-22T19:45:21: Number of temporal par's ..: 1 - 2015-05-22T19:45:21: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) - -The maximum likelihood optimizer required 2 iterations to converge. This -is pretty fast, but recall that we 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 -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:: - - 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:: - - 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: - - .. code-block:: xml - - - - :ref:`ctlike` will then compute the Test Statistic value for that - source and dump the result in the log file: - - .. code-block:: xml - - 2015-05-22T19:58:43: === GModelSky === - 2015-05-22T19:58:43: Name ......................: Crab - 2015-05-22T19:58:43: Instruments ...............: all - 2015-05-22T19:58:43: Test Statistic ............: 18662.6 - - The Test Statistic value will also be added as new attribute - ``ts`` to the XML result file: - - .. code-block:: xml - - - - -.. _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. We 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 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 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 - - *Data selection keywords* - -**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. - -Below some lines of the ``ctselect.log`` file that show the data selection -part: - -.. code-block:: xml - - 2015-05-22T20:01:46: +=================+ - 2015-05-22T20:01:46: | Event selection | - 2015-05-22T20:01:46: +=================+ - 2015-05-22T20:01:46: === Observation === - 2015-05-22T20:01:47: Selected energy range .....: 0.1 - 100 TeV - 2015-05-22T20:01:47: Requested ROI .............: Centre(RA,DEC)=(83.63, 22.01) deg, Radius=3 deg - 2015-05-22T20:01:47: ROI of data ...............: Centre(RA,DEC)=(83.63, 22.01) deg, Radius=5 deg - 2015-05-22T20:01:47: Selected ROI ..............: Centre(RA,DEC)=(83.63, 22.01) deg, Radius=3 deg - 2015-05-22T20:01:47: cfitsio selection .........: ENERGY >= 0.10000000 && ENERGY <= 100.00000000 && ANGSEP(83.630000,22.010000,RA,DEC) <= 3.000000 - 2015-05-22T20:01:47: FITS filename .............: /var/tmp/tmp.0.rRgLpO[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. - -.. warning:: - - :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. **If you use** :ref:`ctselect` **however to make a spatial or energy - selection, make sure that the counts cube is fully enclosed in the - selection intervals.** Otherwise you will get empty zones in the counts - cube of which the ctools are not aware of, and the subsequent analysis - results will be wrong. - -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: - -.. code-block:: bash - - $ ctlike - Event list, counts cube or observation definition file [cntcube.fits] selected_events.fits - Calibration database [prod2] - Instrument response function [South_50h] - Source model [$CTOOLS/share/models/crab.xml] - Source model output 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 15312 events, while the binned counts cube we used before had -200 x 200 x 20 = 800000 pixels. As unbinned maximum likelihood fitting loops -over the events (while binned maximum likelihood loops over the pixels), -there are much less operations to perform in unbinned than in binned mode -(there is some additional overhead in unbinned mode that comes from -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 100 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. - -Below you see the corresponding output from the ctlike.log file. The fitted -parameters are essentially identical to the ones found in binned mode. -The slight difference with respect to the binned analysis may be explained -by the different event sample that has been used for the analysis: while -binned likelihood works on rectangular counts cubes, unbinned likelihood works -on circular event selection regions. It is thus not possible to select exactly -the same events for both analyses. - -.. code-block:: xml - - 2015-05-22T20:05:54: +=================================+ - 2015-05-22T20:05:54: | Maximum likelihood optimisation | - 2015-05-22T20:05:54: +=================================+ - 2015-05-22T20:05:54: >Iteration 0: -logL=93854.543, Lambda=1.0e-03 - 2015-05-22T20:05:54: >Iteration 1: -logL=93848.595, Lambda=1.0e-03, delta=5.948, max(|grad|)=22.480581 [Index:7] - 2015-05-22T20:05:54: >Iteration 2: -logL=93848.587, Lambda=1.0e-04, delta=0.008, max(|grad|)=-0.081846 [Index:3] - 2015-05-22T20:05:54: >Iteration 3: -logL=93848.587, Lambda=1.0e-05, delta=0.000, max(|grad|)=-0.003485 [Index:3] - 2015-05-22T20:05:54: - 2015-05-22T20:05:54: +=========================================+ - 2015-05-22T20:05:54: | Maximum likelihood optimization results | - 2015-05-22T20:05:54: +=========================================+ - 2015-05-22T20:05:54: === GOptimizerLM === - 2015-05-22T20:05:54: Optimized function value ..: 93848.587 - 2015-05-22T20:05:54: Absolute precision ........: 0.005 - 2015-05-22T20:05:54: Acceptable value decrease .: 2 - 2015-05-22T20:05:54: Optimization status .......: converged - 2015-05-22T20:05:54: Number of parameters ......: 10 - 2015-05-22T20:05:54: Number of free parameters .: 4 - 2015-05-22T20:05:54: Number of iterations ......: 3 - 2015-05-22T20:05:54: Lambda ....................: 1e-06 - 2015-05-22T20:05:54: Maximum log likelihood ....: -93848.587 - 2015-05-22T20:05:54: Observed events (Nobs) ...: 15312.000 - 2015-05-22T20:05:54: Predicted events (Npred) ..: 15312.000 (Nobs - Npred = 3.03823e-06) - 2015-05-22T20:05:54: === GModels === - 2015-05-22T20:05:54: Number of models ..........: 2 - 2015-05-22T20:05:54: Number of parameters ......: 10 - 2015-05-22T20:05:54: === GModelSky === - 2015-05-22T20:05:54: Name ......................: Crab - 2015-05-22T20:05:54: Instruments ...............: all - 2015-05-22T20:05:54: Instrument scale factors ..: unity - 2015-05-22T20:05:54: Observation identifiers ...: all - 2015-05-22T20:05:54: Model type ................: PointSource - 2015-05-22T20:05:54: Model components ..........: "SkyDirFunction" * "PowerLaw" * "Constant" - 2015-05-22T20:05:54: Number of parameters ......: 6 - 2015-05-22T20:05:54: Number of spatial par's ...: 2 - 2015-05-22T20:05:54: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1) - 2015-05-22T20:05:54: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1) - 2015-05-22T20:05:54: Number of spectral par's ..: 3 - 2015-05-22T20:05:54: Prefactor ................: 5.80619e-16 +/- 1.11149e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) - 2015-05-22T20:05:54: Index ....................: -2.49282 +/- 0.0163467 [-0,-5] (free,scale=-1,gradient) - 2015-05-22T20:05:54: PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) - 2015-05-22T20:05:54: Number of temporal par's ..: 1 - 2015-05-22T20:05:54: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) - 2015-05-22T20:05:54: === GCTAModelIrfBackground === - 2015-05-22T20:05:54: Name ......................: CTABackgroundModel - 2015-05-22T20:05:54: Instruments ...............: CTA - 2015-05-22T20:05:54: Instrument scale factors ..: unity - 2015-05-22T20:05:54: Observation identifiers ...: all - 2015-05-22T20:05:54: Model type ................: "PowerLaw" * "Constant" - 2015-05-22T20:05:54: Number of parameters ......: 4 - 2015-05-22T20:05:54: Number of spectral par's ..: 3 - 2015-05-22T20:05:54: Prefactor ................: 1.03511 +/- 0.024311 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient) - 2015-05-22T20:05:54: Index ....................: 0.0333176 +/- 0.0134386 [-5,5] (free,scale=1,gradient) - 2015-05-22T20:05:54: PivotEnergy ..............: 1e+06 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) - 2015-05-22T20:05:54: Number of temporal par's ..: 1 - 2015-05-22T20:05:54: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) - - -.. _sec_butterfly_cta: - -Calculating and visualising a butterfly -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -To visualise the analysis results retrieved above, one can calculate the -confidence band of the spectral fit. -The :ref:`ctbutterfly` tool takes the optimised source model as input. -It takes the covariance matrix from the fit to conduct a -Gaussian error propagation for each energy value. It will write the -butterfly information into an ASCII file. -The following example shows how to compute such a butterfly from the command -line. - -.. code-block:: bash - - $ ctbutterfly - Input event list, counts cube or observation definition file [events.fits] - Calibration database [prod2] - Instrument response function [South_50h] - Source model [$CTOOLS/share/models/crab.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 - - python $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 fit* - \ No newline at end of file +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. Finally we will show how +you can create a butterfly diagram that shows the confidence band of the +fitted spectral model. + +.. toctree:: + :maxdepth: 2 + + simulating + binning + fitting + unbinned + butterfly diff --git a/doc/source/user_manual/getting_started/resmap-crab.png b/doc/source/user_manual/getting_started/resmap-crab.png index b98097f..a5840fb 100644 Binary files a/doc/source/user_manual/getting_started/resmap-crab.png and b/doc/source/user_manual/getting_started/resmap-crab.png differ diff --git a/doc/source/user_manual/getting_started/response.rst b/doc/source/user_manual/getting_started/response.rst index 8fcb40a..b0213c8 100644 --- a/doc/source/user_manual/getting_started/response.rst +++ b/doc/source/user_manual/getting_started/response.rst @@ -99,7 +99,7 @@ Specifying the response function as input parameters ctools that require instrument response functions have two parameters to specify the calibration database name and the response function name. The following example shows a ``ctobssim`` run using the ``prod2`` -calibration database and the ``South_50h`` response function: +calibration database and the ``South_0.5h`` response function: .. code-block:: bash @@ -112,8 +112,9 @@ calibration database and the ``South_50h`` response function: Lower energy limit (TeV) [0.1] Upper energy limit (TeV) [100.0] Calibration database [prod2] - Instrument response function [South_50h] - Model [$CTOOLS/share/models/crab.xml] + 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] Running the other tools is equivalent. @@ -131,7 +132,7 @@ definition file: - + @@ -141,9 +142,9 @@ response name. You can then pass this file directly to, e.g., ``ctlike``: .. code-block:: bash $ ctlike - Event list, counts cube or observation definition file [events.fits] obs_irf.xml - Source model [$CTOOLS/share/models/crab.xml] - Source model output file [crab_results.xml] + Input event list, counts cube or observation definition XML file [events.fits] obs_irf.xml + Input model XML file [$CTOOLS/share/models/crab.xml] + Output model XML file [crab_results.xml] Note that ``ctlike`` does not ask for the calibration database and response name as it found the relevant information in the XML file. @@ -181,11 +182,11 @@ and response name from within Python: import gammalib obs = gammalib.GCTAObservation() caldb = gammalib.GCaldb("cta", "prod2") - irf = "South_50h" + irf = "South_0.5h" obs.response(irf, caldb) The calibration database is set by creating a ``GCaldb`` object. The constructor takes as argument the mission (always ``cta``) and the database name, in our case ``prod2``. The response function is then set -by passing the response name (here ``South_50h``) and the calibration +by passing the response name (here ``South_0.5h``) and the calibration database object to the ``response`` method. diff --git a/doc/source/user_manual/getting_started/simulating.rst b/doc/source/user_manual/getting_started/simulating.rst new file mode 100644 index 0000000..32d42c5 --- /dev/null +++ b/doc/source/user_manual/getting_started/simulating.rst @@ -0,0 +1,210 @@ +.. _sec_simulating_cta: + +Simulating CTA data +~~~~~~~~~~~~~~~~~~~ + +You can simulate a CTA observation with the :ref:`ctobssim` tool. + +To invoke the tool, type :ref:`ctobssim` at the console prompt +(which is denoted by ``$``). +:ref:`ctobssim` will query for a number of parameters: + +.. code-block:: bash + + $ ctobssim + 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] + +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). +Otherwise, the specified value will overwrite the default value. +The round brackets ``( )`` indicate the range of possible parameter +values (if applicable). + + +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. + +.. note:: + + 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``. + +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 +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. + +The source and background model is defined by the XML file +``$CTOOLS/share/models/crab.xml``: + +.. code-block:: xml + + + + + + + + + + + + + + + + + + + + + + + +The model consists of a source library that contains two components: +the Crab nebula and an instrumental background model. + +The Crab nebula is modelled by a factorized sky model that has a spectral +and a spatial component (tags ```` and ````, +respectively). The spectrum is modelled by a power law, which is defined by +three parameters: the ``Prefactor``, the ``Index`` and the ``Scale``. +The spatial model has two parameters: Right Ascension in degrees (``RA``), and +Declination in degrees (``DEC``). Each parameter has a value and a scale factor, +the real value of the parameter being the product ``value * scale``. Typically, +``scale`` is chosen so that ``value`` is of the order of 1 (this is relevant for +model fitting). In addition, ``value`` is bound by a minimum (``min``) and +maximum (``max``) value, and a parameter may be free (``free="1"``) or fixed +(``free="0"``). The ``min``, ``max``, and ``free`` attributes are not +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 + +.. math:: + \frac{dN}{dE} = N_0 \left( \frac{E}{E_0} \right)^{\gamma} + +where the parameters in the XML definition have the following mappings: + +* :math:`N_0` = ``Prefactor`` +* :math:`\gamma` = ``Index`` +* :math:`E_0` = ``Scale`` + +.. 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 +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 +in the IRF. + +:ref:`ctobssim` has a couple of hidden parameters, the most important one being +certainly ``seed``. ``seed`` is an integer that specifies the seed value +for the random number generator, and changing this parameter will allow to +generate statistically independent Monte Carlo samples of CTA event data. +To use for example a seed value of 41 you should type: + +.. code-block:: bash + + $ ctobssim seed=41 + +:ref:`ctobssim` will write 2 files in the working directory: ``events.fits`` +and ``ctobssim.log``. The first file contains the simulated events in FITS +format and can be inspected using ``fv`` or ``ds9``. The FITS file will +contain three extensions: an empty primary image, a binary table named +``EVENTS`` that holds the events (one row per event), and a binary table +named ``GTI`` holding the Good Time Intervals (for the moment a single row +with 2 columns providing the start and the stop time of the simulated time +interval). + +The second file produced by :ref:`ctobssim` is a human readable log file that +contains information about the job execution. As example, the last lines +from this file are shown here: + +.. code-block:: xml + + 2015-12-07T20:48:30: +======================+ + 2015-12-07T20:48:30: | Simulate observation | + 2015-12-07T20:48:30: +======================+ + 2015-12-07T20:48:30: === CTA observation === + 2015-12-07T20:48:30: Simulation area ...........: 1.9635e+11 cm2 + 2015-12-07T20:48:30: Simulation cone ...........: RA=83.63 deg, Dec=22.01 deg, r=5.5 deg + 2015-12-07T20:48:30: Time interval .............: 0 - 1800 s + 2015-12-07T20:48:30: Photon energy range .......: 100 GeV - 100 TeV + 2015-12-07T20:48:30: Event energy range ........: 100 GeV - 100 TeV + 2015-12-07T20:48:30: Use model .................: Crab + 2015-12-07T20:48:30: Normalization .............: 1 [Crab] + 2015-12-07T20:48:30: Flux ......................: 5.87297e-10 [Crab] photons/cm2/s + 2015-12-07T20:48:30: Normalized flux ...........: 5.87297e-10 [Crab] photons/cm2/s + 2015-12-07T20:48:30: Photon rate ...............: 115.316 photons/s [Crab] + 2015-12-07T20:48:31: MC source photons .........: 207547 [Crab] + 2015-12-07T20:48:31: MC source events ..........: 3646 [Crab] + 2015-12-07T20:48:31: MC source events ..........: 3646 (all source models) + 2015-12-07T20:48:33: MC events outside ROI .....: 0 + 2015-12-07T20:48:33: MC background events ......: 19039 + 2015-12-07T20:48:33: MC events .................: 22685 (all models) + +Each line starts with the UTC time at which the line has been written. In +this run, 207547 Crab photons have been thrown over an area of 19.6 square +kilometres during a time interval of 1800 seconds. 3646 of these photons have +been registered by CTA as events. In the same time interval, 19039 background +events have been registred by CTA. + +You may change the name of the log file using the hidden parameter +``logfile``: + +.. code-block:: bash + + $ ctobssim logfile=my-private-log-file + +Furthermore, you may decide on the amount of information provided in the +log file (the chattiness of the executable) using the hidden parameter +``chatter``: + +.. code-block:: bash + + $ ctobssim chatter=4 + +``chatter`` can vary between 0 and 4, 0 providing no information while 4 +provides the most detailed information. + +**By default, all ctools have a chatter level of 2.** + +You may also duplicate the log file information into the console by setting +the hidden ``debug`` parameter to yes: + +.. code-block:: bash + + $ ctobssim debug=yes diff --git a/doc/source/user_manual/getting_started/tsmap-crab.png b/doc/source/user_manual/getting_started/tsmap-crab.png index f1ef431..5af71f6 100644 Binary files a/doc/source/user_manual/getting_started/tsmap-crab.png and b/doc/source/user_manual/getting_started/tsmap-crab.png differ diff --git a/doc/source/user_manual/getting_started/unbinned.rst b/doc/source/user_manual/getting_started/unbinned.rst new file mode 100644 index 0000000..636758e --- /dev/null +++ b/doc/source/user_manual/getting_started/unbinned.rst @@ -0,0 +1,200 @@ +.. _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. We 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 + + *Data selection keywords* + +**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. + +Below some lines of the ``ctselect.log`` file that show the data selection +part: + +.. code-block:: xml + + 2015-12-06T23:44:52: +=================+ + 2015-12-06T23:44:52: | Event selection | + 2015-12-06T23:44:52: +=================+ + 2015-12-06T23:44:52: === CTA observation === + 2015-12-06T23:44:52: Selected energy range .....: 0.1 - 100 TeV + 2015-12-06T23:44:52: Requested ROI .............: Centre(RA,DEC)=(83.63, 22.01) deg, Radius=3 deg + 2015-12-06T23:44:52: ROI of data ...............: Centre(RA,DEC)=(83.63, 22.01) deg, Radius=5 deg + 2015-12-06T23:44:52: Selected ROI ..............: Centre(RA,DEC)=(83.63, 22.01) deg, Radius=3 deg + 2015-12-06T23:44:52: cfitsio selection .........: ENERGY >= 0.10000000 && ENERGY <= 100.00000000 && ANGSEP(83.630000,22.010000,RA,DEC) <= 3.000000 + 2015-12-06T23:44:52: FITS filename .............: /var/tmp/tmp.0.jHg5hJ[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. + +.. warning:: + + :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. **If you use** :ref:`ctselect` **however to make a spatial or energy + selection, make sure that the counts cube is fully enclosed in the + selection intervals.** Otherwise you will get empty zones in the counts + cube of which the ctools are not aware of, and the subsequent analysis + results will be wrong. + +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: + +.. 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] + +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 21203 events, while the binned counts cube we used before had +200 x 200 x 20 = 800000 pixels. As unbinned maximum likelihood fitting loops +over the events (while binned maximum likelihood loops over the pixels), +there are much less operations to perform in unbinned than in binned mode +(there is some additional overhead in unbinned mode that comes from +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 100 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. + +Below you see the corresponding output from the ctlike.log file. The fitted +parameters are essentially identical to the ones found in binned mode. +The slight difference with respect to the binned analysis may be explained +by the different event sample that has been used for the analysis: while +binned likelihood works on rectangular counts cubes, unbinned likelihood works +on circular event selection regions. It is thus not possible to select exactly +the same events for both analyses. + +.. code-block:: xml + + 2015-12-07T21:02:17: +=================================+ + 2015-12-07T21:02:17: | Maximum likelihood optimisation | + 2015-12-07T21:02:17: +=================================+ + 2015-12-07T21:02:17: >Iteration 0: -logL=135278.644, Lambda=1.0e-03 + 2015-12-07T21:02:17: >Iteration 1: -logL=135274.475, Lambda=1.0e-03, delta=4.168, max(|grad|)=5.806873 [Index:7] + 2015-12-07T21:02:17: >Iteration 2: -logL=135274.474, Lambda=1.0e-04, delta=0.001, max(|grad|)=-0.009463 [Index:7] + 2015-12-07T21:02:17: + ... + 2015-12-07T21:02:17: +=========================================+ + 2015-12-07T21:02:17: | Maximum likelihood optimisation results | + 2015-12-07T21:02:17: +=========================================+ + 2015-12-07T21:02:17: === GOptimizerLM === + 2015-12-07T21:02:17: Optimized function value ..: 135274.474 + 2015-12-07T21:02:17: Absolute precision ........: 0.005 + 2015-12-07T21:02:17: Acceptable value decrease .: 2 + 2015-12-07T21:02:17: Optimization status .......: converged + 2015-12-07T21:02:17: Number of parameters ......: 10 + 2015-12-07T21:02:17: Number of free parameters .: 4 + 2015-12-07T21:02:17: Number of iterations ......: 2 + 2015-12-07T21:02:17: Lambda ....................: 1e-05 + 2015-12-07T21:02:17: Maximum log likelihood ....: -135274.474 + 2015-12-07T21:02:17: Observed events (Nobs) ...: 21203.000 + 2015-12-07T21:02:17: Predicted events (Npred) ..: 21202.998 (Nobs - Npred = 0.00198937) + 2015-12-07T21:02:17: === GModels === + 2015-12-07T21:02:17: Number of models ..........: 2 + 2015-12-07T21:02:17: Number of parameters ......: 10 + 2015-12-07T21:02:17: === GModelSky === + 2015-12-07T21:02:17: Name ......................: Crab + 2015-12-07T21:02:17: Instruments ...............: all + 2015-12-07T21:02:17: Instrument scale factors ..: unity + 2015-12-07T21:02:17: Observation identifiers ...: all + 2015-12-07T21:02:17: Model type ................: PointSource + 2015-12-07T21:02:17: Model components ..........: "SkyDirFunction" * "PowerLaw" * "Constant" + 2015-12-07T21:02:17: Number of parameters ......: 6 + 2015-12-07T21:02:17: Number of spatial par's ...: 2 + 2015-12-07T21:02:17: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1) + 2015-12-07T21:02:17: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1) + 2015-12-07T21:02:17: Number of spectral par's ..: 3 + 2015-12-07T21:02:17: Prefactor ................: 5.81583e-16 +/- 1.01727e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) + 2015-12-07T21:02:17: Index ....................: -2.50611 +/- 0.0153542 [-0,-5] (free,scale=-1,gradient) + 2015-12-07T21:02:17: PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) + 2015-12-07T21:02:17: Number of temporal par's ..: 1 + 2015-12-07T21:02:17: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) + 2015-12-07T21:02:17: === GCTAModelIrfBackground === + 2015-12-07T21:02:17: Name ......................: CTABackgroundModel + 2015-12-07T21:02:17: Instruments ...............: CTA + 2015-12-07T21:02:17: Instrument scale factors ..: unity + 2015-12-07T21:02:17: Observation identifiers ...: all + 2015-12-07T21:02:17: Model type ................: "PowerLaw" * "Constant" + 2015-12-07T21:02:17: Number of parameters ......: 4 + 2015-12-07T21:02:17: Number of spectral par's ..: 3 + 2015-12-07T21:02:17: Prefactor ................: 0.990866 +/- 0.0136406 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient) + 2015-12-07T21:02:17: Index ....................: 0.00536014 +/- 0.00839916 [-5,5] (free,scale=1,gradient) + 2015-12-07T21:02:17: PivotEnergy ..............: 1e+06 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) + 2015-12-07T21:02:17: Number of temporal par's ..: 1 + 2015-12-07T21:02:17: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) diff --git a/doc/source/user_manual/glossary.rst b/doc/source/user_manual/glossary.rst index f6ccf4f..2a10e51 100644 --- a/doc/source/user_manual/glossary.rst +++ b/doc/source/user_manual/glossary.rst @@ -1,7 +1,7 @@ Glossary ======== -- **GammaLib** (``_) +- **GammaLib** (``_) A versatile toolbox for scientific analysis of astronomical gamma-ray data. diff --git a/doc/source/user_manual/index.rst b/doc/source/user_manual/index.rst index da83bea..6edf832 100644 --- a/doc/source/user_manual/index.rst +++ b/doc/source/user_manual/index.rst @@ -8,8 +8,6 @@ User Manual introduction getting getting_started/index - enhancing - expert plotting glossary diff --git a/doc/source/user_manual/introduction.rst b/doc/source/user_manual/introduction.rst index cec202d..5f2367f 100644 --- a/doc/source/user_manual/introduction.rst +++ b/doc/source/user_manual/introduction.rst @@ -22,6 +22,13 @@ also be specified as arguments separated by a whitespace to the ctool: $ ctobssim ra=83.63 dec=22.01 rad=5.0 tmin=0 tmax=1800 emin=0.1 emax=100.0 caldb=prod2 irf=South_50h inmodel=$CTOOLS/share/models/crab.xml outevents=events.fits +If you need help about the usage of a ctool or a cscript, type the name of +the ctool or cscript followed by the ``--help`` option: + +.. code-block:: bash + + $ ctobssim --help + ctools includes also a Python module that declares every ctool as a Python class. The following example illustrates the usage: @@ -47,7 +54,7 @@ command line or used as Python classes, as illustrated below: In this User Manual you will now learn how to install and use the ctools package. -If you encounter any problems, please read the :ref:`issues` section. -If you're desperate, you may need :ref:`help`. +If you encounter any problems, please read the :ref:`faq` and :ref:`issues` +sections. If you're desperate, you may need :ref:`help`. And you should read the :ref:`develop` section if you would like to contribute to the ctools developments. diff --git a/doc/source/user_manual/plotting.rst b/doc/source/user_manual/plotting.rst index 908038b..bbf6abf 100644 --- a/doc/source/user_manual/plotting.rst +++ b/doc/source/user_manual/plotting.rst @@ -1,46 +1,56 @@ Plotting ======== -Gammalib and ctools currently do not include classes and scripts to make plots. +Gammalib and ctools do not include classes and scripts to make plots. -This is intentional because creating and maintaining plotting tools that work on most user's machines -and cover most user's needs is a big project in itself, and deciding to support a -plotting package that is popular now as part of the CTA science tools could become very problematic in the future. +This is intentional because creating and maintaining plotting tools that +work on most user's machines and cover most user's needs is a big project +in itself, and deciding to support a plotting package that is popular now +as part of the CTA science tools could become very problematic in the +future. + +But we realise that most ctools users need to create plots for their +high-level science results, so on this page we provide a few pointers +to external packages and tools. -But we realise that most ctools users need to create plots for their high-level science results, -so on this page we provide a few pointers to external packages and tools. General ------- -The most common way to create publication-quality astronomy plots at the moment is to use the `matplotlib`_ -Python plotting package. As illustrated by the `matplotlib gallery`_, virtually any plot can be constructed, -but matplotlib doesn't have built-in easy-to-use classes or functions to create the common plots for gamma-ray -astronomy, which are sky maps and spectra, so in the sections below we give some specific pointers for those. +The most common way to create publication-quality astronomy plots at the +moment is to use the `matplotlib`_ Python plotting package. As illustrated +by the `matplotlib gallery`_, virtually any plot can be constructed, +but matplotlib doesn't have built-in easy-to-use classes or functions to +create the common plots for gamma-ray astronomy, which are sky maps and +spectra, so in the sections below we give some specific pointers for those. + Sky maps -------- -For creating publication quality sky map plots via a Python script a good choice is to use `APLPy`_, `WCSAxes`_ -or `kapteyn`_, which are all based on matplotlib. +For creating publication quality sky map plots via a Python script a good +choice is to use `APLPy`_, `WCSAxes`_ or `kapteyn`_, which are all based on +matplotlib. -You can also use interactive astronomy image viewers like `DS9`_, `Aladin`_ or `Ginga`_ -which allow you to save the image to a PNG or PDF file after customising it interactively, -but also contain command line options to create the plot, i.e. you can write a bash or Python script -to make the plot reproducible (DS9 and Ginga also feature direct communication with Python, -but that is a bit harder to use and usually used for interactive image analysis, not plotting scripts). +You can also use interactive astronomy image viewers like `DS9`_, `Aladin`_ +or `Ginga`_ which allow you to save the image to a PNG or PDF file after +customising it interactively, but also contain command line options to create +the plot, i.e. you can write a bash or Python script to make the plot +reproducible (DS9 and Ginga also feature direct communication with Python, +but that is a bit harder to use and usually used for interactive image +analysis, not plotting scripts). Spectra and light curves ------------------------ -To plot spectra or spectral energy distributions (including flux model curves, flux model error bands, -flux points and residuals) and light curves, you can use matplotlib directly, -or use the classes and scripts in `Gammapy`_. - -Another option is to use virtual observatory interactive spectrum tools like `IRIS`_, `CASSIS`_ or `VOSpec`_ -which allow you to save the plot to a PNG or PDF file after creating it interactively. +To plot spectra or spectral energy distributions (including flux model curves, +flux model error bands, flux points and residuals) and light curves, you can +use matplotlib directly, or use the classes and scripts in `Gammapy`_. +Another option is to use virtual observatory interactive spectrum tools +like `IRIS`_, `CASSIS`_ or `VOSpec`_ which allow you to save the plot to a +PNG or PDF file after creating it interactively. .. _DS9: http://ds9.si.edu/ .. _Aladin: http://aladin.u-strasbg.fr/ -- libgit2 0.21.4