How to connect observations to specific models?

In the previous examples there was always a single background model component to describe the residual particle background in the various dataset. This implies that the spatial and spectral shape of the background distribution is assumed to be identical for all observations. This is fine in a simulation, but for a real life situation this assumption will probably not hold.

Start by importing the gammalib, ctools, and cscripts Python modules.

import gammalib
import ctools
import cscripts

Simulating the dataset

Let’s start with creating a pointing definition ASCII file for two 30 min pointings near the Crab.

f = open('pnt.def', 'wb')
f.write('name,ra,dec,duration\n')
f.write('Crab,83.63,21.51,1800.0\n')
f.write('Crab,83.63,22.51,1800.0\n')
f.close()

Inspect the file that you just created. For that purpose let’s create a peek() function that will also be used later to display XML files.

def peek(filename):
    f = open(gammalib.expand_env(filename), 'r')
    for line in f:
        print(line.rstrip())
    f.close()
peek('pnt.def')

name,ra,dec,duration
Crab,83.63,21.51,1800.0
Crab,83.63,22.51,1800.0

A pointing definition file is an ASCII file in Comma Separated Values (CSV) format that specifies one pointing per row. The file provides the name of the observation, the Right Ascension ra and Declination dec of the pointing, and its duration. Additional optional columns are possible (defining for example the energy range or the Instrument Response Function), but for this simulation the provided information is sufficient.

Now transform the pointing definition file into an observation definition XML file using the csobsdef script.

obsdef = cscripts.csobsdef()
obsdef['inpnt']  = 'pnt.def'
obsdef['outobs'] = 'obs.xml'
obsdef['caldb']  = 'prod2'
obsdef['irf']    = 'South_0.5h'
obsdef.execute()

Let’s peek the resulting observation definition XML file

peek('obs.xml')

<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<observation_list title="observation list">
  <observation name="Crab" id="000001" instrument="CTA" statistic="cstat">
    <parameter name="Pointing" ra="83.63" dec="21.51" />
    <parameter name="GoodTimeIntervals" tmin="0" tmax="1800" />
    <parameter name="TimeReference" mjdrefi="51544" mjdreff="0.5" timeunit="s" timesys="TT" timeref="LOCAL" />
    <parameter name="Deadtime" deadc="0.98" />
    <parameter name="Calibration" database="prod2" response="South_0.5h" />
  </observation>
  <observation name="Crab" id="000002" instrument="CTA" statistic="cstat">
    <parameter name="Pointing" ra="83.63" dec="22.51" />
    <parameter name="GoodTimeIntervals" tmin="1800" tmax="3600" />
    <parameter name="TimeReference" mjdrefi="51544" mjdreff="0.5" timeunit="s" timesys="TT" timeref="LOCAL" />
    <parameter name="Deadtime" deadc="0.98" />
    <parameter name="Calibration" database="prod2" response="South_0.5h" />
  </observation>
</observation_list>

The file contains two observations which distinguish by the id attributes 000001 and 000002, which are unique to each observation for a given instrument. The id attribute can therefore be used to uniquely identify a CTA observation.

The value of the id attributes can be controlled by adding specific values to the pointing definition file, but if the values are missing - which is the case in the example - they simply count from 000001 upwards.

Feed now the observation definition XML file into the ctobssim tool to simulate the event data.

obssim = ctools.ctobssim()
obssim['inobs']     = 'obs.xml'
obssim['rad']       = 5.0
obssim['emin']      = 0.1
obssim['emax']      = 100.0
obssim['inmodel']   = '$CTOOLS/share/models/crab_2bkg.xml'
obssim['outevents'] = 'obs_2bkg.xml'
obssim.execute()

This will produce the two event files sim_events_000001.fits and sim_events_000002.fits on disk. Note that a specific model was used to simulate the data and peeking that model definiton file shows that it contains two different background components with a different power law Prefactor and Index. Both background components also have an id attribute which is used to tie them to the two observations. In other words, Background_000001 will be used for the observation with identifier 000001 and Background_000002 will be used for the observation with identifier 000002.

peek('$CTOOLS/share/models/crab_2bkg.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="Background_000001" type="CTAIrfBackground" instrument="CTA" id="000001">
    <spectrum type="PowerLaw">
      <parameter name="Prefactor"   scale="1.0"  value="0.5"  min="1e-3" max="1e+3"   free="1"/>
      <parameter name="Index"       scale="1.0"  value="0.2"  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 name="Background_000002" type="CTAIrfBackground" instrument="CTA" id="000002">
    <spectrum type="PowerLaw">
      <parameter name="Prefactor"   scale="1.0"  value="2.0"  min="1e-3" max="1e+3"   free="1"/>
      <parameter name="Index"       scale="1.0"  value="-0.2" 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>

Analysing the data

Now run a maximum likelihood fit of the model to the simulated data

like = ctools.ctlike()
like['inobs']    = 'obs_2bkg.xml'
like['inmodel']  = '$CTOOLS/share/models/crab_2bkg.xml'
like['outmodel'] = 'crab_results.xml'
like.run()

and inspect the model fitting results

print(like.obs().models())

=== GModels ===
 Number of models ..........: 3
 Number of parameters ......: 14
=== GModelSky ===
 Name ......................: Crab
 Instruments ...............: all
 Observation identifiers ...: all
 Model type ................: PointSource
 Model components ..........: "PointSource" * "PowerLaw" * "Constant"
 Number of parameters ......: 6
 Number of spatial par's ...: 2
  RA .......................: 83.6331 [-360,360] deg (fixed,scale=1)
  DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1)
 Number of spectral par's ..: 3
  Prefactor ................: 5.76556140039476e-16 +/- 7.29118825882603e-18 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
  Index ....................: -2.46505796935392 +/- 0.0107196310388945 [-0,-5]  (free,scale=-1,gradient)
  PivotEnergy ..............: 300000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)
 Number of temporal par's ..: 1
  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
 Number of scale par's .....: 0
=== GCTAModelIrfBackground ===
 Name ......................: Background_000001
 Instruments ...............: CTA
 Observation identifiers ...: 000001
 Model type ................: "PowerLaw" * "Constant"
 Number of parameters ......: 4
 Number of spectral par's ..: 3
  Prefactor ................: 0.504862658259721 +/- 0.0077214563116661 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient)
  Index ....................: 0.197685073818273 +/- 0.00982058375963094 [-5,5]  (free,scale=1,gradient)
  PivotEnergy ..............: 1000000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)
 Number of temporal par's ..: 1
  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
=== GCTAModelIrfBackground ===
 Name ......................: Background_000002
 Instruments ...............: CTA
 Observation identifiers ...: 000002
 Model type ................: "PowerLaw" * "Constant"
 Number of parameters ......: 4
 Number of spectral par's ..: 3
  Prefactor ................: 2.01944300996949 +/- 0.0173357023683058 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient)
  Index ....................: -0.19787756674158 +/- 0.00509621480931608 [-5,5]  (free,scale=1,gradient)
  PivotEnergy ..............: 1000000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)
 Number of temporal par's ..: 1
  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)

Obviously the two background models have different best fitting spectral parameters which correspond to the simulated values (see above).