Johan Bregeon / ctools

		1	+Getting started
		2	+===============
		3	+
		4	+A quick ctools tutorial
		5	+-----------------------
		6	+
		7	+First steps with ctools
		8	+~~~~~~~~~~~~~~~~~~~~~~~
		9	+
		10	+This section provides a step-by-step description of the CTA data analysis
		11	+procedure implemented in ctools. In this example we will simulate data for
		12	+the Crab nebula, observed on-axis, as a single source in the CTA field of
		13	+view.
		14	+
		15	+The example consists of the following analysis steps:
		16	+
		17	+- :ref:`sec_simulating_cta`
		18	+
		19	+- :ref:`sec_binning_cta`
		20	+
		21	+- :ref:`sec_fitting_cta`
		22	+
		23	+- :ref:`sec_unbinned_cta`
		24	+
		25	+
		26	+.. _sec_simulating_cta:
		27	+
		28	+Simulating CTA data
		29	+~~~~~~~~~~~~~~~~~~~
		30	+
		31	+CTA data are simulated by the executable ctobssim. To start the executable,
		32	+type ctobssim at the console prompt (which is denoted by ``$``). ctobssim
		33	+will query for a number of parameters:
		34	+
		35	+.. code-block:: bash
		36	+
		37	+ $ ctobssim
		38	+ Model [test/data/crab.xml] $GAMMALIB/share/models/crab.xml
		39	+ Calibration database [test/irf] $GAMMALIB/share/caldb/cta
		40	+ Instrument response function [cta_dummy_irf]
		41	+ RA of pointing (degrees) (0-360) [83.63]
		42	+ Dec of pointing (degrees) (-90-90) [22.01]
		43	+ Radius of FOV (degrees) (0-180) [10.0]
		44	+ Start time (MET in s) (0) [0.0]
		45	+ End time (MET in s) (0) [1800.0]
		46	+ Lower energy limit (TeV) (0) [0.1]
		47	+ Upper energy limit (TeV) (0) [100.0]
		48	+ Output event data file or observation definition file [events.fits]
		49	+
		50	+Each line represents a query for one parameter. The line starts with a short
		51	+description of the parameter, followed by the default parameter value proposed
		52	+by ctobssim in squared brackets ``[ ]``. If no parameter is entered (which
		53	+is the case for the majority of parameters shown here), the default parameter
		54	+will be used. Otherwise, the specified parameter will overwrite the
		55	+default parameter.
		56	+
		57	+You may have recognised that the environment variable ``$GAMMALIB`` has
		58	+been used in the path names of the first two parameters. ctools will
		59	+automatically expand the environment variables.
		60	+
		61	+The CTA instrument properties (effective area, PSF width) are taken for
		62	+the moment from Konrad Bernloehrs sensitivity files which are stored in
		63	+``$GAMMALIB/share/caldb/cta``.
		64	+
		65	+Events are simulated based on the instrument properties and based on a
		66	+source and background model. Only events that fall within the specified
		67	+region of interest (ROI), defined as a circle around a sky position in
		68	+Right Ascension and Declination (in degrees), will be stored in the output
		69	+event data file. The duration of the simulation is taken here to 30 minutes
		70	+(or 1800 seconds). Events are simulated for energies between 0.1 and 100 TeV.
		71	+
		72	+The source and background model is defined by the XML file
		73	+``$GAMMALIB/share/models/crab.xml``:
		74	+
		75	+.. code-block:: xml
		76	+
		77	+ <?xml version="1.0" standalone="no"?>
		78	+ <source_library title="source library">
		79	+ <source name="MAGIC Crab (Albert et al. 2008, ApJ, 674, 1037)" type="PointSource">
		80	+ <spectrum type="PowerLaw">
		81	+ <parameter name="Prefactor" scale="1e-16" value="5.7" min="1e-07" max="1000.0" free="1"/>
		82	+ <parameter name="Index" scale="-1" value="2.48" min="0.0" max="+5.0" free="1"/>
		83	+ <parameter name="Scale" scale="1e6" value="0.3" min="0.01" max="1000.0" free="0"/>
		84	+ </spectrum>
		85	+ <spatialModel type="SkyDirFunction">
		86	+ <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="0"/>
		87	+ <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="0"/>
		88	+ </spatialModel>
		89	+ </source>
		90	+ <source name="Background" type="RadialAcceptance" instrument="CTA">
		91	+ <spectrum type="FileFunction" file="$GAMMALIB/share/models/bkg_dummy.txt">
		92	+ <parameter scale="1.0" name="Normalization" min="0.0" max="1000.0" value="1.0" free="1"/>
		93	+ </spectrum>
		94	+ <radialModel type="Gaussian">
		95	+ <parameter name="Sigma" scale="1.0" value="3.0" min="0.01" max="10.0" free="1"/>
		96	+ </radialModel>
		97	+ </source>
		98	+ </source_library>
		99	+
		100	+The model consists of a source library that contains 2 "sources": the Crab
		101	+nebula and an instrumental background model.
		102	+
		103	+The Crab nebula is modelled by a factorized sky model that has a spectral
		104	+and a spatial component (tags ``<spectrum>`` and ``<spatialModel>``,
		105	+respectively). The spectrum is modelled by a power law, which is defined by
		106	+3 parameters: the ``Prefactor``, the ``Index`` and the ``Scale``.
		107	+The spatial model has 2 parameters: Right Ascension in degrees (RA), and
		108	+Declination in degrees (DEC). Each parameter has a value and a scale factor,
		109	+the real value of the parameter being the product value * scale. Typically,
		110	+scale is chosen so that value is of the order of 1 (this is relevant for
		111	+model fitting later). In addition, value is bound by a minimum (min) and
		112	+maximum (max) value, and a parameter may be free (``free="1"``) or fixed
		113	+(``free="0"``). The min, max, and free attributes are not relevant here for
		114	+the simulations, but they will be important for the model fitting later.
		115	+
		116	+The spectral intensity I(E) (in units of photons/cm2/s/MeV) of the power
		117	+law is given by
		118	+
		119	+.. math::
		120	+ \frac{dN}{dE} = N_0 \left( \frac{E}{E_0} \right)^{\gamma}
		121	+
		122	+where the parameters in the XML definition have the following mappings:
		123	+
		124	+* :math:`N_0` = ``Prefactor``
		125	+* :math:`\gamma` = ``Index``
		126	+* :math:`E_0` = ``Scale``
		127	+
		128	+Note that energies are given in MeV.
		129	+
		130	+The instrumental background of CTA is modelled by a factorized data model
		131	+that has a spectral and a radial component (tags ``<spectrum>`` and
		132	+``<radialModel>``, respectively). The spectral component describes the
		133	+on-axis background counting rate of CTA as function of energy in units of
		134	+counts/s/sr/TeV. The radial component describes the variation of the
		135	+background rate with offset angle squared, (i.e. square of the offset angle
		136	+with respect to the pointing direction) which is modelled here by a Gaussian.
		137	+The only parameter of the radial component is the width of the Gaussian Sigma,
		138	+which is here set to 3 degrees squared.
		139	+
		140	+ctobssim has a couple of hidden parameters, the most important one being
		141	+certainly ``seed``. ``seed`` is an integer that specifies the seed value
		142	+for the random number generator, and changing this parameter will allow to
		143	+generate statistically independent Monte Carlo samples of CTA event data.
		144	+To use for example a seed value of 41 you should type:
		145	+
		146	+.. code-block:: bash
		147	+
		148	+ $ ctobssim seed=41
		149	+
		150	+ctobssim will write 2 files in the working directory: ``events.fits`` and
		151	+``ctobssim.log``. The first file contains the simulated events in FITS
		152	+format and can be inspected using ``fv`` or ``ds9``. The FITS file will
		153	+contain 3 extensions: an empty primary image, a binary table named
		154	+``EVENTS`` that holds the events (one row per event), and a binary table
		155	+named ``GTI`` holding the Good Time Intervals (for the moment a single row
		156	+with 2 columns providing the start and the stop time of the simulated time
		157	+interval).
		158	+
		159	+The second file produced by ctobssim is a human readable log file that
		160	+contains information about the job execution. As example, the last lines
		161	+from this file are shown here:
		162	+
		163	+.. code-block:: xml
		164	+
		165	+ 2014-01-09T20:51:36: +======================+
		166	+ 2014-01-09T20:51:36: \| Simulate observation \|
		167	+ 2014-01-09T20:51:36: +======================+
		168	+ 2014-01-09T20:51:36: === Observation ===
		169	+ 2014-01-09T20:51:36: Simulation area ...........: 1.9635e+11 cm2
		170	+ 2014-01-09T20:51:36: Simulation cone ...........: RA=83.63 deg, Dec=22.01 deg, r=10.5 deg
		171	+ 2014-01-09T20:51:36: Time interval .............: 0 - 1800 s
		172	+ 2014-01-09T20:51:36: Energy range ..............: 100 GeV - 100 TeV
		173	+ 2014-01-09T20:51:36: MC source photons .........: 207547 [Crab]
		174	+ 2014-01-09T20:51:36: MC source events ..........: 995 [Crab]
		175	+ 2014-01-09T20:51:36: MC source events ..........: 995 (all source models)
		176	+ 2014-01-09T20:51:36: MC background events ......: 5146
		177	+ 2014-01-09T20:51:36:
		178	+ 2014-01-09T20:51:36: +==================+
		179	+ 2014-01-09T20:51:36: \| Save observation \|
		180	+ 2014-01-09T20:51:36: +==================+
		181	+ 2014-01-09T20:51:37:
		182	+ 2014-01-09T20:51:37: Application "ctobssim" terminated after 24 wall clock seconds, consuming 0.307116 seconds of CPU time.
		183	+
		184	+Each line starts with the UTC time at which the line has been written. In
		185	+this run, 207547 Crab photons have been thrown over an area of 19.6 square
		186	+kilometres during a time interval of 1800 seconds. 995 of these photons have
		187	+been registered by CTA as events. In the same time interval, 5146 background
		188	+events have been registred by CTA.
		189	+
		190	+You may change the name of the log file using the hidden parameter
		191	+``logfile``:
		192	+
		193	+.. code-block:: bash
		194	+
		195	+ $ ctobssim logfile=my-private-log-file
		196	+
		197	+Furthermore, you may decide on the amount of information provided in the
		198	+log file (the chattiness of the executable) using the hidden parameter
		199	+``chatter``:
		200	+
		201	+.. code-block:: bash
		202	+
		203	+ $ ctobssim chatter=4
		204	+
		205	+``chatter`` can vary between 0 and 4, 0 providing no information while 4
		206	+provides the most detailed information. You may also copy the log file
		207	+information to the console using
		208	+
		209	+.. code-block:: bash
		210	+
		211	+ $ ctobssim debug=yes
		212	+
		213	+.. _sec_binning_cta:
		214	+
		215	+Binning CTA data
		216	+~~~~~~~~~~~~~~~~
		217	+
		218	+As next analysis step you will bin the data in a counts map using the
		219	+executable ctbin. A counts map is a 3 dimensional data cube, spanned by
		220	+Right Ascension (or Galactic longitude), Declination (or Galactic latitude),
		221	+and the logarithm (base 10) of energy.
		222	+
		223	+ctbin is executed by typing:
		224	+
		225	+.. code-block:: bash
		226	+
		227	+ $ ctbin
		228	+ Input event list or observation definition file [test/data/crab_events.fits.gz] events.fits
		229	+ First coordinate of image center in degrees (RA or galactic l) [83.63]
		230	+ Second coordinate of image center in degrees (DEC or galactic b) [22.01]
		231	+ Algorithm for defining energy bins (FILE\|LIN\|LOG) [LOG]
		232	+ Start value for first energy bin in TeV [0.1]
		233	+ Stop value for last energy bin in TeV [100.0]
		234	+ Number of energy bins [20]
		235	+ Projection method e.g. AIT\|AZP\|CAR\|MER\|STG\|TAN (AIT\|AZP\|CAR\|MER\|STG\|TAN) [CAR]
		236	+ Coordinate system (CEL - celestial, GAL - galactic) (CEL\|GAL) [CEL]
		237	+ Image scale (in degrees/pixel) [0.02]
		238	+ Size of the X axis in pixels [200]
		239	+ Size of the Y axis in pixels [200]
		240	+ Output counts map or observation definition file [cntmap.fits]
		241	+
		242	+In this example we adjust the event data file name and accept all the
		243	+remaining parameter defaults as they perfectly satisfy our needs. The counts
		244	+map will be centred on the location of the Crab (Right Ascension 83.63 degrees,
		245	+Declination 22.01 degrees) and will be aligned in celestial coordinates. A
		246	+cartesian projection has been selected. The counts map has 200 x 200 spatial
		247	+pixels of 0.02 x 0.02 degrees in size, hence it covers a total area of 4 x 4
		248	+degrees.
		249	+
		250	+A set of 20 counts maps will be created, which are logarithmically spaced
		251	+in energy, and which cover the energy range from 0.1 TeV to 100 TeV. In this
		252	+example, the counts map will be saved as ``cntmap.fits`` in the working
		253	+directory. In addition to the counts map, that is stored as the primary
		254	+image extension, the FITS file also contains an extension named ``EBOUNDS``
		255	+that defines the energy boundaries that were used, and an extension ``GTI``
		256	+that defines the Good Time Intervals that have been used. The following
		257	+image shows the resulting FITS file. The ``EBOUNDS`` table has 20 rows, one
		258	+for each energy bin, while the ``GTI`` table has just a single row, indicating
		259	+the start and stop time of the simulated data.
		260	+
		261	+.. figure:: cntmap-fits.jpg
		262	+ :width: 100%
		263	+
		264	+ Counts map FITS file
		265	+
		266	+
		267	+An image of the first bin, covering the energy range 100 - 141 GeV, is
		268	+shown below:
		269	+
		270	+.. figure:: cntmap-map.jpg
		271	+ :width: 50%
		272	+
		273	+ Counts map for first energy bin
		274	+
		275	+
		276	+For illustration, the last few lines of the log file ``ctbin.log`` are
		277	+reproduced below:
		278	+
		279	+.. code-block:: xml
		280	+
		281	+ 2014-01-09T21:16:18: +====================+
		282	+ 2014-01-09T21:16:18: \| Binned observation \|
		283	+ 2014-01-09T21:16:18: +====================+
		284	+ 2014-01-09T21:16:18: === GObservations ===
		285	+ 2014-01-09T21:16:18: Number of observations ....: 1
		286	+ 2014-01-09T21:16:18: Number of predicted events : 0
		287	+ 2014-01-09T21:16:18: === GCTAObservation ===
		288	+ 2014-01-09T21:16:18: Name ......................:
		289	+ 2014-01-09T21:16:18: Identifier ................:
		290	+ 2014-01-09T21:16:18: Instrument ................: CTA
		291	+ 2014-01-09T21:16:18: Statistics ................: Poisson
		292	+ 2014-01-09T21:16:18: Ontime ....................: 1800 s
		293	+ 2014-01-09T21:16:18: Livetime ..................: 1710 s
		294	+ 2014-01-09T21:16:18: Deadtime correction .......: 0.95
		295	+ 2014-01-09T21:16:18: === GCTAPointing ===
		296	+ 2014-01-09T21:16:18: Pointing direction ........:
		297	+ 2014-01-09T21:16:18: === GCTAResponse ===
		298	+ 2014-01-09T21:16:18: Calibration database ......:
		299	+ 2014-01-09T21:16:18: Response name .............:
		300	+ 2014-01-09T21:16:18: RMF file name .............:
		301	+ 2014-01-09T21:16:18: === GCTAEventCube ===
		302	+ 2014-01-09T21:16:18: Number of events ..........: 5542
		303	+ 2014-01-09T21:16:18: Number of elements ........: 800000
		304	+ 2014-01-09T21:16:18: Number of pixels ..........: 40000
		305	+ 2014-01-09T21:16:18: Number of energy bins .....: 20
		306	+ 2014-01-09T21:16:18: Time interval .............: -3.15576e+08 - -3.15574e+08 sec
		307	+ 2014-01-09T21:16:18:
		308	+ 2014-01-09T21:16:18:
		309	+ 2014-01-09T21:16:18:
		310	+ 2014-01-09T21:16:18: +==================+
		311	+ 2014-01-09T21:16:18: \| Save observation \|
		312	+ 2014-01-09T21:16:18: +==================+
		313	+ 2014-01-09T21:16:18:
		314	+ 2014-01-09T21:16:18: Application "ctbin" terminated after 37 wall clock seconds, consuming 0.242495 seconds of CPU time.
		315	+
		316	+From the 6141 events that have been simulated and stored in the
		317	+``events.fits`` file, 5542 lie within the map boundaries and are thus put
		318	+into the resulting counts map. The counts map is stored in a cartesian
		319	+projection in a World Coordinate System (WCS) compliant format.
		320	+
		321	+
		322	+.. _sec_fitting_cta:
		323	+
		324	+Fitting CTA data
		325	+~~~~~~~~~~~~~~~~
		326	+
		327	+Now we are ready to fit the simulated data with a model. For simplicity
		328	+we use in this example the same model that we used to simulate the data
		329	+with ctobssim. Model fitting is done using the executable ctlike, and we
		330	+do the fit by typing:
		331	+
		332	+.. code-block:: bash
		333	+
		334	+ $ ctlike
		335	+ Event list, counts map or observation definition file [test/data/crab_events.fits.gz] cntmap.fits
		336	+ Calibration database [test/irf] $GAMMALIB/share/caldb/cta
		337	+ Instrument response function [cta_dummy_irf]
		338	+ Source model [test/data/crab.xml] $GAMMALIB/share/models/crab.xml
		339	+ Source model output file [crab_results.xml]
		340	+
		341	+Fitting of the data is done in binned mode, which means that the events
		342	+have been binned into a counts map and the fit computes the log-likelihood
		343	+function by summing over all 200 x 200 x 20 bins of the counts map. There is
		344	+an alternative method, the so called unbinned mode, where the events are
		345	+not binned into a counts map and the log-likelihood is computed directly by
		346	+summing over all events. We will explore the unbinned mode later.
		347	+
		348	+One of the parameters given to ctlike is a source model output file
		349	+(we specified ``crab_results.xml`` in the example), and this file will be
		350	+a copy of the model XML file where the parameter values have been replaced
		351	+by the fit results. In addition, the statistical uncertainties are added
		352	+for each fitted parameter using the attribute error. Below we show the XML
		353	+result file that has been produced by the run:
		354	+
		355	+.. code-block:: xml
		356	+
		357	+ <?xml version="1.0" encoding="UTF-8" standalone="no"?>
		358	+ <source_library title="source library">
		359	+ <source name="Crab" type="PointSource">
		360	+ <spectrum type="PowerLaw">
		361	+ <parameter name="Prefactor" value="6.07927" error="0.204581" scale="1e-16" min="1e-07" max="1000" free="1" />
		362	+ <parameter name="Index" value="2.5009" error="0.0252056" scale="-1" min="0" max="5" free="1" />
		363	+ <parameter name="Scale" value="0.3" scale="1e+06" min="0.01" max="1000" free="0" />
		364	+ </spectrum>
		365	+ <spatialModel type="SkyDirFunction">
		366	+ <parameter name="RA" value="83.6331" scale="1" min="-360" max="360" free="0" />
		367	+ <parameter name="DEC" value="22.0145" scale="1" min="-90" max="90" free="0" />
		368	+ </spatialModel>
		369	+ </source>
		370	+ <source name="Background" type="RadialAcceptance" instrument="CTA">
		371	+ <spectrum type="FileFunction" file="$GAMMALIB/share/models/bkg_dummy.txt">
		372	+ <parameter name="Normalization" value="0.99069" error="0.0192887" scale="1" min="0" max="1000" free="1" />
		373	+ </spectrum>
		374	+ <radialModel type="Gaussian">
		375	+ <parameter name="Sigma" value="3.02527" error="0.0594142" scale="1" min="0.01" max="10" free="1" />
		376	+ </radialModel>
		377	+ </source>
		378	+ </source_library>
		379	+
		380	+In this example, the ``Prefactor`` and ``Index`` of the spectral model for the
		381	+Crab as well as the ``Normalization`` and ``Sigma`` parameter of the radial
		382	+acceptance model have been fitted (all parameters having the attribute
		383	+``free="1"`` are fitted).
		384	+
		385	+To get more details about the model fitting you can inspect the log file.
		386	+Below the last lines of the ctlike.log log file that has been produced by
		387	+this run:
		388	+
		389	+.. code-block:: xml
		390	+
		391	+ 2014-01-09T21:26:56: +=================================+
		392	+ 2014-01-09T21:26:56: \| Maximum likelihood optimisation \|
		393	+ 2014-01-09T21:26:56: +=================================+
		394	+ 2014-01-09T21:26:57: Initial iteration: func=22752.5, Lambda=0.001
		395	+ 2014-01-09T21:26:59: Iteration 1: func=22750.6, Lambda=0.0001, delta=1.91297, max(grad)=0 [-1]
		396	+ 2014-01-09T21:27:00: Iteration 2: func=22750.6, Lambda=1e-05, delta=0.00782419, max(grad)=0.000463786 [6]
		397	+ 2014-01-09T21:27:01: Iteration 3: func=22750.6, Lambda=1e-06, delta=2.49478e-06, max(grad)=0 [-1]
		398	+ 2014-01-09T21:27:03: Iteration 4: func=22750.6, Lambda=1e-05, delta=0, max(grad)=7.93613e-08 [6] (stalled)
		399	+ 2014-01-09T21:27:04: Iteration 5: func=22750.6, Lambda=0.0001, delta=0, max(grad)=7.91125e-08 [6] (stalled)
		400	+ 2014-01-09T21:27:05: Iteration 6: func=22750.6, Lambda=0.001, delta=0, max(grad)=7.66376e-08 [6] (stalled)
		401	+ 2014-01-09T21:27:07: Iteration 7: func=22750.6, Lambda=0.01, delta=0, max(grad)=1.28271e-07 [7] (stalled)
		402	+ 2014-01-09T21:27:08: Iteration 8: func=22750.6, Lambda=0.001, delta=2.98678e-09, max(grad)=1.23881e-06 [7]
		403	+ 2014-01-09T21:27:09:
		404	+ 2014-01-09T21:27:09: +=========================================+
		405	+ 2014-01-09T21:27:09: \| Maximum likelihood optimization results \|
		406	+ 2014-01-09T21:27:09: +=========================================+
		407	+ 2014-01-09T21:27:09: === GOptimizerLM ===
		408	+ 2014-01-09T21:27:09: Optimized function value ..: 22750.6
		409	+ 2014-01-09T21:27:09: Absolute precision ........: 1e-06
		410	+ 2014-01-09T21:27:09: Optimization status .......: converged
		411	+ 2014-01-09T21:27:09: Number of parameters ......: 9
		412	+ 2014-01-09T21:27:09: Number of free parameters .: 4
		413	+ 2014-01-09T21:27:09: Number of iterations ......: 8
		414	+ 2014-01-09T21:27:09: Lambda ....................: 0.001
		415	+ 2014-01-09T21:27:09: Maximum log likelihood ....: -22750.6
		416	+ 2014-01-09T21:27:09: Observed events (Nobs) ...: 5542
		417	+ 2014-01-09T21:27:09: Predicted events (Npred) ..: 5542 (Nobs - Npred = 5.52808e-06)
		418	+ 2014-01-09T21:27:09: === GModels ===
		419	+ 2014-01-09T21:27:09: Number of models ..........: 2
		420	+ 2014-01-09T21:27:09: Number of parameters ......: 9
		421	+ 2014-01-09T21:27:09: === GModelSky ===
		422	+ 2014-01-09T21:27:09: Name ......................: Crab
		423	+ 2014-01-09T21:27:09: Instruments ...............: all
		424	+ 2014-01-09T21:27:09: Instrument scale factors ..: unity
		425	+ 2014-01-09T21:27:09: Observation identifiers ...: all
		426	+ 2014-01-09T21:27:09: Model type ................: PointSource
		427	+ 2014-01-09T21:27:09: Model components ..........: "SkyDirFunction" * "PowerLaw" * "Constant"
		428	+ 2014-01-09T21:27:09: Number of parameters ......: 6
		429	+ 2014-01-09T21:27:09: Number of spatial par's ...: 2
		430	+ 2014-01-09T21:27:09: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1)
		431	+ 2014-01-09T21:27:09: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1)
		432	+ 2014-01-09T21:27:09: Number of spectral par's ..: 3
		433	+ 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)
		434	+ 2014-01-09T21:27:09: Index ....................: -2.5009 +/- 0.0252056 [-0,-5] (free,scale=-1,gradient)
		435	+ 2014-01-09T21:27:09: PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient)
		436	+ 2014-01-09T21:27:09: Number of temporal par's ..: 1
		437	+ 2014-01-09T21:27:09: Constant .................: 1 (relative value) (fixed,scale=1,gradient)
		438	+ 2014-01-09T21:27:09: === GCTAModelRadialAcceptance ===
		439	+ 2014-01-09T21:27:09: Name ......................: Background
		440	+ 2014-01-09T21:27:09: Instruments ...............: CTA
		441	+ 2014-01-09T21:27:09: Instrument scale factors ..: unity
		442	+ 2014-01-09T21:27:09: Observation identifiers ...: all
		443	+ 2014-01-09T21:27:09: Model type ................: "Gaussian" * "FileFunction" * "Constant"
		444	+ 2014-01-09T21:27:09: Number of parameters ......: 3
		445	+ 2014-01-09T21:27:09: Number of radial par's ....: 1
		446	+ 2014-01-09T21:27:09: Sigma ....................: 3.02527 +/- 0.0594142 [0.01,10] deg2 (free,scale=1,gradient)
		447	+ 2014-01-09T21:27:09: Number of spectral par's ..: 1
		448	+ 2014-01-09T21:27:09: Normalization ............: 0.99069 +/- 0.0192887 [0,1000] (free,scale=1,gradient)
		449	+ 2014-01-09T21:27:09: Number of temporal par's ..: 1
		450	+ 2014-01-09T21:27:09: Constant .................: 1 (relative value) (fixed,scale=1,gradient)
		451	+ 2014-01-09T21:27:09:
		452	+ 2014-01-09T21:27:09: +==============+
		453	+ 2014-01-09T21:27:09: \| Save results \|
		454	+ 2014-01-09T21:27:09: +==============+
		455	+ 2014-01-09T21:27:09:
		456	+ 2014-01-09T21:27:09: Application "ctlike" terminated after 55 wall clock seconds, consuming 13.4985 seconds of CPU time.
		457	+
		458	+The maximum likelihood optimizer required 8 iterations to converge. This
		459	+is pretty fast, but recall that we used the same model file for the simulation
		460	+and for fitting, hence the initial parameter values were already very close
		461	+to the best fitting values. To see the impact of the initial parameters on
		462	+the fit result, you may re-run ctlike using another copy of the model XML file
		463	+where you change the value attributes of the parameters that should be
		464	+fitted. You will see that the optimizer requires a couple of more iterations,
		465	+but it should converge to the same solution (provided that the initial values
		466	+are not too far of the best fitting values).
		467	+
		468	+
		469	+.. _sec_unbinned_cta:
		470	+
		471	+Doing an unbinned analysis
		472	+~~~~~~~~~~~~~~~~~~~~~~~~~~
		473	+
		474	+As gamma-ray events are rare, the counts maps generated by ctbin will in
		475	+general be sparse, having many empty pixels, in particular at high energies.
		476	+An alternative analysis technique consists of working directly on the event
		477	+list without binning the events in a counts map. We will see the benefit of
		478	+such an analysis later once you re-run ctlike in unbinned mode.
		479	+
		480	+For unbinned analysis you first have to define the data space region over
		481	+which the analysis is done. This is similiar to the ctbin step in binned
		482	+analysis where you defined the size of the counts map, the energy range, and
		483	+the time interval. For unbinned analysis you have no such thing as a counts
		484	+map, but you have to define over which region of the data space the selected
		485	+events are spread (because the ctools have to integrate over this region to
		486	+compute the total number of predicted events in the data space that you
		487	+analyse). Furthermore, you have to define what energy range is covered, and
		488	+what time interval is spanned by the data. All this is done by the executable
		489	+ctselect, which replaces the ctbin step in an unbinned analysis.
		490	+
		491	+ctselect performs an event selection by choosing only events within a given
		492	+region-of-interest (ROI), within a given energy band, and within a given
		493	+time interval from the input event list. The ROI is a circular region on
		494	+the sky, for which you define the centre (in celestial coordinates) and the
		495	+radius. Such a circular ROI is sometimes also called an acceptance cone. The
		496	+following example shows how to run ctselect:
		497	+
		498	+
		499	+.. code-block:: bash
		500	+
		501	+ $ ctselect
		502	+ Input event list or observation definition file [test/data/crab_events.fits.gz] events.fits
		503	+ RA for ROI centre (degrees) (0-360) [83.63]
		504	+ Dec for ROI centre (degrees) (-90-90) [22.01]
		505	+ Radius of ROI (degrees) (0-180) [3.0]
		506	+ Start time (CTA MET in seconds) (0) [0.0]
		507	+ End time (CTA MET in seconds) (0) [0.0]
		508	+ Lower energy limit (TeV) (0) [0.1]
		509	+ Upper energy limit (TeV) (0) [100.0]
		510	+ Output event list or observation definition file [selected_events.fits]
		511	+
		512	+ctselect takes the input event list ``events.fits``, performs an event
		513	+selection, and writes the selected event into the file
		514	+``selected_events.fits``. The parameters it will query for are the centre
		515	+of the ROI, the radius of the ROI, the start and stop time (in seconds),
		516	+and the energy lower and upper limits (in TeV). The event selection information
		517	+is also written as a set of data selection keywords to the output events
		518	+file ``selected_events.fits``, by respecting the same syntax that has been
		519	+implemented for Fermi/LAT. The following image is a screen dump of the data
		520	+selection keywords that have been written to the ``EVENTS`` header in the
		521	+file ``selected_events.fits``:
		522	+
		523	+.. figure:: dskeys.jpg
		524	+ :width: 60%
		525	+
		526	+ Data selection keywords
		527	+
		528	+It is important for unbinned analysis that these data selection keywords
		529	+exist. If they don't exist, ctlike will not execute in unbinned mode.
		530	+
		531	+Below some lines of the ``ctselect.log`` file that show the data selection
		532	+part:
		533	+
		534	+.. code-block:: xml
		535	+
		536	+ 2014-01-09T21:36:20: +=================+
		537	+ 2014-01-09T21:36:20: \| Event selection \|
		538	+ 2014-01-09T21:36:20: +=================+
		539	+ 2014-01-09T21:36:20: === Observation ===
		540	+ 2014-01-09T21:36:21: Energy range ..............: 0.1 - 100 TeV
		541	+ 2014-01-09T21:36:21: Acceptance cone centre ....: RA=83.63, DEC=22.01 deg
		542	+ 2014-01-09T21:36:21: Acceptance cone radius ....: 3 deg
		543	+ 2014-01-09T21:36:21: cfitsio selection .........: ENERGY >= 0.10000000 && ENERGY <= 100.00000000 && ANGSEP(83.630000,22.010000,RA,DEC) <= 3.000000
		544	+ 2014-01-09T21:36:21: FITS filename .............: ctselectKKHEiF[EVENTS][ENERGY >= 0.10000000 && ENERGY <= 100.00000000 && ANGSEP(83.630000,22.010000,RA,DEC) <= 3.000000]
		545	+
		546	+Now that you have selected the events of interest, you can run ctlike in
		547	+unbinned mode. To do this you have to specify the selected event list
		548	+instead of the counts map:
		549	+
		550	+.. code-block:: bash
		551	+
		552	+ $ ctlike
		553	+ Event list, counts map or observation definition file [cntmap.fits] selected_events.fits
		554	+ Calibration database [$GAMMALIB/share/caldb/cta]
		555	+ Instrument response function [cta_dummy_irf]
		556	+ Source model [$GAMMALIB/share/models/crab.xml]
		557	+ Source model output file [crab_results.xml]
		558	+
		559	+You will recognise that ctlike runs much faster in unbinned mode compared
		560	+to binned mode. This is understandable as the selected event list contains
		561	+only 6127 events, while the binned counts map we used before had
		562	+200 x 200 x 20 = 800000 pixels. As unbinned maximum likelihood fitting loops
		563	+over the events (while binned maximum likelihood loops over the pixels),
		564	+there are much less operations to perform in unbinned than in binned mode
		565	+(there is some additional overhead in unbinned mode that comes from
		566	+integrating the models over the region of interest, yet this is negligible
		567	+compared to the operations needed when looping over the pixels). So as long
		568	+as you work with short event lists, unbinned mode is faster. Unbinned ctlike
		569	+should also be more precise as no binning is performed, hence there is no
		570	+loss of information due to histogramming.
		571	+
		572	+Below you see the corresponding output from the ctlike.log file. The fitted
		573	+parameters are essentially identical to the ones found in binned mode.
		574	+The slight difference with respect to the binned analysis may be explained
		575	+by the different event sample that has been used for the analysis: while
		576	+binned likelihood works on rectangular counts maps, unbinned likelihood works
		577	+on circular event selection regions. It is thus not possible to select exactly
		578	+the same events for both analyses.
		579	+
		580	+.. code-block:: xml
		581	+
		582	+ 2014-01-09T21:47:54: +=================================+
		583	+ 2014-01-09T21:47:54: \| Maximum likelihood optimisation \|
		584	+ 2014-01-09T21:47:54: +=================================+
		585	+ 2014-01-09T21:47:54: Initial iteration: func=44393.9, Lambda=0.001
		586	+ 2014-01-09T21:47:54: Iteration 1: func=44392.2, Lambda=0.001, delta=1.7034, max(grad)=9.60142 [6]
		587	+ 2014-01-09T21:47:54: Iteration 2: func=44392.1, Lambda=0.0001, delta=0.099482, max(grad)=0 [-1]
		588	+ 2014-01-09T21:47:54: Iteration 3: func=44392.1, Lambda=1e-05, delta=6.63429e-05, max(grad)=0.00559314 [6]
		589	+ 2014-01-09T21:47:54: Iteration 4: func=44392.1, Lambda=1e-06, delta=4.7803e-08, max(grad)=0 [-1]
		590	+ 2014-01-09T21:47:54:
		591	+ 2014-01-09T21:47:54: +=========================================+
		592	+ 2014-01-09T21:47:54: \| Maximum likelihood optimization results \|
		593	+ 2014-01-09T21:47:54: +=========================================+
		594	+ 2014-01-09T21:47:54: === GOptimizerLM ===
		595	+ 2014-01-09T21:47:54: Optimized function value ..: 44392.1
		596	+ 2014-01-09T21:47:54: Absolute precision ........: 1e-06
		597	+ 2014-01-09T21:47:54: Optimization status .......: converged
		598	+ 2014-01-09T21:47:54: Number of parameters ......: 9
		599	+ 2014-01-09T21:47:54: Number of free parameters .: 4
		600	+ 2014-01-09T21:47:54: Number of iterations ......: 4
		601	+ 2014-01-09T21:47:54: Lambda ....................: 1e-07
		602	+ 2014-01-09T21:47:54: Maximum log likelihood ....: -44392.1
		603	+ 2014-01-09T21:47:54: Observed events (Nobs) ...: 6127
		604	+ 2014-01-09T21:47:54: Predicted events (Npred) ..: 6127 (Nobs - Npred = 1.01613e-07)
		605	+ 2014-01-09T21:47:54: === GModels ===
		606	+ 2014-01-09T21:47:54: Number of models ..........: 2
		607	+ 2014-01-09T21:47:54: Number of parameters ......: 9
		608	+ 2014-01-09T21:47:54: === GModelSky ===
		609	+ 2014-01-09T21:47:54: Name ......................: Crab
		610	+ 2014-01-09T21:47:54: Instruments ...............: all
		611	+ 2014-01-09T21:47:54: Instrument scale factors ..: unity
		612	+ 2014-01-09T21:47:54: Observation identifiers ...: all
		613	+ 2014-01-09T21:47:54: Model type ................: PointSource
		614	+ 2014-01-09T21:47:54: Model components ..........: "SkyDirFunction" * "PowerLaw" * "Constant"
		615	+ 2014-01-09T21:47:54: Number of parameters ......: 6
		616	+ 2014-01-09T21:47:54: Number of spatial par's ...: 2
		617	+ 2014-01-09T21:47:54: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1)
		618	+ 2014-01-09T21:47:54: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1)
		619	+ 2014-01-09T21:47:54: Number of spectral par's ..: 3
		620	+ 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)
		621	+ 2014-01-09T21:47:54: Index ....................: -2.50585 +/- 0.0250821 [-0,-5] (free,scale=-1,gradient)
		622	+ 2014-01-09T21:47:54: PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient)
		623	+ 2014-01-09T21:47:54: Number of temporal par's ..: 1
		624	+ 2014-01-09T21:47:54: Constant .................: 1 (relative value) (fixed,scale=1,gradient)
		625	+ 2014-01-09T21:47:54: === GCTAModelRadialAcceptance ===
		626	+ 2014-01-09T21:47:54: Name ......................: Background
		627	+ 2014-01-09T21:47:54: Instruments ...............: CTA
		628	+ 2014-01-09T21:47:54: Instrument scale factors ..: unity
		629	+ 2014-01-09T21:47:54: Observation identifiers ...: all
		630	+ 2014-01-09T21:47:54: Model type ................: "Gaussian" * "FileFunction" * "Constant"
		631	+ 2014-01-09T21:47:54: Number of parameters ......: 3
		632	+ 2014-01-09T21:47:54: Number of radial par's ....: 1
		633	+ 2014-01-09T21:47:54: Sigma ....................: 3.04429 +/- 0.0330229 [0.01,10] deg2 (free,scale=1,gradient)
		634	+ 2014-01-09T21:47:54: Number of spectral par's ..: 1
		635	+ 2014-01-09T21:47:54: Normalization ............: 0.997516 +/- 0.0175361 [0,1000] (free,scale=1,gradient)
		636	+ 2014-01-09T21:47:54: Number of temporal par's ..: 1
		637	+ 2014-01-09T21:47:54: Constant .................: 1 (relative value) (fixed,scale=1,gradient)
		638	+ 2014-01-09T21:47:54:
		639	+ 2014-01-09T21:47:54: +==============+
		640	+ 2014-01-09T21:47:54: \| Save results \|
		641	+ 2014-01-09T21:47:54: +==============+
		642	+ 2014-01-09T21:47:54:
		643	+ 2014-01-09T21:47:54: Application "ctlike" terminated after 6 wall clock seconds, consuming 0.163384 seconds of CPU time.
		644	+
		645	+
		646	+Models available in ctools
		647	+--------------------------
		648	+
		649	+Using ctools from Python
		650	+------------------------
		651	+
		652	+Frequently asked questions
		653	+--------------------------
		654	+