Simulating event data¶

What you will learn

You will learn how to use the ctobssim tool from the command line to simulate event data.

ctobssim draws random events using a model of the celestial gamma-ray intensity distribution that was convolved with the instrument response function. The tool adds random events from the expected distribution of the residual background to the data.

Let’s assume that you want to generate a simulated CTA event list. You do this by using the ctobssim tool from the command line.

To invoke the tool, type ctobssim at the console prompt (which is denoted by $). ctobssim will then query for several parameters that define the characteristics of the simulation:

$ ctobssim
RA of pointing (degrees) (0-360) [83.63]
Dec of pointing (degrees) (-90-90) [22.51]
Radius of FOV (degrees) (0-180) [5.0]
Start time (UTC string, JD, MJD or MET in seconds) [2020-01-01T00:00:00]
Stop time (UTC string, JD, MJD or MET in seconds) [2020-01-01T00:30:00] 2020-01-01T01:00:00
Lower energy limit (TeV) [0.1] 0.03
Upper energy limit (TeV) [100.0] 150.0
Calibration database [prod2]
Instrument response function [South_0.5h]
Input model definition XML file [$CTOOLS/share/models/crab.xml]
Output event data file or observation definition XML file [events.fits]

Each line represents a query for one parameter value. The line starts with a short description of the parameter, followed by the default value proposed by ctobssim in squared brackets [ ].

If no parameter value is entered the default value will be used. Otherwise, the specified value will overwrite the default value. The round brackets ( ) indicate the range of possible parameter values (if applicable).

In the example above the stop time and the energy range were specified, for the remaining parameters the default value was used.

Note

Times can be entered in various formats in ctools. Times can be provided as

UTC date strings (e.g. "2020-01-01T00:00:00")
Modified Julian days (e.g. "MJD 58849.0")
Julian days (e.g. "JD 2458849.5")
Mission elapsed time in seconds, counted from 2000-01-01T12:00:00 (e.g. "631108869.18")

Usually times are given in Terrestial Time (TT) but may also be specified as TAI or UTC (e.g. "MJD 58849.0 (TAI)" or "MJD 58849.0 (UTC)").

You may have recognised that the environment variable $CTOOLS has been used in the path name of the model. ctools will automatically expand the environment variables in parameter values.

The CTA instrument response function (IRF) is taken from the prod2 database. The response for the southern array using the cuts optimised for 0.5 hours of observing time is used.

Note

ctools comes bundled with CTA instrument response functions for the northern and the southern array. The IRFs are based on a prod2 analysis with cuts optimised for 0.5 hours, 5 hours and 50 hours of observing time. Be aware that these times do not need to correspond to the actual observing time that is simulated. IRFs optimised for short observation times correspond to loose cuts that keep more photons but also more background events, while IRFs optimised for long observation times correspond to hard cuts that limit the number of background events at the expense of loosing some photons. The following IRFs are available: North_0.5h, North_5h, North_50h, South_0.5h, South_5h, and South_50h.

Events are simulated based on the instrument response function 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 one hour. Events are simulated for energies between 30 GeV and 200 TeV.

The source and background model is defined by the model definition XML file $CTOOLS/share/models/crab.xml:

<?xml version="1.0" standalone="no"?>
<source_library title="source library">
  <source name="Crab" type="PointSource">
    <spectrum type="PowerLaw">
       <parameter name="Prefactor"   scale="1e-16" value="5.7"  min="1e-07" max="1000.0" free="1"/>
       <parameter name="Index"       scale="-1"    value="2.48" min="0.0"   max="+5.0"   free="1"/>
       <parameter name="PivotEnergy" scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="0"/>
    </spectrum>
    <spatialModel type="PointSource">
      <parameter name="RA"  scale="1.0" value="83.6331" min="-360" max="360" free="0"/>
      <parameter name="DEC" scale="1.0" value="22.0145" min="-90"  max="90"  free="0"/>
    </spatialModel>
  </source>
  <source name="CTABackgroundModel" type="CTAIrfBackground" instrument="CTA">
    <spectrum type="PowerLaw">
      <parameter name="Prefactor"   scale="1.0"  value="1.0"  min="1e-3" max="1e+3"   free="1"/>
      <parameter name="Index"       scale="1.0"  value="0.0"  min="-5.0" max="+5.0"   free="1"/>
      <parameter name="PivotEnergy" scale="1e6"  value="1.0"  min="0.01" max="1000.0" free="0"/>
    </spectrum>
  </source>
</source_library>

The model consists of a source library that contains two components: the Crab nebula and an instrumental background model.

The Crab nebula is modelled by a factorized sky model that has a spectral and a spatial component (tags <spectrum> and <spatialModel>, respectively). The spectrum is modelled by a power law, which is defined by three parameters: the Prefactor, the Index and the Scale. The spatial model has two parameters: Right Ascension in degrees (RA), and Declination in degrees (DEC). Each parameter has a value and a scale factor, the real value of the parameter being the product value * scale. Typically, scale is chosen so that value is of the order of 1 (this is relevant for model fitting). In addition, value is bound by a minimum (min) and maximum (max) value, and a parameter may be free (free="1") or fixed (free="0"). The min, max, and free attributes are not relevant here for the simulations, but they will be important for the model fitting later.

The spectral intensity I(E) (in units of ${\rm photons} \, {\rm cm}^{-2} \, {\rm s}^{-1} \, {\rm MeV}^{-1}$) of the power law is given by

\[\frac{dN}{dE} = N_0 \left( \frac{E}{E_0} \right)^{\gamma}\]

where the parameters in the XML definition have the following mappings:

$N_0$ = Prefactor
$\gamma$ = Index
$E_0$ = PivotEnergy

Warning

Energies are given in the XML file in MeV units. This is a GammaLib convention that can not be modified. So make sure you always use MeV as energy unit in an XML file.

Note

As customary for IACT observations, the pointing direction was slightly offset from the source of interest, i.e., the Crab. This makes it possible to better handle systematics due to the limited knowledge of the instrumental background.

The instrumental background of CTA is modelled using the background information provided in the instrument response function (CTAIrfBackground), where the energy dependence of the background model is multipled by a power law. As it is defined here, the power law represents a constant of 1, hence the background IRF will be used without any modification. The power law will become active when fitting the data later and allows a spectral adjustment of the background model that may account for uncertainties in the background information provided in the IRF.

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:

$ ctobssim seed=41

Note

Hidden parameters are parameters that are not queried by a tool since in general their values is not expected to change frequently. To change hidden parameters they have to be given as arguments on the command line. Multiple hidden parameters need to be separated by a white space.

ctobssim will write two files in the working directory: events.fits and ctobssim.log. The first file contains the simulated events in FITS format and can be inspected using fv or ds9. The FITS file will contain three extensions: an empty primary image, a binary table named EVENTS that holds the events (one row per event), and a binary table named GTI holding the Good Time Intervals (for the moment a single row with two 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:

2019-04-02T12:47:42: === CTA observation ===
2019-04-02T12:47:42:  Simulation cone ...........: RA=83.63 deg, Dec=22.51 deg, radius=5.5 deg
2019-04-02T12:47:42:  Time interval .............: 6.31109e+08 - 6.31112e+08 s
2019-04-02T12:47:42:  Photon energy range .......: 30 GeV - 70.310187347763 GeV
2019-04-02T12:47:42:  Event energy range ........: 30 GeV - 70.310187347763 GeV
2019-04-02T12:47:42:   Simulation area ..........: 1.90195e+09 cm2
2019-04-02T12:47:42:   Use model ................: Crab
2019-04-02T12:47:42:   Normalization ............: 1 [Crab]
2019-04-02T12:47:42:   Flux .....................: 2.50006e-09 [Crab] photons/cm2/s
2019-04-02T12:47:42:   Normalized flux ..........: 2.50006e-09 [Crab] photons/cm2/s
2019-04-02T12:47:42:   Photon rate ..............: 4.75499 photons/s [Crab]
2019-04-02T12:47:42:   MC source photons ........: 17169 [Crab]
2019-04-02T12:47:42:   MC source events .........: 3870 [Crab]
2019-04-02T12:47:42:   MC source events .........: 3870 (all source models)
2019-04-02T12:47:42:  Photon energy range .......: 70.310187347763 GeV - 164.784081495918 GeV
...
2019-04-02T12:47:42:  MC source photons .........: 46368 [Crab]
2019-04-02T12:47:42:  MC source events ..........: 12749 [Crab]
2019-04-02T12:47:50:  MC events outside ROI .....: 0
2019-04-02T12:47:50:  MC background events ......: 193299
2019-04-02T12:47:50:  MC identifier 1 ...........: Crab
2019-04-02T12:47:50:  MC identifier 2 ...........: CTABackgroundModel
2019-04-02T12:47:50:  MC events .................: 206048 (all models)

Each line starts with the UTC time at which the line has been written. In this run, 46368 Crab photons have been thrown. 12749 of these photons have been registered by CTA as events. In the same time interval, 193299 background events have been registred by CTA.

Note

ctobssim will split the simulated energy range into a number of slices, controlled via the hidden eslices parameter (ten energy slices are used by default). For each energy slice, the simulation area will be adapted to the effective area of the array in that energy slice, which helps to keep the computing time low. The log file will provide information about the simulation in each slice. In the example above, the simulation results for the first energy slice are shown, followed by a summary of the results for all slices.

You may change the name of the log file using the hidden parameter logfile:

$ 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:

$ ctobssim chatter=4

chatter can vary between 0 and 4, 0 providing no information while 4 provides the most detailed information.

By default, all ctools have a chatter level of 2.

You may also duplicate the log file information into the console by setting the hidden debug parameter to yes:

$ ctobssim debug=yes

Note

All tools have the hidden parameters logfile, chatter, and debug and you can use these parameters to control the log file output. In addition, all tools have the hidden parameter clobber that allows to overwrite existing files (set to yes by default) and mode that defines the mode of automatic parameters (set to ql for query and learn by default).

Previous topic

Next topic

Simulating event data¶