From 0197267a917e68e694dad41bcd11b221f321ff73 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Thu, 9 Jan 2014 23:21:00 +0100 Subject: [PATCH] Copy over quickstart from ctools web page and update slightly the examples for the actual ctools version --- doc/source/user_manual/cntmap-fits.jpg | Bin 0 -> 48122 bytes doc/source/user_manual/cntmap-map.jpg | Bin 0 -> 23229 bytes doc/source/user_manual/dskeys.jpg | Bin 0 -> 102007 bytes doc/source/user_manual/quickstart.rst | 654 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 654 insertions(+), 0 deletions(-) create mode 100644 doc/source/user_manual/cntmap-fits.jpg create mode 100644 doc/source/user_manual/cntmap-map.jpg create mode 100644 doc/source/user_manual/dskeys.jpg create mode 100644 doc/source/user_manual/quickstart.rst diff --git a/doc/source/user_manual/cntmap-fits.jpg b/doc/source/user_manual/cntmap-fits.jpg new file mode 100644 index 0000000..f9c73aa Binary files /dev/null and b/doc/source/user_manual/cntmap-fits.jpg differ diff --git a/doc/source/user_manual/cntmap-map.jpg b/doc/source/user_manual/cntmap-map.jpg new file mode 100644 index 0000000..b22de89 Binary files /dev/null and b/doc/source/user_manual/cntmap-map.jpg differ diff --git a/doc/source/user_manual/dskeys.jpg b/doc/source/user_manual/dskeys.jpg new file mode 100644 index 0000000..0185466 Binary files /dev/null and b/doc/source/user_manual/dskeys.jpg differ diff --git a/doc/source/user_manual/quickstart.rst b/doc/source/user_manual/quickstart.rst new file mode 100644 index 0000000..28e73cb --- /dev/null +++ b/doc/source/user_manual/quickstart.rst @@ -0,0 +1,654 @@ +Getting started +=============== + +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` + + +.. _sec_simulating_cta: + +Simulating CTA data +~~~~~~~~~~~~~~~~~~~ + +CTA data are simulated by the executable ctobssim. To start the executable, +type ctobssim at the console prompt (which is denoted by ``$``). ctobssim +will query for a number of parameters: + +.. code-block:: bash + + $ ctobssim + Model [test/data/crab.xml] $GAMMALIB/share/models/crab.xml + Calibration database [test/irf] $GAMMALIB/share/caldb/cta + Instrument response function [cta_dummy_irf] + RA of pointing (degrees) (0-360) [83.63] + Dec of pointing (degrees) (-90-90) [22.01] + Radius of FOV (degrees) (0-180) [10.0] + Start time (MET in s) (0) [0.0] + End time (MET in s) (0) [1800.0] + Lower energy limit (TeV) (0) [0.1] + Upper energy limit (TeV) (0) [100.0] + Output event data file or observation definition file [events.fits] + +Each line represents a query for one parameter. The line starts with a short +description of the parameter, followed by the default parameter value proposed +by ctobssim in squared brackets ``[ ]``. If no parameter is entered (which +is the case for the majority of parameters shown here), the default parameter +will be used. Otherwise, the specified parameter will overwrite the +default parameter. + +You may have recognised that the environment variable ``$GAMMALIB`` has +been used in the path names of the first two parameters. ctools will +automatically expand the environment variables. + +The CTA instrument properties (effective area, PSF width) are taken for +the moment from Konrad Bernloehrs sensitivity files which are stored in +``$GAMMALIB/share/caldb/cta``. + +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 +``$GAMMALIB/share/models/crab.xml``: + +.. code-block:: xml + + + + + + + + + + + + + + + + + + + + + + + + +The model consists of a source library that contains 2 "sources": 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 +3 parameters: the ``Prefactor``, the ``Index`` and the ``Scale``. +The spatial model has 2 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 later). 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`` + +Note that energies are given in MeV. + +The instrumental background of CTA is modelled by a factorized data model +that has a spectral and a radial component (tags ```` and +````, respectively). The spectral component describes the +on-axis background counting rate of CTA as function of energy in units of +counts/s/sr/TeV. The radial component describes the variation of the +background rate with offset angle squared, (i.e. square of the offset angle +with respect to the pointing direction) which is modelled here by a Gaussian. +The only parameter of the radial component is the width of the Gaussian Sigma, +which is here set to 3 degrees squared. + +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 + +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 3 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 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 + + 2014-01-09T20:51:36: +======================+ + 2014-01-09T20:51:36: | Simulate observation | + 2014-01-09T20:51:36: +======================+ + 2014-01-09T20:51:36: === Observation === + 2014-01-09T20:51:36: Simulation area ...........: 1.9635e+11 cm2 + 2014-01-09T20:51:36: Simulation cone ...........: RA=83.63 deg, Dec=22.01 deg, r=10.5 deg + 2014-01-09T20:51:36: Time interval .............: 0 - 1800 s + 2014-01-09T20:51:36: Energy range ..............: 100 GeV - 100 TeV + 2014-01-09T20:51:36: MC source photons .........: 207547 [Crab] + 2014-01-09T20:51:36: MC source events ..........: 995 [Crab] + 2014-01-09T20:51:36: MC source events ..........: 995 (all source models) + 2014-01-09T20:51:36: MC background events ......: 5146 + 2014-01-09T20:51:36: + 2014-01-09T20:51:36: +==================+ + 2014-01-09T20:51:36: | Save observation | + 2014-01-09T20:51:36: +==================+ + 2014-01-09T20:51:37: + 2014-01-09T20:51:37: Application "ctobssim" terminated after 24 wall clock seconds, consuming 0.307116 seconds of CPU time. + +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. 995 of these photons have +been registered by CTA as events. In the same time interval, 5146 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. You may also copy the log file +information to the console using + +.. code-block:: bash + + $ ctobssim debug=yes + +.. _sec_binning_cta: + +Binning CTA data +~~~~~~~~~~~~~~~~ + +As next analysis step you will bin the data in a counts map using the +executable ctbin. A counts map is a 3 dimensional data cube, spanned by +Right Ascension (or Galactic longitude), Declination (or Galactic latitude), +and the logarithm (base 10) of energy. + +ctbin is executed by typing: + +.. code-block:: bash + + $ ctbin + Input event list or observation definition file [test/data/crab_events.fits.gz] events.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] + 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] + Projection method e.g. AIT|AZP|CAR|MER|STG|TAN (AIT|AZP|CAR|MER|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] + Output counts map or observation definition file [cntmap.fits] + +In this example we adjust the event data file name and accept all the +remaining parameter defaults as they perfectly satisfy our needs. The counts +map 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 map 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. + +A set of 20 counts maps will be created, which are logarithmically spaced +in energy, and which cover the energy range from 0.1 TeV to 100 TeV. In this +example, the counts map will be saved as ``cntmap.fits`` in the working +directory. In addition to the counts map, 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: 100% + + *Counts map FITS file* + + +An image of the first bin, covering the energy range 100 - 141 GeV, is +shown below: + +.. figure:: cntmap-map.jpg + :width: 50% + + *Counts map for first energy bin* + + +For illustration, the last few lines of the log file ``ctbin.log`` are +reproduced below: + +.. code-block:: xml + + 2014-01-09T21:16:18: +====================+ + 2014-01-09T21:16:18: | Binned observation | + 2014-01-09T21:16:18: +====================+ + 2014-01-09T21:16:18: === GObservations === + 2014-01-09T21:16:18: Number of observations ....: 1 + 2014-01-09T21:16:18: Number of predicted events : 0 + 2014-01-09T21:16:18: === GCTAObservation === + 2014-01-09T21:16:18: Name ......................: + 2014-01-09T21:16:18: Identifier ................: + 2014-01-09T21:16:18: Instrument ................: CTA + 2014-01-09T21:16:18: Statistics ................: Poisson + 2014-01-09T21:16:18: Ontime ....................: 1800 s + 2014-01-09T21:16:18: Livetime ..................: 1710 s + 2014-01-09T21:16:18: Deadtime correction .......: 0.95 + 2014-01-09T21:16:18: === GCTAPointing === + 2014-01-09T21:16:18: Pointing direction ........: + 2014-01-09T21:16:18: === GCTAResponse === + 2014-01-09T21:16:18: Calibration database ......: + 2014-01-09T21:16:18: Response name .............: + 2014-01-09T21:16:18: RMF file name .............: + 2014-01-09T21:16:18: === GCTAEventCube === + 2014-01-09T21:16:18: Number of events ..........: 5542 + 2014-01-09T21:16:18: Number of elements ........: 800000 + 2014-01-09T21:16:18: Number of pixels ..........: 40000 + 2014-01-09T21:16:18: Number of energy bins .....: 20 + 2014-01-09T21:16:18: Time interval .............: -3.15576e+08 - -3.15574e+08 sec + 2014-01-09T21:16:18: + 2014-01-09T21:16:18: + 2014-01-09T21:16:18: + 2014-01-09T21:16:18: +==================+ + 2014-01-09T21:16:18: | Save observation | + 2014-01-09T21:16:18: +==================+ + 2014-01-09T21:16:18: + 2014-01-09T21:16:18: Application "ctbin" terminated after 37 wall clock seconds, consuming 0.242495 seconds of CPU time. + +From the 6141 events that have been simulated and stored in the +``events.fits`` file, 5542 lie within the map boundaries and are thus put +into the resulting counts map. The counts map 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 ctobssim. Model fitting is done using the executable ctlike, and we +do the fit by typing: + +.. code-block:: bash + + $ ctlike + Event list, counts map or observation definition file [test/data/crab_events.fits.gz] cntmap.fits + Calibration database [test/irf] $GAMMALIB/share/caldb/cta + Instrument response function [cta_dummy_irf] + Source model [test/data/crab.xml] $GAMMALIB/share/models/crab.xml + Source model output file [crab_results.xml] + +Fitting of the data is done in *binned* mode, which means that the events +have been binned into a counts map and the fit computes the log-likelihood +function by summing over all 200 x 200 x 20 bins of the counts map. There is +an alternative method, the so called *unbinned* mode, where the events are +not binned into a counts map and the log-likelihood is computed directly by +summing over all events. We will explore the *unbinned* mode later. + +One of the parameters given to ctlike is a source model output file +(we specified ``crab_results.xml`` in the example), and this file will be +a copy of the model 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 ``Normalization`` and ``Sigma`` parameter of the radial +acceptance model have been fitted (all parameters having the attribute +``free="1"`` are fitted). + +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 + + 2014-01-09T21:26:56: +=================================+ + 2014-01-09T21:26:56: | Maximum likelihood optimisation | + 2014-01-09T21:26:56: +=================================+ + 2014-01-09T21:26:57: Initial iteration: func=22752.5, Lambda=0.001 + 2014-01-09T21:26:59: Iteration 1: func=22750.6, Lambda=0.0001, delta=1.91297, max(grad)=0 [-1] + 2014-01-09T21:27:00: Iteration 2: func=22750.6, Lambda=1e-05, delta=0.00782419, max(grad)=0.000463786 [6] + 2014-01-09T21:27:01: Iteration 3: func=22750.6, Lambda=1e-06, delta=2.49478e-06, max(grad)=0 [-1] + 2014-01-09T21:27:03: Iteration 4: func=22750.6, Lambda=1e-05, delta=0, max(grad)=7.93613e-08 [6] (stalled) + 2014-01-09T21:27:04: Iteration 5: func=22750.6, Lambda=0.0001, delta=0, max(grad)=7.91125e-08 [6] (stalled) + 2014-01-09T21:27:05: Iteration 6: func=22750.6, Lambda=0.001, delta=0, max(grad)=7.66376e-08 [6] (stalled) + 2014-01-09T21:27:07: Iteration 7: func=22750.6, Lambda=0.01, delta=0, max(grad)=1.28271e-07 [7] (stalled) + 2014-01-09T21:27:08: Iteration 8: func=22750.6, Lambda=0.001, delta=2.98678e-09, max(grad)=1.23881e-06 [7] + 2014-01-09T21:27:09: + 2014-01-09T21:27:09: +=========================================+ + 2014-01-09T21:27:09: | Maximum likelihood optimization results | + 2014-01-09T21:27:09: +=========================================+ + 2014-01-09T21:27:09: === GOptimizerLM === + 2014-01-09T21:27:09: Optimized function value ..: 22750.6 + 2014-01-09T21:27:09: Absolute precision ........: 1e-06 + 2014-01-09T21:27:09: Optimization status .......: converged + 2014-01-09T21:27:09: Number of parameters ......: 9 + 2014-01-09T21:27:09: Number of free parameters .: 4 + 2014-01-09T21:27:09: Number of iterations ......: 8 + 2014-01-09T21:27:09: Lambda ....................: 0.001 + 2014-01-09T21:27:09: Maximum log likelihood ....: -22750.6 + 2014-01-09T21:27:09: Observed events (Nobs) ...: 5542 + 2014-01-09T21:27:09: Predicted events (Npred) ..: 5542 (Nobs - Npred = 5.52808e-06) + 2014-01-09T21:27:09: === GModels === + 2014-01-09T21:27:09: Number of models ..........: 2 + 2014-01-09T21:27:09: Number of parameters ......: 9 + 2014-01-09T21:27:09: === GModelSky === + 2014-01-09T21:27:09: Name ......................: Crab + 2014-01-09T21:27:09: Instruments ...............: all + 2014-01-09T21:27:09: Instrument scale factors ..: unity + 2014-01-09T21:27:09: Observation identifiers ...: all + 2014-01-09T21:27:09: Model type ................: PointSource + 2014-01-09T21:27:09: Model components ..........: "SkyDirFunction" * "PowerLaw" * "Constant" + 2014-01-09T21:27:09: Number of parameters ......: 6 + 2014-01-09T21:27:09: Number of spatial par's ...: 2 + 2014-01-09T21:27:09: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1) + 2014-01-09T21:27:09: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1) + 2014-01-09T21:27:09: Number of spectral par's ..: 3 + 2014-01-09T21:27:09: Prefactor ................: 6.07927e-16 +/- 2.04581e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) + 2014-01-09T21:27:09: Index ....................: -2.5009 +/- 0.0252056 [-0,-5] (free,scale=-1,gradient) + 2014-01-09T21:27:09: PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) + 2014-01-09T21:27:09: Number of temporal par's ..: 1 + 2014-01-09T21:27:09: Constant .................: 1 (relative value) (fixed,scale=1,gradient) + 2014-01-09T21:27:09: === GCTAModelRadialAcceptance === + 2014-01-09T21:27:09: Name ......................: Background + 2014-01-09T21:27:09: Instruments ...............: CTA + 2014-01-09T21:27:09: Instrument scale factors ..: unity + 2014-01-09T21:27:09: Observation identifiers ...: all + 2014-01-09T21:27:09: Model type ................: "Gaussian" * "FileFunction" * "Constant" + 2014-01-09T21:27:09: Number of parameters ......: 3 + 2014-01-09T21:27:09: Number of radial par's ....: 1 + 2014-01-09T21:27:09: Sigma ....................: 3.02527 +/- 0.0594142 [0.01,10] deg2 (free,scale=1,gradient) + 2014-01-09T21:27:09: Number of spectral par's ..: 1 + 2014-01-09T21:27:09: Normalization ............: 0.99069 +/- 0.0192887 [0,1000] (free,scale=1,gradient) + 2014-01-09T21:27:09: Number of temporal par's ..: 1 + 2014-01-09T21:27:09: Constant .................: 1 (relative value) (fixed,scale=1,gradient) + 2014-01-09T21:27:09: + 2014-01-09T21:27:09: +==============+ + 2014-01-09T21:27:09: | Save results | + 2014-01-09T21:27:09: +==============+ + 2014-01-09T21:27:09: + 2014-01-09T21:27:09: Application "ctlike" terminated after 55 wall clock seconds, consuming 13.4985 seconds of CPU time. + +The maximum likelihood optimizer required 8 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 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 of the best fitting values). + + +.. _sec_unbinned_cta: + +Doing an unbinned analysis +~~~~~~~~~~~~~~~~~~~~~~~~~~ + +As gamma-ray events are rare, the counts maps generated by ctbin will in +general be sparse, 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 map. We will see the benefit of +such an analysis later once you re-run 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 ctbin step in binned +analysis where you defined the size of the counts map, the energy range, and +the time interval. For unbinned analysis you have no such thing as a counts +map, 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 executable +ctselect, which replaces the ctbin step in an unbinned analysis. + +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 [test/data/crab_events.fits.gz] 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.0] + End time (CTA MET in seconds) (0) [0.0] + Lower energy limit (TeV) (0) [0.1] + Upper energy limit (TeV) (0) [100.0] + Output event list or observation definition file [selected_events.fits] + +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: 60% + + *Data selection keywords* + +It is important for unbinned analysis that these data selection keywords +exist. If they don't exist, ctlike will not execute in unbinned mode. + +Below some lines of the ``ctselect.log`` file that show the data selection +part: + +.. code-block:: xml + + 2014-01-09T21:36:20: +=================+ + 2014-01-09T21:36:20: | Event selection | + 2014-01-09T21:36:20: +=================+ + 2014-01-09T21:36:20: === Observation === + 2014-01-09T21:36:21: Energy range ..............: 0.1 - 100 TeV + 2014-01-09T21:36:21: Acceptance cone centre ....: RA=83.63, DEC=22.01 deg + 2014-01-09T21:36:21: Acceptance cone radius ....: 3 deg + 2014-01-09T21:36:21: cfitsio selection .........: ENERGY >= 0.10000000 && ENERGY <= 100.00000000 && ANGSEP(83.630000,22.010000,RA,DEC) <= 3.000000 + 2014-01-09T21:36:21: FITS filename .............: ctselectKKHEiF[EVENTS][ENERGY >= 0.10000000 && ENERGY <= 100.00000000 && ANGSEP(83.630000,22.010000,RA,DEC) <= 3.000000] + +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 map: + +.. code-block:: bash + + $ ctlike + Event list, counts map or observation definition file [cntmap.fits] selected_events.fits + Calibration database [$GAMMALIB/share/caldb/cta] + Instrument response function [cta_dummy_irf] + Source model [$GAMMALIB/share/models/crab.xml] + Source model output file [crab_results.xml] + +You will recognise that ctlike runs much faster in unbinned mode compared +to binned mode. This is understandable as the selected event list contains +only 6127 events, while the binned counts map 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 the pixels). So as long +as you work with short event lists, unbinned mode is faster. Unbinned 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 maps, 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 + + 2014-01-09T21:47:54: +=================================+ + 2014-01-09T21:47:54: | Maximum likelihood optimisation | + 2014-01-09T21:47:54: +=================================+ + 2014-01-09T21:47:54: Initial iteration: func=44393.9, Lambda=0.001 + 2014-01-09T21:47:54: Iteration 1: func=44392.2, Lambda=0.001, delta=1.7034, max(grad)=9.60142 [6] + 2014-01-09T21:47:54: Iteration 2: func=44392.1, Lambda=0.0001, delta=0.099482, max(grad)=0 [-1] + 2014-01-09T21:47:54: Iteration 3: func=44392.1, Lambda=1e-05, delta=6.63429e-05, max(grad)=0.00559314 [6] + 2014-01-09T21:47:54: Iteration 4: func=44392.1, Lambda=1e-06, delta=4.7803e-08, max(grad)=0 [-1] + 2014-01-09T21:47:54: + 2014-01-09T21:47:54: +=========================================+ + 2014-01-09T21:47:54: | Maximum likelihood optimization results | + 2014-01-09T21:47:54: +=========================================+ + 2014-01-09T21:47:54: === GOptimizerLM === + 2014-01-09T21:47:54: Optimized function value ..: 44392.1 + 2014-01-09T21:47:54: Absolute precision ........: 1e-06 + 2014-01-09T21:47:54: Optimization status .......: converged + 2014-01-09T21:47:54: Number of parameters ......: 9 + 2014-01-09T21:47:54: Number of free parameters .: 4 + 2014-01-09T21:47:54: Number of iterations ......: 4 + 2014-01-09T21:47:54: Lambda ....................: 1e-07 + 2014-01-09T21:47:54: Maximum log likelihood ....: -44392.1 + 2014-01-09T21:47:54: Observed events (Nobs) ...: 6127 + 2014-01-09T21:47:54: Predicted events (Npred) ..: 6127 (Nobs - Npred = 1.01613e-07) + 2014-01-09T21:47:54: === GModels === + 2014-01-09T21:47:54: Number of models ..........: 2 + 2014-01-09T21:47:54: Number of parameters ......: 9 + 2014-01-09T21:47:54: === GModelSky === + 2014-01-09T21:47:54: Name ......................: Crab + 2014-01-09T21:47:54: Instruments ...............: all + 2014-01-09T21:47:54: Instrument scale factors ..: unity + 2014-01-09T21:47:54: Observation identifiers ...: all + 2014-01-09T21:47:54: Model type ................: PointSource + 2014-01-09T21:47:54: Model components ..........: "SkyDirFunction" * "PowerLaw" * "Constant" + 2014-01-09T21:47:54: Number of parameters ......: 6 + 2014-01-09T21:47:54: Number of spatial par's ...: 2 + 2014-01-09T21:47:54: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1) + 2014-01-09T21:47:54: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1) + 2014-01-09T21:47:54: Number of spectral par's ..: 3 + 2014-01-09T21:47:54: Prefactor ................: 6.13416e-16 +/- 2.05772e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) + 2014-01-09T21:47:54: Index ....................: -2.50585 +/- 0.0250821 [-0,-5] (free,scale=-1,gradient) + 2014-01-09T21:47:54: PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) + 2014-01-09T21:47:54: Number of temporal par's ..: 1 + 2014-01-09T21:47:54: Constant .................: 1 (relative value) (fixed,scale=1,gradient) + 2014-01-09T21:47:54: === GCTAModelRadialAcceptance === + 2014-01-09T21:47:54: Name ......................: Background + 2014-01-09T21:47:54: Instruments ...............: CTA + 2014-01-09T21:47:54: Instrument scale factors ..: unity + 2014-01-09T21:47:54: Observation identifiers ...: all + 2014-01-09T21:47:54: Model type ................: "Gaussian" * "FileFunction" * "Constant" + 2014-01-09T21:47:54: Number of parameters ......: 3 + 2014-01-09T21:47:54: Number of radial par's ....: 1 + 2014-01-09T21:47:54: Sigma ....................: 3.04429 +/- 0.0330229 [0.01,10] deg2 (free,scale=1,gradient) + 2014-01-09T21:47:54: Number of spectral par's ..: 1 + 2014-01-09T21:47:54: Normalization ............: 0.997516 +/- 0.0175361 [0,1000] (free,scale=1,gradient) + 2014-01-09T21:47:54: Number of temporal par's ..: 1 + 2014-01-09T21:47:54: Constant .................: 1 (relative value) (fixed,scale=1,gradient) + 2014-01-09T21:47:54: + 2014-01-09T21:47:54: +==============+ + 2014-01-09T21:47:54: | Save results | + 2014-01-09T21:47:54: +==============+ + 2014-01-09T21:47:54: + 2014-01-09T21:47:54: Application "ctlike" terminated after 6 wall clock seconds, consuming 0.163384 seconds of CPU time. + + +Models available in ctools +-------------------------- + +Using ctools from Python +------------------------ + +Frequently asked questions +-------------------------- + -- libgit2 0.21.4