Performing a classical On/Off analysis

In this tutorial you will learn to perform a classical On/Off analysis of the data.

We start by importing the gammalib, ctools, and cscripts Python modules.

import gammalib
import ctools
import cscripts

We will also use matplotlib to display the results.

%matplotlib inline
import matplotlib.pyplot as plt

And finally we add to our path the directory containing the example plotting scripts provided with the ctools installation.

import sys
import os
sys.path.append(os.environ['CTOOLS']+'/share/examples/python/')

Classical analysis without background model

Preparation of the On/Off binned data

The first step for a classical analysis is to bin the events in bins of energy for an On (signal) region, and one or several Off (background regions).

We will use the event selection performed in the previous tutorial.

obsfile = 'obs_crab_selected.xml'
phagen  = cscripts.csphagen()
phagen['inobs']         = obsfile
phagen['inmodel']       = 'NONE' # assume that the source is pointlike
phagen['ebinalg']       = 'LOG'
phagen['emin']          = 0.66
phagen['emax']          = 30.0
phagen['enumbins']      = 20
phagen['coordsys']      = 'CEL'
phagen['ra']            = 83.63
phagen['dec']           = 22.01
phagen['rad']           = 0.2
phagen['bkgmethod']     = 'REFLECTED'
phagen['use_model_bkg'] = False
phagen['maxoffset']     = 2.0
phagen['stack']         = True
phagen['outobs']        = 'obs_crab_onoff.xml'
phagen['outmodel']      = 'onoff_model.xml'
phagen.execute()

Note that we have used the hidden parameter use_model_bkg to specify that no background model should be used in the computation of the background normalisation factor, which is simply assumed to be the ratio of the solid angles of On/Off regions.

Let us peek at the resulting set of On/Off observations.

print(phagen.obs()[0])

=== GCTAOnOffObservation ===
 Name ......................:
 Identifier ................:
 Instrument ................: HESSOnOff
 Statistic .................: wstat
 Ontime ....................: 6313.8117676 s
 Livetime ..................: 6313.8117676 s
 Deadtime correction .......: 1
=== GPha ===
 Exposure ..................: 6313.8117676 s
 Number of bins ............: 20
 Energy range ..............: 660 GeV - 30 TeV
 Observation energy range ..: 660.693448007596 GeV - 100 TeV
 Total number of counts ....: 576
 Underflow counts ..........: 0
 Overflow counts ...........: 1
 Outflow counts ............: 0
=== GPha ===
 Exposure ..................: 6313.8117676 s
 Number of bins ............: 20
 Energy range ..............: 660 GeV - 30 TeV
 Observation energy range ..: 660.693448007596 GeV - 100 TeV
 Total number of counts ....: 772
 Underflow counts ..........: 0
 Overflow counts ...........: 4
 Outflow counts ............: 0
=== GArf ===
 Number of bins ............: 61
 Energy range ..............: 330 GeV - 36.000002048 TeV
=== GRmf ===
 Number of true energy bins : 61
 Number of measured bins ...: 20
 True energy range .........: 330 GeV - 36.000002048 TeV
 Measured energy range .....: 660 GeV - 30 TeV

Note that the instrument is now HESSOnOff, to signify that we are dealing with On/Off observations. Since we have no background model the statistic has been set to WSTAT. The script has also created a model for the On region.

print(phagen.obs().models())

=== GModels ===
 Number of models ..........: 1
 Number of parameters ......: 6
=== GModelSky ===
 Name ......................: Dummy
 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.63 deg (fixed,scale=1)
  DEC ......................: 22.01 deg (fixed,scale=1)
 Number of spectral par's ..: 3
  Prefactor ................: 1e-18 +/- 0 [0,infty[ ph/cm2/s/MeV (free,scale=1e-18,gradient)
  Index ....................: -2 +/- 0 [10,-10]  (free,scale=-2,gradient)
  PivotEnergy ..............: 1000000 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

Since we do not have specified any input model the script has placed a point source named Dummy, at the center of the On region with a power law spectrum. We will rename this source as Crab. At this point you can also modify the spectral model, e.g., to use more appropriate parameter values or a different spectral shape.

phagen.obs().models()['Dummy'].name('Crab')

Likelihood fit and residuals

We can now fit this model to the data.

like = ctools.ctlike(phagen.obs())
like.run()

Let’s look at the results from the optimisation and the fitted model.

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

=== GOptimizerLM ===
 Optimized function value ..: 15.746
 Absolute precision ........: 0.005
 Acceptable value decrease .: 2
 Optimization status .......: converged
 Number of parameters ......: 6
 Number of free parameters .: 2
 Number of iterations ......: 11
 Lambda ....................: 0.0001
=== GModels ===
 Number of models ..........: 1
 Number of parameters ......: 6
=== 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.63 deg (fixed,scale=1)
  DEC ......................: 22.01 deg (fixed,scale=1)
 Number of spectral par's ..: 3
  Prefactor ................: 4.61852950191309e-17 +/- 2.97687577951453e-18 [0,infty[ ph/cm2/s/MeV (free,scale=1e-18,gradient)
  Index ....................: -2.63298373450669 +/- 0.0795120102122515 [10,-10]  (free,scale=-2,gradient)
  PivotEnergy ..............: 1000000 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

The fit has properly converged and the results are consistent with those obtained with the unbinned analysis.

We will also check the spectral residuals.

residuals = 'residuals_classical_nobkg.fits'
res = cscripts.csresspec(like.obs())
res['components'] = True
res['algorithm']  = 'SIGNIFICANCE'
res['outfile']    = residuals
res.execute()
from show_residuals import plot_residuals
plot_residuals(residuals,'',0)

Classical analysis with a background model

Preparation of the On/Off binned data

If a background model is available this can be used to compute the expected background rates for different observations as a function of energy and normalise the number of background counts between On and Off regions. Here we will use the background model derived using csbkgmodel.

srcmodel   = os.environ['CTOOLS']+'/share/models/crab_nobkg.xml'
bkgmodel   = 'bkgmodel.xml'
startmodel = 'model_classical_bkg.xml'
merge = cscripts.csmodelmerge()
merge['inmodels'] = srcmodel + ' ' + bkgmodel
merge['outmodel'] = startmodel
merge.execute()

Now we can run csphagen using this model including the background in input. Note that since the background model is defined per observation we will not stack the observations together.

phagen = cscripts.csphagen()
phagen['inobs']         = obsfile
phagen['inmodel']       = startmodel
phagen['srcname']       = 'Crab'
phagen['ebinalg']       = 'LOG'
phagen['emin']          = 0.66
phagen['emax']          = 30.0
phagen['enumbins']      = 20
phagen['coordsys']      = 'CEL'
phagen['ra']            = 83.63
phagen['dec']           = 22.01
phagen['rad']           = 0.2
phagen['bkgmethod']     = 'REFLECTED'
phagen['use_model_bkg'] = True
phagen['maxoffset']     = 2.0
phagen['stack']         = False # treat observations separately in joint likelihood fit
phagen.run()

Likelihood fit and residuals

We will now fit the model to the data and look at the results.

like = ctools.ctlike(phagen.obs())
like.run()
print(like.opt())
print(like.obs().models())

=== GOptimizerLM ===
 Optimized function value ..: -2417.944
 Absolute precision ........: 0.005
 Acceptable value decrease .: 2
 Optimization status .......: converged
 Number of parameters ......: 86
 Number of free parameters .: 42
 Number of iterations ......: 29
 Lambda ....................: 0.0001
=== GModels ===
 Number of models ..........: 5
 Number of parameters ......: 86
=== 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 ................: 4.44907776007907e-17 +/- 2.98566064292929e-18 [1e-24,1e-14] ph/cm2/s/MeV (free,scale=1e-17,gradient)
  Index ....................: -2.64019099792601 +/- 0.0844654299892131 [-0,-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)
 Number of scale par's .....: 0
=== GCTAModelBackground ===
 Name ......................: Background_023523
 Instruments ...............: HESSOnOff
 Observation identifiers ...: 023523
 Model type ................: CTABackground
 Model components ..........: "Multiplicative" * "NodeFunction" * "Constant"
 Number of parameters ......: 20
 Number of spatial par's ...: 3
  1:Normalization ..........: 1 [0.001,1000]  (fixed,scale=1,gradient)
  2:Grad_DETX ..............: 0.0881748792818854 +/- 0 deg^-1 (free,scale=1,gradient)
  2:Grad_DETY ..............: -0.000701218995455618 +/- 0 deg^-1 (free,scale=1,gradient)
 Number of spectral par's ..: 16
  Energy0 ..................: 700000 [7e-05,infty[ MeV (fixed,scale=700000)
  Intensity0 ...............: 1.89896487985208 +/- 3.09950981249023 [1.00182541969977e-13,infty[ ph/cm2/s/MeV (free,scale=0.00100182541969977,gradient)
  Energy1 ..................: 1197413.67463685 [0.000119741367463685,infty[ MeV (fixed,scale=1197413.67463685)
  Intensity1 ...............: 0.754995013333362 +/- 0.267674756209559 [3.21240099713853e-14,infty[ ph/cm2/s/MeV (free,scale=0.000321240099713853,gradient)
  Energy2 ..................: 2048285.01172474 [0.000204828501172474,infty[ MeV (fixed,scale=2048285.01172474)
  Intensity2 ...............: 1.17347310940125 +/- 0.327498259113055 [1.03007170346199e-14,infty[ ph/cm2/s/MeV (free,scale=0.000103007170346199,gradient)
  Energy3 ..................: 3503777.83227556 [0.000350377783227556,infty[ MeV (fixed,scale=3503777.83227556)
  Intensity3 ...............: 1.39420069744212 +/- 0.404449186550592 [3.30297405342054e-15,infty[ ph/cm2/s/MeV (free,scale=3.30297405342054e-05,gradient)
  Energy4 ..................: 5993530.69893743 [0.000599353069893743,infty[ MeV (fixed,scale=5993530.69893743)
  Intensity4 ...............: 1.32748213047069 +/- 0.551998524503437 [1.05911438600856e-15,infty[ ph/cm2/s/MeV (free,scale=1.05911438600856e-05,gradient)
  Energy5 ..................: 10252479.454662 [0.0010252479454662,infty[ MeV (fixed,scale=10252479.454662)
  Intensity5 ...............: 1.53590487400556 +/- 0.862929853325373 [3.39610080039421e-16,infty[ ph/cm2/s/MeV (free,scale=3.39610080039421e-06,gradient)
  Energy6 ..................: 17537798.7113509 [0.00175377987113509,infty[ MeV (fixed,scale=17537798.7113509)
  Intensity6 ...............: 1.26422194838399 +/- 1.34653218966622 [1.08897592165697e-16,infty[ ph/cm2/s/MeV (free,scale=1.08897592165697e-06,gradient)
  Energy7 ..................: 30000000 [0.003,infty[ MeV (fixed,scale=30000000)
  Intensity7 ...............: 6.77476052324258 +/- 8.56245789647827 [3.49185323889972e-17,infty[ ph/cm2/s/MeV (free,scale=3.49185323889972e-07,gradient)
 Number of temporal par's ..: 1
  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
=== GCTAModelBackground ===
 Name ......................: Background_023526
 Instruments ...............: HESSOnOff
 Observation identifiers ...: 023526
 Model type ................: CTABackground
 Model components ..........: "Multiplicative" * "NodeFunction" * "Constant"
 Number of parameters ......: 20
 Number of spatial par's ...: 3
  1:Normalization ..........: 1 [0.001,1000]  (fixed,scale=1,gradient)
  2:Grad_DETX ..............: 0.0153166143265463 +/- 0 deg^-1 (free,scale=1,gradient)
  2:Grad_DETY ..............: -0.00908019181469828 +/- 0 deg^-1 (free,scale=1,gradient)
 Number of spectral par's ..: 16
  Energy0 ..................: 700000 [7e-05,infty[ MeV (fixed,scale=700000)
  Intensity0 ...............: 1.30442623250239 +/- 0.615959779856916 [9.62519069451052e-14,infty[ ph/cm2/s/MeV (free,scale=0.000962519069451052,gradient)
  Energy1 ..................: 1197413.67463685 [0.000119741367463685,infty[ MeV (fixed,scale=1197413.67463685)
  Intensity1 ...............: 0.566084707224135 +/- 0.189571117408801 [2.95441986174209e-14,infty[ ph/cm2/s/MeV (free,scale=0.000295441986174209,gradient)
  Energy2 ..................: 2048285.01172474 [0.000204828501172474,infty[ MeV (fixed,scale=2048285.01172474)
  Intensity2 ...............: 0.955930666245828 +/- 0.272410595459555 [9.068492247571e-15,infty[ ph/cm2/s/MeV (free,scale=9.068492247571e-05,gradient)
  Energy3 ..................: 3503777.83227556 [0.000350377783227556,infty[ MeV (fixed,scale=3503777.83227556)
  Intensity3 ...............: 1.91071527979174 +/- 0.589677088986724 [2.78354314866283e-15,infty[ ph/cm2/s/MeV (free,scale=2.78354314866283e-05,gradient)
  Energy4 ..................: 5993530.69893743 [0.000599353069893743,infty[ MeV (fixed,scale=5993530.69893743)
  Intensity4 ...............: 0.796340097360303 +/- 0.511186601597472 [8.54399193266454e-16,infty[ ph/cm2/s/MeV (free,scale=8.54399193266454e-06,gradient)
  Energy5 ..................: 10252479.454662 [0.0010252479454662,infty[ MeV (fixed,scale=10252479.454662)
  Intensity5 ...............: 1.31629202714156 +/- 0.980251656539197 [2.62254954375342e-16,infty[ ph/cm2/s/MeV (free,scale=2.62254954375342e-06,gradient)
  Energy6 ..................: 17537798.7113509 [0.00175377987113509,infty[ MeV (fixed,scale=17537798.7113509)
  Intensity6 ...............: 2.27556128029012 +/- 3.06865755885414 [8.04982748537824e-17,infty[ ph/cm2/s/MeV (free,scale=8.04982748537824e-07,gradient)
  Energy7 ..................: 30000000 [0.003,infty[ MeV (fixed,scale=30000000)
  Intensity7 ...............: 2.4708674312253e-17 +/- 1.05696632580831e-13 [2.4708674312253e-17,infty[ ph/cm2/s/MeV (free,scale=2.4708674312253e-07,gradient)
 Number of temporal par's ..: 1
  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
=== GCTAModelBackground ===
 Name ......................: Background_023559
 Instruments ...............: HESSOnOff
 Observation identifiers ...: 023559
 Model type ................: CTABackground
 Model components ..........: "Multiplicative" * "NodeFunction" * "Constant"
 Number of parameters ......: 20
 Number of spatial par's ...: 3
  1:Normalization ..........: 1 [0.001,1000]  (fixed,scale=1,gradient)
  2:Grad_DETX ..............: 0.018749400342 +/- 0 deg^-1 (free,scale=1,gradient)
  2:Grad_DETY ..............: -0.126640342510319 +/- 0 deg^-1 (free,scale=1,gradient)
 Number of spectral par's ..: 16
  Energy0 ..................: 700000 [7e-05,infty[ MeV (fixed,scale=700000)
  Intensity0 ...............: 0.885669706791242 +/- 0.219147786869672 [8.7418891142236e-14,infty[ ph/cm2/s/MeV (free,scale=0.00087418891142236,gradient)
  Energy1 ..................: 1197413.67463685 [0.000119741367463685,infty[ MeV (fixed,scale=1197413.67463685)
  Intensity1 ...............: 1.12011463147092 +/- 0.162106623373511 [2.69181956755491e-14,infty[ ph/cm2/s/MeV (free,scale=0.000269181956755491,gradient)
  Energy2 ..................: 2048285.01172474 [0.000204828501172474,infty[ MeV (fixed,scale=2048285.01172474)
  Intensity2 ...............: 0.747992399165037 +/- 0.145060618647741 [8.28870337932099e-15,infty[ ph/cm2/s/MeV (free,scale=8.28870337932099e-05,gradient)
  Energy3 ..................: 3503777.83227556 [0.000350377783227556,infty[ MeV (fixed,scale=3503777.83227556)
  Intensity3 ...............: 0.750550500933036 +/- 0.184075541493125 [2.55227373106486e-15,infty[ ph/cm2/s/MeV (free,scale=2.55227373106486e-05,gradient)
  Energy4 ..................: 5993530.69893743 [0.000599353069893743,infty[ MeV (fixed,scale=5993530.69893743)
  Intensity4 ...............: 0.804148948547538 +/- 0.258602777342412 [7.85901111449516e-16,infty[ ph/cm2/s/MeV (free,scale=7.85901111449516e-06,gradient)
  Energy5 ..................: 10252479.454662 [0.0010252479454662,infty[ MeV (fixed,scale=10252479.454662)
  Intensity5 ...............: 0.947268446501844 +/- 0.33090320971996 [2.4199620497598e-16,infty[ ph/cm2/s/MeV (free,scale=2.4199620497598e-06,gradient)
  Energy6 ..................: 17537798.7113509 [0.00175377987113509,infty[ MeV (fixed,scale=17537798.7113509)
  Intensity6 ...............: 0.510924531144191 +/- 0.330349647333246 [7.45159440158635e-17,infty[ ph/cm2/s/MeV (free,scale=7.45159440158635e-07,gradient)
  Energy7 ..................: 30000000 [0.003,infty[ MeV (fixed,scale=30000000)
  Intensity7 ...............: 1.98644760587573 +/- 1.95978734590705 [2.29450949990163e-17,infty[ ph/cm2/s/MeV (free,scale=2.29450949990163e-07,gradient)
 Number of temporal par's ..: 1
  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
=== GCTAModelBackground ===
 Name ......................: Background_023592
 Instruments ...............: HESSOnOff
 Observation identifiers ...: 023592
 Model type ................: CTABackground
 Model components ..........: "Multiplicative" * "NodeFunction" * "Constant"
 Number of parameters ......: 20
 Number of spatial par's ...: 3
  1:Normalization ..........: 1 [0.001,1000]  (fixed,scale=1,gradient)
  2:Grad_DETX ..............: 0.0290076925442388 +/- 0 deg^-1 (free,scale=1,gradient)
  2:Grad_DETY ..............: 0.165459453396994 +/- 0 deg^-1 (free,scale=1,gradient)
 Number of spectral par's ..: 16
  Energy0 ..................: 700000 [7e-05,infty[ MeV (fixed,scale=700000)
  Intensity0 ...............: 1.15673248048702 +/- 0.961419402688181 [9.64800822522881e-14,infty[ ph/cm2/s/MeV (free,scale=0.000964800822522881,gradient)
  Energy1 ..................: 1197413.67463685 [0.000119741367463685,infty[ MeV (fixed,scale=1197413.67463685)
  Intensity1 ...............: 0.995573238520915 +/- 0.167640057372438 [3.05783271724063e-14,infty[ ph/cm2/s/MeV (free,scale=0.000305783271724063,gradient)
  Energy2 ..................: 2048285.01172474 [0.000204828501172474,infty[ MeV (fixed,scale=2048285.01172474)
  Intensity2 ...............: 0.97430603880028 +/- 0.15455196754381 [9.6914727976462e-15,infty[ ph/cm2/s/MeV (free,scale=9.6914727976462e-05,gradient)
  Energy3 ..................: 3503777.83227556 [0.000350377783227556,infty[ MeV (fixed,scale=3503777.83227556)
  Intensity3 ...............: 0.748497898296835 +/- 0.15953787473352 [3.07160834724385e-15,infty[ ph/cm2/s/MeV (free,scale=3.07160834724385e-05,gradient)
  Energy4 ..................: 5993530.69893743 [0.000599353069893743,infty[ MeV (fixed,scale=5993530.69893743)
  Intensity4 ...............: 0.725418936580643 +/- 0.203356640612422 [9.735133179293e-16,infty[ ph/cm2/s/MeV (free,scale=9.735133179293e-06,gradient)
  Energy5 ..................: 10252479.454662 [0.0010252479454662,infty[ MeV (fixed,scale=10252479.454662)
  Intensity5 ...............: 1.11101986414082 +/- 0.366485282801602 [3.08544603688196e-16,infty[ ph/cm2/s/MeV (free,scale=3.08544603688196e-06,gradient)
  Energy6 ..................: 17537798.7113509 [0.00175377987113509,infty[ MeV (fixed,scale=17537798.7113509)
  Intensity6 ...............: 1.39697833190385 +/- 0.57861299733536 [9.77899025229569e-17,infty[ ph/cm2/s/MeV (free,scale=9.77899025229569e-07,gradient)
  Energy7 ..................: 30000000 [0.003,infty[ MeV (fixed,scale=30000000)
  Intensity7 ...............: 1.25607862208781 +/- 0.917860648895514 [3.09934606573553e-17,infty[ ph/cm2/s/MeV (free,scale=3.09934606573553e-07,gradient)
 Number of temporal par's ..: 1
  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)

Results for the spectrum of the Crab nebula are once again consistent with those obtained before. Note how the spatial parameters of the background model, although formally free, have not been optimised (their errors are 0), since we are running a spectral analysis only.

Since we have provided a background model to csphagen the statistic was set by default to CSTAT, i.e. the background model was also used in the fit. You may also specify as statistic in ctlike WSTAT. In this case the spectral model for the background would be derived from the data during the fit, and the impact of the model used above would extend only to the computation of the background normalisation factor between On and Off regions. Using an a-priori background model on the other hand can be advantageous to avoid intrinsic biases of WSTAT when the counting statistics are limited.

Let’s check the residuals for the first of the 4 observations.

residuals = 'residuals_classical_bkgmod.fits'
res = cscripts.csresspec(like.obs())
res['components'] = True
res['algorithm']  = 'SIGNIFICANCE'
res['outfile']    = residuals
res['stack']      = False
res.execute()
from show_residuals import plot_residuals
plot_residuals(residuals,'',0)

As expected the spectral residuals are not as good as in the first analysis in which the model was derived from the data during the fit itself, but the models still satisfactorily reproduce the data. Of course the total background coincides with the background model of the specific observation considered.