Thierry Reposeur / ctools

Commit e4831aad988cf7c8df2fa8ddd081b27f9d5017e5

Authored by Jürgen Knödlseder
1 parent ed8025d2

Write analysis tutorial that goes beyond the initial steps

Inline Side-by-side
Showing 3 changed files with 613 additions and 8 deletions
doc/source/user_manual/getting_started/beyond.rst
  View file @e4831aa
  1 +.. _beyond:
  2 +
1 3 Beyond the first steps
2 4 ----------------------
3 5  
  6 +Assuming that you have read the :ref:`quickstart <quickstart>` tutorial,
  7 +here now some information on how to make the following more complex
  8 +analyses:
  9 +
  10 +- :ref:`sec_combining_obs`
  11 +
  12 +- :ref:`sec_stacked`
  13 +
  14 +- :ref:`sec_tsmap`
  15 +
  16 +- :ref:`sec_residual`
  17 +
  18 +
  19 +.. _sec_combining_obs:
  20 +
4 21 Combining observations
5 22 ~~~~~~~~~~~~~~~~~~~~~~
6 23  
  24 +Generally, the CTA data you may want to analyse will not only be composed of
  25 +a single observation (a.k.a. run) but of a list of observations that should
  26 +be combined in a joint analysis.
  27 +ctools has the capability to collect individual observations in a list and
  28 +to perform for example a joint maximum likelihood fit of all observations
  29 +in a single shot.
  30 +Here is an example that illustrates how to do that.
  31 +
  32 +Let's start with the simulation of two 30 min long observations of the Crab
  33 +nebula, each offset by 0.5 deg from the nebula in opposite directions:
  34 +
  35 +.. code-block:: bash
  36 +
  37 + $ ctobssim
  38 + RA of pointing (degrees) (0-360) [83.63]
  39 + Dec of pointing (degrees) (-90-90) [22.01] 21.51
  40 + Radius of FOV (degrees) (0-180) [5.0]
  41 + Start time (MET in s) (0) [0.0]
  42 + End time (MET in s) (0) [1800.0]
  43 + Lower energy limit (TeV) (0) [0.1]
  44 + Upper energy limit (TeV) (0) [100.0]
  45 + Calibration database [dummy]
  46 + Instrument response function [cta_dummy_irf]
  47 + Model [$CTOOLS/share/models/crab.xml]
  48 + Output event data file or observation definition file [events.fits] events1.fits
  49 +
  50 +.. code-block:: bash
  51 +
  52 + $ ctobssim
  53 + RA of pointing (degrees) (0-360) [83.63]
  54 + Dec of pointing (degrees) (-90-90) [21.51] 22.51
  55 + Radius of FOV (degrees) (0-180) [5.0]
  56 + Start time (MET in s) (0) [0.0]
  57 + End time (MET in s) (0) [1800.0]
  58 + Lower energy limit (TeV) (0) [0.1]
  59 + Upper energy limit (TeV) (0) [100.0]
  60 + Calibration database [dummy]
  61 + Instrument response function [cta_dummy_irf]
  62 + Model [$CTOOLS/share/models/crab.xml]
  63 + Output event data file or observation definition file [events1.fits] events2.fits
  64 +
  65 +This will produce the two event files ``events1.fits`` and ``events2.fits``
  66 +on disk.
  67 +
  68 +As next step we have to create an **observation definition XML file** that
  69 +collects both observations in a list.
  70 +Here is how such a file would look like:
  71 +
  72 +.. code-block:: xml
  73 +
  74 + <?xml version="1.0" standalone="no"?>
  75 + <observation_list title="observation library">
  76 + <observation name="Crab" id="00001" instrument="CTA">
  77 + <parameter name="EventList" file="events1.fits"/>
  78 + </observation>
  79 + <observation name="Crab" id="00002" instrument="CTA">
  80 + <parameter name="EventList" file="events2.fits"/>
  81 + </observation>
  82 + </observation_list>
  83 +
  84 +The file contains a single ``<observation_list>`` tag that contains
  85 +two ``<observation>`` tags that each define an observation.
  86 +Each observation has a ``name``, an ``id`` and an ``instrument``
  87 +attribute.
  88 +The ``name`` attribute can have any arbitrary value, and may be the
  89 +same for all observations.
  90 +**However, the** ``id`` **attribute needs to be a unique character
  91 +string for any given instrument.**
  92 +The ``instrument`` attribute is a case-sensitive string that identifies
  93 +the instrument with which the observation was taken.
  94 +**Please make sure that the instrument string is set correctly so that
  95 +ctools knows which instrument specific functions need to be called.**
  96 +
  97 +.. note::
  98 +
  99 + The instrument string for a CTA observation is obviously ``CTA``.
  100 + In case that data from any existing Imaging Air Cherenkov Telescope
  101 + are provided in the correct format, the instrument strings
  102 + ``HESS``, ``MAGIC`` and ``VERITAS`` are also recognised. This opens
  103 + the obvious possibility to combine data from different
  104 + Imaging Air Cherenkov Telescope in a joint analysis.
  105 + **Please recall that instrument strings are case sensitive.**
  106 +
  107 +Now we are ready to do a joint maximum likelihood analysis using
  108 +:ref:`ctlike`:
  109 +
  110 +.. code-block:: bash
  111 +
  112 + $ ctlike
  113 + Event list, counts cube or observation definition file [events.fits] obs.xml
  114 + Calibration database [dummy]
  115 + Instrument response function [cta_dummy_irf]
  116 + Source model [$CTOOLS/share/models/crab.xml]
  117 + Source model output file [crab_results.xml]
  118 +
  119 +Instead of providing an event list or a counts cube, we now provided the
  120 +filename of the observation definition XML file (here ``obs.xml``) as input
  121 +parameter.
  122 +:ref:`ctlike` recognises this format and automatically performs a
  123 +joint analysis, as indicated in the log file:
  124 +
  125 +.. code-block:: xml
  126 +
  127 + 2015-02-04T20:33:26: +==============+
  128 + 2015-02-04T20:33:26: | Observations |
  129 + 2015-02-04T20:33:26: +==============+
  130 + 2015-02-04T20:33:26: === GObservations ===
  131 + 2015-02-04T20:33:26: Number of observations ....: 2
  132 + 2015-02-04T20:33:26: Number of predicted events : 0
  133 + 2015-02-04T20:33:26: === GCTAObservation ===
  134 + 2015-02-04T20:33:26: Name ......................: Crab
  135 + 2015-02-04T20:33:26: Identifier ................: 00001
  136 + 2015-02-04T20:33:26: Instrument ................: CTA
  137 + 2015-02-04T20:33:26: Event file ................: events1.fits
  138 + 2015-02-04T20:33:26: Event type ................: EventList
  139 + 2015-02-04T20:33:26: Statistics ................: Poisson
  140 + 2015-02-04T20:33:26: Ontime ....................: 1800 s
  141 + 2015-02-04T20:33:26: Livetime ..................: 1710 s
  142 + 2015-02-04T20:33:26: Deadtime correction .......: 0.95
  143 + 2015-02-04T20:33:26: User energy range .........: undefined
  144 + 2015-02-04T20:33:26: === GCTAPointing ===
  145 + 2015-02-04T20:33:26: Pointing direction ........: (RA,Dec)=(83.63,21.51)
  146 + ...
  147 + 2015-02-04T20:33:26: === GCTAObservation ===
  148 + 2015-02-04T20:33:26: Name ......................: Crab
  149 + 2015-02-04T20:33:26: Identifier ................: 00002
  150 + 2015-02-04T20:33:26: Instrument ................: CTA
  151 + 2015-02-04T20:33:26: Event file ................: events2.fits
  152 + 2015-02-04T20:33:26: Event type ................: EventList
  153 + 2015-02-04T20:33:26: Statistics ................: Poisson
  154 + 2015-02-04T20:33:26: Ontime ....................: 1800 s
  155 + 2015-02-04T20:33:26: Livetime ..................: 1710 s
  156 + 2015-02-04T20:33:26: Deadtime correction .......: 0.95
  157 + 2015-02-04T20:33:26: User energy range .........: undefined
  158 + 2015-02-04T20:33:26: === GCTAPointing ===
  159 + 2015-02-04T20:33:26: Pointing direction ........: (RA,Dec)=(83.63,22.51)
  160 + ...
  161 + 2015-02-04T20:33:26: +=================================+
  162 + 2015-02-04T20:33:26: | Maximum likelihood optimisation |
  163 + 2015-02-04T20:33:26: +=================================+
  164 + 2015-02-04T20:33:26: >Iteration 0: -logL=89099.788, Lambda=1.0e-03
  165 + 2015-02-04T20:33:26: >Iteration 1: -logL=89094.435, Lambda=1.0e-03, delta=5.353, max(|grad|)=-4.388872 [Normalization:7]
  166 + 2015-02-04T20:33:26: >Iteration 2: -logL=89094.420, Lambda=1.0e-04, delta=0.015, max(|grad|)=-0.116432 [Index:3]
  167 + 2015-02-04T20:33:26: >Iteration 3: -logL=89094.420, Lambda=1.0e-05, delta=0.000, max(|grad|)=-0.003989 [Index:3]
  168 + 2015-02-04T20:33:26:
  169 + 2015-02-04T20:33:26: +=========================================+
  170 + 2015-02-04T20:33:26: | Maximum likelihood optimization results |
  171 + 2015-02-04T20:33:26: +=========================================+
  172 + 2015-02-04T20:33:26: === GOptimizerLM ===
  173 + 2015-02-04T20:33:26: Optimized function value ..: 89094.420
  174 + 2015-02-04T20:33:26: Absolute precision ........: 0.005
  175 + 2015-02-04T20:33:26: Acceptable value decrease .: 2
  176 + 2015-02-04T20:33:26: Optimization status .......: converged
  177 + 2015-02-04T20:33:26: Number of parameters ......: 9
  178 + 2015-02-04T20:33:26: Number of free parameters .: 4
  179 + 2015-02-04T20:33:26: Number of iterations ......: 3
  180 + 2015-02-04T20:33:26: Lambda ....................: 1e-06
  181 + 2015-02-04T20:33:26: Maximum log likelihood ....: -89094.420
  182 + 2015-02-04T20:33:26: Observed events (Nobs) ...: 12260.000
  183 + 2015-02-04T20:33:26: Predicted events (Npred) ..: 12260.000 (Nobs - Npred = 6.95764e-06)
  184 + 2015-02-04T20:33:26: === GModels ===
  185 + 2015-02-04T20:33:26: Number of models ..........: 2
  186 + 2015-02-04T20:33:26: Number of parameters ......: 9
  187 + 2015-02-04T20:33:26: === GModelSky ===
  188 + 2015-02-04T20:33:26: Name ......................: Crab
  189 + 2015-02-04T20:33:26: Instruments ...............: all
  190 + 2015-02-04T20:33:26: Instrument scale factors ..: unity
  191 + 2015-02-04T20:33:26: Observation identifiers ...: all
  192 + 2015-02-04T20:33:26: Model type ................: PointSource
  193 + 2015-02-04T20:33:26: Model components ..........: "SkyDirFunction" * "PowerLaw" * "Constant"
  194 + 2015-02-04T20:33:26: Number of parameters ......: 6
  195 + 2015-02-04T20:33:26: Number of spatial par's ...: 2
  196 + 2015-02-04T20:33:26: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1)
  197 + 2015-02-04T20:33:26: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1)
  198 + 2015-02-04T20:33:26: Number of spectral par's ..: 3
  199 + 2015-02-04T20:33:26: Prefactor ................: 6.06322e-16 +/- 1.44948e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
  200 + 2015-02-04T20:33:26: Index ....................: -2.49814 +/- 0.0178524 [-0,-5] (free,scale=-1,gradient)
  201 + 2015-02-04T20:33:26: PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient)
  202 + 2015-02-04T20:33:26: Number of temporal par's ..: 1
  203 + 2015-02-04T20:33:26: Constant .................: 1 (relative value) (fixed,scale=1,gradient)
  204 + 2015-02-04T20:33:26: === GCTAModelRadialAcceptance ===
  205 + 2015-02-04T20:33:26: Name ......................: Background
  206 + 2015-02-04T20:33:26: Instruments ...............: CTA
  207 + 2015-02-04T20:33:26: Instrument scale factors ..: unity
  208 + 2015-02-04T20:33:26: Observation identifiers ...: all
  209 + 2015-02-04T20:33:26: Model type ................: "Gaussian" * "FileFunction" * "Constant"
  210 + 2015-02-04T20:33:26: Number of parameters ......: 3
  211 + 2015-02-04T20:33:26: Number of radial par's ....: 1
  212 + 2015-02-04T20:33:26: Sigma ....................: 3.0356 +/- 0.0215259 [0.01,10] deg2 (free,scale=1,gradient)
  213 + 2015-02-04T20:33:26: Number of spectral par's ..: 1
  214 + 2015-02-04T20:33:26: Normalization ............: 0.99936 +/- 0.0122024 [0,1000] (free,scale=1,gradient)
  215 + 2015-02-04T20:33:26: Number of temporal par's ..: 1
  216 + 2015-02-04T20:33:26: Constant .................: 1 (relative value) (fixed,scale=1,gradient)
  217 + 2015-02-04T20:33:26:
  218 + 2015-02-04T20:33:26: +==============+
  219 + 2015-02-04T20:33:26: | Save results |
  220 + 2015-02-04T20:33:26: +==============+
  221 + 2015-02-04T20:33:26:
  222 + 2015-02-04T20:33:26: Application "ctlike" terminated after 10 wall clock seconds, consuming 0.385432 seconds of CPU time.
  223 +
  224 +The log file indicates that the fit converged quickly, the spectral
  225 +parameters of the Crab nebula have now been constrained using the events
  226 +from both observations.
  227 +The computation time increases roughly linearly with the number of
  228 +observations that are combined, although ctools implements parallel
  229 +multi-core processing which will spread the likelihood computation for
  230 +the different observations over all CPU cores that are available.
  231 +**Doing a joint unbinned analysis is thus an efficient solution if
  232 +data from multiple observations should be combined.**
  233 +
  234 +Combining observations is not limited to unbinned data (i.e. event lists)
  235 +but may also be applied to binned data (i.e. counts cubes).
  236 +Using :ref:`ctbin` we can create counts cubes from both event lists which
  237 +may then be combined in an observation definition XML file:
  238 +
  239 +.. code-block:: xml
  240 +
  241 + <?xml version="1.0" standalone="no"?>
  242 + <observation_list title="observation library">
  243 + <observation name="Crab" id="00001" instrument="CTA">
  244 + <parameter name="CountsCube" file="cntcube1.fits"/>
  245 + </observation>
  246 + <observation name="Crab" id="00002" instrument="CTA">
  247 + <parameter name="CountsCube" file="cntcube2.fits"/>
  248 + </observation>
  249 + </observation_list>
  250 +
  251 +Feeding the observation definition XML file to :ref:`ctlike` will then
  252 +lead to a joint binned analysis.
  253 +In the joint binned analysis, the events of individual observations are
  254 +not combined, but are kept separate in distinct counts cubes.
  255 +This is not very efficient, as generally counts cubes for short duration
  256 +observations are only sparsly populated and the likelihood computation
  257 +has to loop over a hugh number of data space bins (though also here
  258 +:ref:`ctlike` benefits from multi-core parallel processing).
  259 +**Though possible, a joint binned analysis is thus not the recommended
  260 +method for combining observations.**
  261 +An alternative is to stack the events of all observations into a single
  262 +counts cube.
  263 +The :ref:`following section <sec_stacked>` describes how such a stacked
  264 +analysis is done with ctools.
  265 +
  266 +.. note::
  267 +
  268 + Given that logic, unbinned and binned observations may also be combined
  269 + in a joint analysis, although this Use Case may be a bit academic:
  270 +
  271 + .. code-block:: xml
  272 +
  273 + <?xml version="1.0" standalone="no"?>
  274 + <observation_list title="observation library">
  275 + <observation name="Crab" id="00001" instrument="CTA">
  276 + <parameter name="EventList" file="events1.fits"/>
  277 + </observation>
  278 + <observation name="Crab" id="00002" instrument="CTA">
  279 + <parameter name="CountsCube" file="cntcube2.fits"/>
  280 + </observation>
  281 + </observation_list>
  282 +
  283 +
  284 +.. _sec_stacked:
7 285  
8 286 Performing a stacked analysis
9   -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
10 287 \ No newline at end of file
  288 +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  289 +
  290 +A stacked analysis is a binned analysis where all data from multiple
  291 +observations are stacked into a single counts cube.
  292 +The event stacking is done using the :ref:`ctbin` tool.
  293 +Instead of providing to :ref:`ctbin` an event list we
  294 +now specify the observation definition XML file ``obs.xml``
  295 +on input.
  296 +:ref:`ctbin` will then loop over all observations and collect all events
  297 +into a single counts cube:
  298 +
  299 +.. code-block:: bash
  300 +
  301 + $ ctbin
  302 + Event list or observation definition file [events2.fits] obs.xml
  303 + First coordinate of image center in degrees (RA or galactic l) [83.63]
  304 + Second coordinate of image center in degrees (DEC or galactic b) [22.01]
  305 + Projection method e.g. AIT|AZP|CAR|MER|STG|TAN (AIT|AZP|CAR|MER|STG|TAN) [CAR]
  306 + Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL]
  307 + Image scale (in degrees/pixel) [0.02]
  308 + Size of the X axis in pixels [200]
  309 + Size of the Y axis in pixels [200]
  310 + Algorithm for defining energy bins (FILE|LIN|LOG) [LOG]
  311 + Start value for first energy bin in TeV [0.1]
  312 + Stop value for last energy bin in TeV [100.0]
  313 + Number of energy bins [20]
  314 + Output counts cube [cntcube2.fits] cntcube.fits
  315 +
  316 +We now have a stacked counts cube ``cntcube.fits`` on disk.
  317 +Before we can use that counts cube in a maximum likelihood
  318 +analysis, we have to compute the instrument response and the
  319 +background model that are needed to describe the stacked data.
  320 +For the former, we have to compute the total exposure for the stacked
  321 +cube (i.e. the sum of the effective areas for each observation multiplied
  322 +by the corresponding lifetimes) and an effective point spread function
  323 +(i.e. the point spread function of the different observations weighted by
  324 +the corresponding exposures).
  325 +To get both informations we use the :ref:`ctexpcube` and
  326 +:ref:`ctpsfcube` tools:
  327 +
  328 +.. code-block:: bash
  329 +
  330 + $ ctexpcube
  331 + Event list or observation definition file [NONE] obs.xml
  332 + Calibration database [dummy]
  333 + Instrument response function [cta_dummy_irf]
  334 + Counts cube for exposure cube definition [NONE] cntcube.fits
  335 + Output exposure cube file [expcube.fits]
  336 +
  337 +.. code-block:: bash
  338 +
  339 + $ ctpsfcube
  340 + Event list or observation definition file [NONE] obs.xml
  341 + Calibration database [dummy]
  342 + Instrument response function [cta_dummy_irf]
  343 + Counts cube for psf cube definition [NONE]
  344 + First coordinate of image center in degrees (RA or galactic l) [83.63]
  345 + Second coordinate of image center in degrees (DEC or galactic b) [22.01]
  346 + Projection method e.g. AIT|AZP|CAR|MER|MOL|STG|TAN (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR]
  347 + Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL]
  348 + Image scale (in degrees/pixel) [1.0]
  349 + Size of the X axis in pixels [10]
  350 + Size of the Y axis in pixels [10]
  351 + Start value for first energy bin in TeV [0.1]
  352 + Stop value for last energy bin in TeV [100.0]
  353 + Number of energy bins [20]
  354 + Output psf cube file [psfcube.fits]
  355 +
  356 +We provide the ``obs.xml`` file on input to inform both tools which
  357 +observations have been combined.
  358 +For :ref:`ctexpcube` we further provide the counts cube so that the
  359 +tool can copy the exposure cube definition (number of spatial pixels
  360 +and pixel size, number of energy bins) from the counts cube.
  361 +This minimises the number of further user parameters that need to be
  362 +provided and assures an exposure cube that is compatible with the counts
  363 +cube.
  364 +For :ref:`ctpsfcube` we do not use the counts cube for the PSF cube
  365 +definition as this would lead to a hugh file owing to the fine spatial
  366 +pixelisation of the counts cube.
  367 +Since the PSF evolves only slowly over the field of fiew, we provide a
  368 +rather coarse spatial binning of 1 degree covering a grid of 10 x 10
  369 +degrees around the centre of the counts cube.
  370 +For the energy binning, we use the same logarithmic binning that has
  371 +also been used for the counts cube.
  372 +
  373 +As final step of the analysis preparation, we need to generate a
  374 +background cube using the :ref:`ctbkgcube` tool:
  375 +
  376 +.. code-block:: bash
  377 +
  378 + $ ctbkgcube
  379 + Input event list or observation definition file [NONE] obs.xml
  380 + Calibration database [dummy]
  381 + Instrument response function [cta_dummy_irf]
  382 + Input (background) model XML file [NONE] $CTOOLS/share/models/crab.xml
  383 + Counts cube for background cube definition [NONE] cntcube.fits
  384 + Output background cube file [bkgcube.fits]
  385 + Output (background) model XML file [NONE] model.xml
  386 +
  387 +The usage of :ref:`ctbkgcube` is very similar to that of :ref:`ctexpcube`,
  388 +yet it takes the model XML file as an additional input parameter.
  389 +We here use the usual ``$CTOOLS/share/models/crab.xml`` Crab plus
  390 +background model file that is shipped with the ctools.
  391 +:ref:`ctbkgcube` provides on output the background cube file
  392 +``bkgcube.fits`` and the model XML file ``model.xml`` that can
  393 +be used for further analysis.
  394 +Having a look at ``model.xml`` illustrates how the background
  395 +modelling works:
  396 +
  397 +.. code-block:: xml
  398 +
  399 + <?xml version="1.0" encoding="UTF-8" standalone="no"?>
  400 + <source_library title="source library">
  401 + <source name="Crab" type="PointSource" tscalc="0">
  402 + <spectrum type="PowerLaw">
  403 + <parameter name="Prefactor" value="5.7" error="0" scale="1e-16" min="1e-07" max="1000" free="1" />
  404 + <parameter name="Index" value="2.48" error="0" scale="-1" min="0" max="5" free="1" />
  405 + <parameter name="Scale" value="0.3" scale="1e+06" min="0.01" max="1000" free="0" />
  406 + </spectrum>
  407 + <spatialModel type="SkyDirFunction">
  408 + <parameter name="RA" value="83.6331" scale="1" min="-360" max="360" free="0" />
  409 + <parameter name="DEC" value="22.0145" scale="1" min="-90" max="90" free="0" />
  410 + </spatialModel>
  411 + </source>
  412 + <source name="ctbkgcube default background model" type="CTACubeBackground">
  413 + <spectrum type="PowerLaw">
  414 + <parameter name="Prefactor" value="1" error="0" scale="1" min="0" free="1" />
  415 + <parameter name="Index" value="0" error="0" scale="1" min="-10" max="10" free="1" />
  416 + <parameter name="Scale" value="1" scale="1e+06" free="0" />
  417 + </spectrum>
  418 + <spatialModel type="MapCubeFunction" file="bkgcube.fits">
  419 + <parameter name="Normalization" value="1" scale="1" min="0.001" max="1000" free="0" />
  420 + </spatialModel>
  421 + </source>
  422 + </source_library>
  423 +
  424 +The Crab source component is the same that is also present in
  425 +``$CTOOLS/share/models/crab.xml`` and is not modified.
  426 +The background component, however, has been replaced and now is
  427 +the ``ctbkgcube default background model``.
  428 +This model is of type ``CTACubeBackground`` which is a 3-dimensional
  429 +data cube that describes the expected background rate as function
  430 +of spatial position and energy.
  431 +The data cube is multiplied by a power law spectrum that allows to adjust
  432 +the normalization and slope of the background spectrum in the fit.
  433 +This power law could be replaced by any spectral model that is found
  434 +as an appropriate multiplicator to the background cube.
  435 +
  436 +.. note::
  437 +
  438 + There is no constraint on providing the same spatial binning or
  439 + the same energy binning for an exposure cube, a PSF cube,
  440 + a background cube and a counts cube.
  441 + ctools interpolates internally the exposure cube, PSF cube and
  442 + background cube values, hence any arbitrary appropriate binning may
  443 + be used.
  444 + Using the same binning for the exposure cube, the background cube and
  445 + the counts cube is only a convenience.
  446 +
  447 +Now we have all files at hand to launch a stacked maximum likelihood
  448 +analysis using the :ref:`ctlike` tool:
  449 +
  450 +.. code-block:: bash
  451 +
  452 + $ ctlike
  453 + Event list, counts cube or observation definition file [events.fits] cntcube.fits
  454 + Exposure cube file [NONE] expcube.fits
  455 + PSF cube file [NONE] psfcube.fits
  456 + Source model [$CTOOLS/share/models/crab.xml] model.xml
  457 + Source model output file [crab_results.xml]
  458 +
  459 +:ref:`ctlike` recognises that a counts cube should be analysed and queries
  460 +for the exposure cube and PSF cube file names.
  461 +We specified the names of the files produced by the :ref:`ctexpcube` and
  462 +:ref:`ctpsfcube` tools and furthermore provided the ``model.xml`` file
  463 +that is generated by the :ref:`ctbkgcube` tool as source model.
  464 +The log file of the :ref:`ctlike` run is shown below.
  465 +Note that the spectral model that is multiplied with the background
  466 +cube has a Prefactor of 1.06 +/- 0.02 and an Index of 0.004 +/- 0.009,
  467 +indicating a very small correction to the actual spectrum of the background
  468 +cube.
  469 +Real life situations may of course require larger correction factors.
  470 +
  471 +.. code-block:: xml
  472 +
  473 + 2015-02-04T23:00:25: +=================================+
  474 + 2015-02-04T23:00:25: | Maximum likelihood optimisation |
  475 + 2015-02-04T23:00:25: +=================================+
  476 + 2015-02-04T23:00:29: >Iteration 0: -logL=36748.112, Lambda=1.0e-03
  477 + 2015-02-04T23:00:33: >Iteration 1: -logL=36734.170, Lambda=1.0e-03, delta=13.942, max(|grad|)=34.291911 [Index:8]
  478 + 2015-02-04T23:00:38: >Iteration 2: -logL=36734.127, Lambda=1.0e-04, delta=0.044, max(|grad|)=0.178088 [Index:8]
  479 + 2015-02-04T23:00:42: >Iteration 3: -logL=36734.127, Lambda=1.0e-05, delta=0.000, max(|grad|)=-0.001993 [Index:8]
  480 + 2015-02-04T23:00:46:
  481 + 2015-02-04T23:00:46: +=========================================+
  482 + 2015-02-04T23:00:46: | Maximum likelihood optimization results |
  483 + 2015-02-04T23:00:46: +=========================================+
  484 + 2015-02-04T23:00:46: === GOptimizerLM ===
  485 + 2015-02-04T23:00:46: Optimized function value ..: 36734.127
  486 + 2015-02-04T23:00:46: Absolute precision ........: 0.005
  487 + 2015-02-04T23:00:46: Acceptable value decrease .: 2
  488 + 2015-02-04T23:00:46: Optimization status .......: converged
  489 + 2015-02-04T23:00:46: Number of parameters ......: 11
  490 + 2015-02-04T23:00:46: Number of free parameters .: 4
  491 + 2015-02-04T23:00:46: Number of iterations ......: 3
  492 + 2015-02-04T23:00:46: Lambda ....................: 1e-06
  493 + 2015-02-04T23:00:46: Maximum log likelihood ....: -36734.127
  494 + 2015-02-04T23:00:46: Observed events (Nobs) ...: 10685.000
  495 + 2015-02-04T23:00:46: Predicted events (Npred) ..: 10685.000 (Nobs - Npred = 2.06522e-06)
  496 + 2015-02-04T23:00:46: === GModels ===
  497 + 2015-02-04T23:00:46: Number of models ..........: 2
  498 + 2015-02-04T23:00:46: Number of parameters ......: 11
  499 + 2015-02-04T23:00:46: === GModelSky ===
  500 + 2015-02-04T23:00:46: Name ......................: Crab
  501 + 2015-02-04T23:00:46: Instruments ...............: all
  502 + 2015-02-04T23:00:46: Instrument scale factors ..: unity
  503 + 2015-02-04T23:00:46: Observation identifiers ...: all
  504 + 2015-02-04T23:00:46: Model type ................: PointSource
  505 + 2015-02-04T23:00:46: Model components ..........: "SkyDirFunction" * "PowerLaw" * "Constant"
  506 + 2015-02-04T23:00:46: Number of parameters ......: 6
  507 + 2015-02-04T23:00:46: Number of spatial par's ...: 2
  508 + 2015-02-04T23:00:46: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1)
  509 + 2015-02-04T23:00:46: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1)
  510 + 2015-02-04T23:00:46: Number of spectral par's ..: 3
  511 + 2015-02-04T23:00:46: Prefactor ................: 6.01471e-16 +/- 1.44183e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
  512 + 2015-02-04T23:00:46: Index ....................: -2.49533 +/- 0.0179769 [-0,-5] (free,scale=-1,gradient)
  513 + 2015-02-04T23:00:46: PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient)
  514 + 2015-02-04T23:00:46: Number of temporal par's ..: 1
  515 + 2015-02-04T23:00:46: Constant .................: 1 (relative value) (fixed,scale=1,gradient)
  516 + 2015-02-04T23:00:46: === GCTAModelCubeBackground ===
  517 + 2015-02-04T23:00:46: Name ......................: ctbkgcube default background model
  518 + 2015-02-04T23:00:46: Instruments ...............: all
  519 + 2015-02-04T23:00:46: Instrument scale factors ..: unity
  520 + 2015-02-04T23:00:46: Observation identifiers ...: all
  521 + 2015-02-04T23:00:46: Model type ................: "MapCubeFunction" * "PowerLaw" * "Constant"
  522 + 2015-02-04T23:00:46: Number of parameters ......: 5
  523 + 2015-02-04T23:00:46: Number of spatial par's ...: 1
  524 + 2015-02-04T23:00:46: Normalization ............: 1 [0.001,1000] (fixed,scale=1,gradient)
  525 + 2015-02-04T23:00:46: Number of spectral par's ..: 3
  526 + 2015-02-04T23:00:46: Prefactor ................: 1.05876 +/- 0.0178132 [0,infty[ ph/cm2/s/MeV (free,scale=1,gradient)
  527 + 2015-02-04T23:00:46: Index ....................: 0.00396962 +/- 0.00881842 [-10,10] (free,scale=1,gradient)
  528 + 2015-02-04T23:00:46: PivotEnergy ..............: 1e+06 MeV (fixed,scale=1e+06,gradient)
  529 + 2015-02-04T23:00:46: Number of temporal par's ..: 1
  530 + 2015-02-04T23:00:46: Constant .................: 1 (relative value) (fixed,scale=1,gradient)
  531 + 2015-02-04T23:00:46:
  532 + 2015-02-04T23:00:46: +==============+
  533 + 2015-02-04T23:00:46: | Save results |
  534 + 2015-02-04T23:00:46: +==============+
  535 + 2015-02-04T23:00:46:
  536 + 2015-02-04T23:00:46: Application "ctlike" terminated after 43 wall clock seconds, consuming 21.5624 seconds of CPU time.
  537 +
  538 +
  539 +.. _sec_tsmap:
  540 +
  541 +Generating a Test Statistics map
  542 +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  543 +
  544 +In the above examples we were always considering that we knew exactly
  545 +where the source is located, but in real life situations this may not
  546 +necessarily be the case.
  547 +We thus need a mean to search for significant excess emission in our data.
  548 +One possibility to accomplish this task is to compute the Test Statistics
  549 +for a specific source at a grid of trial positions.
  550 +The :ref:`cttsmap` tool will exactly do that job.
  551 +The Test Statistics is defined as twice the log-likelihood difference
  552 +between fitting a source at a given position on top of a model or fitting
  553 +no source.
  554 +Roughly speaken, the square root of the Test Statistics value gives
  555 +the source detection significance in Gaussian sigmas, although the
  556 +exact releation depends somewhat on the formulation of the statistical
  557 +problem.
  558 +
  559 +Here an example where we apply the :ref:`cttsmap` tool to our two
  560 +simulated observations of the Crab in unbinned mode.
  561 +On input we provide a source model XML file and the name of the source
  562 +in the XML file that should be moved on a grid of test positions.
  563 +This means that the position of the source given in the XML file is
  564 +in fact ignored and replaced by the grid positions that are defined by
  565 +the user parameters.
  566 +In the example we defined a grid of 10 x 10 positions around the
  567 +nominal position of the Crab nebula with a grid spacing of 0.02 degrees.
  568 +
  569 +.. code-block:: bash
  570 +
  571 + $ cttsmap
  572 + Input event list, counts cube or observation definition file [events.fits] obs.xml
  573 + Calibration database [dummy]
  574 + Instrument response function [cta_dummy_irf]
  575 + Source model [$CTOOLS/share/models/crab.xml]
  576 + Test source [Crab]
  577 + First coordinate of image center in degrees (RA or galactic l) [83.63]
  578 + Second coordinate of image center in degrees (DEC or galactic b) [22.01]
  579 + Projection method e.g. AIT|AZP|CAR|MER|STG|TAN (AIT|AZP|CAR|MER|STG|TAN) [CAR]
  580 + Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL]
  581 + Image scale (in degrees/pixel) [0.02]
  582 + Size of the X axis in pixels [200] 10
  583 + Size of the Y axis in pixels [200] 10
  584 + Output Test Statistics map [tsmap.fits]
  585 +
  586 +:ref:`cttsmap` writes the Test Statistics map in the ``tsmap.fits`` file
  587 +that contains one extension for the Test Statistics value and further
  588 +extensions for the spectral parameters that have been fitted for the
  589 +source at each position of the grid.
  590 +The figure below show the Test Statistics map which reaches a maximum
  591 +value of 13782 near the centre of the map.
  592 +
  593 +.. figure:: tsmap-crab.png
  594 + :width: 60%
  595 +
  596 + *Test Statistics map of the Crab region*
  597 +
  598 +.. note::
  599 +
  600 + Alike :ref:`ctlike`, the :ref:`cttsmap` tool works either for unbinned,
  601 + binned or stacked analysis. On input it takes either a single event
  602 + list, a single counts cube, or an observation definition XML file that
  603 + allows a joint analysis of multiple observations.
  604 +
  605 +
  606 +.. _sec_residual:
  607 +
  608 +Generating a residual map
  609 +~~~~~~~~~~~~~~~~~~~~~~~~~
  610 +
  611 +TBW
... ...
doc/source/user_manual/getting_started/quickstart.rst
  View file @e4831aa
  1 +.. _quickstart:
  2 +
1 3 A quick ctools tutorial
2 4 -----------------------
3 5  
... ... @@ -242,7 +244,7 @@ control).
242 244 Start value for first energy bin in TeV [0.1]
243 245 Stop value for last energy bin in TeV [100.0]
244 246 Number of energy bins [20]
245   - Output counts cube [cntmap.fits]
  247 + Output counts cube [cntcube.fits]
246 248  
247 249 The counts cube will be centred on the location of the Crab (Right Ascension
248 250 83.63 degrees, Declination 22.01 degrees) and will be aligned in celestial
... ... @@ -252,7 +254,7 @@ total area of 4 x 4 degrees.
252 254  
253 255 The counts cube will contain 20 maps, which are logarithmically spaced
254 256 in energy, and which cover the energy range from 0.1 TeV to 100 TeV. In this
255   -example, the counts cube will be saved as ``cntmap.fits`` in the working
  257 +example, the counts cube will be saved as ``cntcube.fits`` in the working
256 258 directory. In addition to the counts cube, that is stored as the primary
257 259 image extension, the FITS file also contains an extension named ``EBOUNDS``
258 260 that defines the energy boundaries that were used, and an extension ``GTI``
... ... @@ -291,8 +293,8 @@ reproduced below:
291 293 2015-02-04T15:33:34: Name ......................:
292 294 2015-02-04T15:33:34: Identifier ................:
293 295 2015-02-04T15:33:34: Instrument ................: CTA
294   - 2015-02-04T15:33:34: Event file ................: events.fits
295   - 2015-02-04T15:33:34: Event type ................: CountsMap
  296 + 2015-02-04T15:33:34: Event file ................:
  297 + 2015-02-04T15:33:34: Event type ................: CountsCube
296 298 2015-02-04T15:33:34: Statistics ................: Poisson
297 299 2015-02-04T15:33:34: Ontime ....................: 1800 s
298 300 2015-02-04T15:33:34: Livetime ..................: 1710 s
... ... @@ -333,7 +335,9 @@ and we do the fit by typing:
333 335 .. code-block:: bash
334 336  
335 337 $ ctlike
336   - Event list, counts cube or observation definition file [events.fits] cntmap.fits
  338 + Event list, counts cube or observation definition file [events.fits] cntcube.fits
  339 + Exposure cube file [NONE]
  340 + PSF cube file [NONE]
337 341 Calibration database [dummy]
338 342 Instrument response function [cta_dummy_irf]
339 343 Source model [$CTOOLS/share/models/crab.xml]
... ... @@ -695,8 +699,8 @@ the same events for both analyses.
695 699  
696 700 .. _sec_butterfly_cta:
697 701  
698   -Calculate and visualise butterfly
699   -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  702 +Calculating and visualising a butterfly
  703 +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
700 704  
701 705 To visualise the analysis results retrieved above, one can calculate the
702 706 confidence band of the spectral fit.
... ...
doc/source/user_manual/getting_started/tsmap-crab.png 0 → 100644
View file @e4831aa

75.2 KB