As gamma-ray events are rare, the counts cubes generated by ctbin will in general be sparsly populated, having many empty pixels, in particular at high energies. An alternative analysis technique consists of working directly on the event list without binning the events in a counts cube. We will see the benefit of such an analysis later once you re-run 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 for a binned analysis where you defined the size of the counts cube, the energy range, and the time interval. For unbinned analysis you have no such thing as a counts cube, but you have to define over which region of the data space the selected events are spread (because the ctools have to integrate over this region to compute the total number of predicted events in the data space that you analyse). Furthermore, you have to define what energy range is covered, and what time interval is spanned by the data. All this is done by the ctselect tool, 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:
$ ctselect
Input event list or observation definition XML file [events.fits]
RA for ROI centre (degrees) (0-360) [83.63]
Dec for ROI centre (degrees) (-90-90) [22.01]
Radius of ROI (degrees) (0-180) [3.0]
Start time (CTA MET in seconds) [0.0]
End time (CTA MET in seconds) [0.0]
Lower energy limit (TeV) [0.1]
Upper energy limit (TeV) [100.0]
Output event list or observation definition XML file [selected_events.fits]
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:
Data selection keywords
It is mandatory for an unbinned analysis that these data selection keywords exist in the FITS file. If they don’t exist, ctlike will not execute in unbinned mode.
Below some lines of the ctselect.log file that show the data selection part:
2015-12-06T23:44:52: +=================+
2015-12-06T23:44:52: | Event selection |
2015-12-06T23:44:52: +=================+
2015-12-06T23:44:52: === CTA observation ===
2015-12-06T23:44:52: Selected energy range .....: 0.1 - 100 TeV
2015-12-06T23:44:52: Requested ROI .............: Centre(RA,DEC)=(83.63, 22.01) deg, Radius=3 deg
2015-12-06T23:44:52: ROI of data ...............: Centre(RA,DEC)=(83.63, 22.01) deg, Radius=5 deg
2015-12-06T23:44:52: Selected ROI ..............: Centre(RA,DEC)=(83.63, 22.01) deg, Radius=3 deg
2015-12-06T23:44:52: cfitsio selection .........: ENERGY >= 0.10000000 && ENERGY <= 100.00000000 && ANGSEP(83.630000,22.010000,RA,DEC) <= 3.000000
2015-12-06T23:44:52: FITS filename .............: /var/tmp/tmp.0.jHg5hJ[EVENTS][ENERGY >= 0.10000000 && ENERGY <= 100.00000000 && ANGSEP(83.630000,22.010000,RA,DEC) <= 3.000000]
Note
ctobssim will also write data selection keywords in the event list FITS file, hence you can run ctlike directly on a FITS file produced by ctobssim. Any selection performed by ctselect needs to be fully enclosed within any previous selection, e.g. the ROI needs to be fully enclosed in the acceptance cone used for event simulation, the energy selection must be fully comprised in the range of simulated energies, the same applies for the temporal selection. ctselect will automatically adjust the selection parameters to guarantee full enclosure. To keep track of this adjustment, the ctselect log file quotes the requested selection, any existing selections, and the selection that was finally applied.
Warning
ctselect may of course also be used for event selection prior to binned analysis, for example to select events for a given period in time. If you use ctselect however to make a spatial or energy selection, make sure that the counts cube is fully enclosed in the selection intervals. Otherwise you will get empty zones in the counts cube of which the ctools are not aware of, and the subsequent analysis results will be wrong.
Now that you have selected the events of interest, you can run ctlike in unbinned mode. To do this you have to specify the selected event list instead of the counts cube:
$ ctlike
Input event list, counts cube or observation definition XML file [cntcube.fits] selected_events.fits
Calibration database [prod2]
Instrument response function [South_0.5h]
Input model XML file [$CTOOLS/share/models/crab.xml]
Output model XML file [crab_results.xml]
You will recognise that ctlike runs much faster in unbinned mode compared to binned mode. This is understandable as the selected event list contains only 21203 events, while the binned counts cube we used before had 200 x 200 x 20 = 800000 pixels. As unbinned maximum likelihood fitting loops over the events (while binned maximum likelihood loops over the pixels), there are much less operations to perform in unbinned than in binned mode (there is some additional overhead in unbinned mode that comes from integrating the models over the region of interest, yet this is negligible compared to the operations needed when looping over all pixels). So as long as you work with small event lists, unbinned mode is faster (this typically holds up to 100 hours of observing time). Unbinned ctlike should also be more precise as no binning is performed, hence there is no loss of information due to histogramming.
Below you see the corresponding output from the ctlike.log file. The fitted parameters are essentially identical to the ones found in binned mode. The slight difference with respect to the binned analysis may be explained by the different event sample that has been used for the analysis: while binned likelihood works on rectangular counts cubes, unbinned likelihood works on circular event selection regions. It is thus not possible to select exactly the same events for both analyses.
2015-12-07T21:02:17: +=================================+
2015-12-07T21:02:17: | Maximum likelihood optimisation |
2015-12-07T21:02:17: +=================================+
2015-12-07T21:02:17: >Iteration 0: -logL=135278.644, Lambda=1.0e-03
2015-12-07T21:02:17: >Iteration 1: -logL=135274.475, Lambda=1.0e-03, delta=4.168, max(|grad|)=5.806873 [Index:7]
2015-12-07T21:02:17: >Iteration 2: -logL=135274.474, Lambda=1.0e-04, delta=0.001, max(|grad|)=-0.009463 [Index:7]
2015-12-07T21:02:17:
...
2015-12-07T21:02:17: +=========================================+
2015-12-07T21:02:17: | Maximum likelihood optimisation results |
2015-12-07T21:02:17: +=========================================+
2015-12-07T21:02:17: === GOptimizerLM ===
2015-12-07T21:02:17: Optimized function value ..: 135274.474
2015-12-07T21:02:17: Absolute precision ........: 0.005
2015-12-07T21:02:17: Acceptable value decrease .: 2
2015-12-07T21:02:17: Optimization status .......: converged
2015-12-07T21:02:17: Number of parameters ......: 10
2015-12-07T21:02:17: Number of free parameters .: 4
2015-12-07T21:02:17: Number of iterations ......: 2
2015-12-07T21:02:17: Lambda ....................: 1e-05
2015-12-07T21:02:17: Maximum log likelihood ....: -135274.474
2015-12-07T21:02:17: Observed events (Nobs) ...: 21203.000
2015-12-07T21:02:17: Predicted events (Npred) ..: 21202.998 (Nobs - Npred = 0.00198937)
2015-12-07T21:02:17: === GModels ===
2015-12-07T21:02:17: Number of models ..........: 2
2015-12-07T21:02:17: Number of parameters ......: 10
2015-12-07T21:02:17: === GModelSky ===
2015-12-07T21:02:17: Name ......................: Crab
2015-12-07T21:02:17: Instruments ...............: all
2015-12-07T21:02:17: Instrument scale factors ..: unity
2015-12-07T21:02:17: Observation identifiers ...: all
2015-12-07T21:02:17: Model type ................: PointSource
2015-12-07T21:02:17: Model components ..........: "SkyDirFunction" * "PowerLaw" * "Constant"
2015-12-07T21:02:17: Number of parameters ......: 6
2015-12-07T21:02:17: Number of spatial par's ...: 2
2015-12-07T21:02:17: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1)
2015-12-07T21:02:17: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1)
2015-12-07T21:02:17: Number of spectral par's ..: 3
2015-12-07T21:02:17: Prefactor ................: 5.81583e-16 +/- 1.01727e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
2015-12-07T21:02:17: Index ....................: -2.50611 +/- 0.0153542 [-0,-5] (free,scale=-1,gradient)
2015-12-07T21:02:17: PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient)
2015-12-07T21:02:17: Number of temporal par's ..: 1
2015-12-07T21:02:17: Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
2015-12-07T21:02:17: === GCTAModelIrfBackground ===
2015-12-07T21:02:17: Name ......................: CTABackgroundModel
2015-12-07T21:02:17: Instruments ...............: CTA
2015-12-07T21:02:17: Instrument scale factors ..: unity
2015-12-07T21:02:17: Observation identifiers ...: all
2015-12-07T21:02:17: Model type ................: "PowerLaw" * "Constant"
2015-12-07T21:02:17: Number of parameters ......: 4
2015-12-07T21:02:17: Number of spectral par's ..: 3
2015-12-07T21:02:17: Prefactor ................: 0.990866 +/- 0.0136406 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient)
2015-12-07T21:02:17: Index ....................: 0.00536014 +/- 0.00839916 [-5,5] (free,scale=1,gradient)
2015-12-07T21:02:17: PivotEnergy ..............: 1e+06 [10000,1e+09] MeV (fixed,scale=1e+06,gradient)
2015-12-07T21:02:17: Number of temporal par's ..: 1
2015-12-07T21:02:17: Normalization ............: 1 (relative value) (fixed,scale=1,gradient)