Bug #2312

csspec gives upper limits only for statistic=wstat and absent bkg model

Added by Specovius Andreas about 6 years ago. Updated almost 6 years ago.

Status:ClosedStart date:02/12/2018
Priority:HighDue date:
Assigned To:Knödlseder Jürgen% Done:

100%

Category:-
Target version:1.6.0
Duration:

Description

I used csspec to generate a spectrum for a classical on/off analysis (following Luigis tutorial on CasA).
In principle when using statistic=wstat no background information has to be provided.
Hence, in principle, I would assume that it should be possible to use models which only contain a source model, no background model.

Indeed this works fine for ctlike.
For csspec things are different:
If no background model (only src model) is provided when using csspec (statistic=wstat) the resulting spectrum consists of only upper limits.
If a background model is provided (even with senseless parameters) the spectrum shows nice spectral points.

As statistic=wstat implies no need of a description of the background I would assume that csspec should also work for models that do not contain a bkg component?

spectrum_issue.py Magnifier (3.39 KB) Specovius Andreas, 06/15/2018 02:57 PM

Recurrence

No recurrence.

History

#1 Updated by Tibaldo Luigi almost 6 years ago

  • Assigned To set to Knödlseder Jürgen

Jürgen, would you be willing to look into this? I just found it is still a problem ...

#2 Updated by Knödlseder Jürgen almost 6 years ago

  • Priority changed from Normal to High
  • Target version set to 1.6.0

#3 Updated by Knödlseder Jürgen almost 6 years ago

  • Status changed from New to In Progress
  • % Done changed from 0 to 10

I was not able to reproduce the problem with the current software version.

I simulated an event file using ctobssim:

$ 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] 
Lower energy limit (TeV) [0.1] 
Upper energy limit (TeV) [100.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]

I then copied over the input model and produced an XML file with a On/Off background model and another background file without a background model:

$ cp $CTOOLS/share/models/crab.xml crab_onoff.xml
nano crab_onoff.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="CTAOnOff">
    <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>


$ cp crab_onoff.xml crab_onoff_nobgd.xml
$ nano crab_onoff_nobgd.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_library>

Then I ran csphagen:
$ csphagen
Input event list or observation definition XML file [obs.xml] events.fits
Calibration database [prod2] 
Instrument response function [South_0.5h] 
Input model definition XML file (if NONE, use point source) [NONE] 
Algorithm for defining energy bins (FILE|LIN|LOG) [LOG] 
Start value for first energy bin in TeV [0.1] 
Stop value for last energy bin in TeV [100.0] 
Number of energy bins [120] 20
Stack multiple observations into single PHA, ARF and RMF files? [no] 
Output observation definition XML file [onoff_obs.xml] 
Method for background estimation (REFLECTED|CUSTOM) [REFLECTED] 
Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] 
Right Ascension of source region centre (deg) (0-360) [83.63] 
Declination of source region centre (deg) (-90-90) [22.01] 
Radius of source region circle (deg) (0-180) [0.2]

And finally csspec, one with the XML file with background model, and once without the background model:
$ csspec debug=yes statistic=wstat
Input event list, counts cube, or observation definition XML file [events.fits] onoff_obs.xml
Input model definition XML file [$CTOOLS/share/models/crab.xml] crab_onoff.xml
Source name [Crab] 
Spectrum generation method (SLICE|NODES|AUTO) [AUTO] 
Algorithm for defining energy bins (FILE|LIN|LOG) [LOG] 
Start value for first energy bin in TeV [0.1] 
Stop value for last energy bin in TeV [100.0] 
Number of energy bins (1-200) [20] 
Output spectrum file [spectrum.fits] 
2018-06-15T08:44:45: +===================+
2018-06-15T08:44:45: | Generate spectrum |
2018-06-15T08:44:45: +===================+
2018-06-15T08:44:45: === GEbounds ===
2018-06-15T08:44:45:  Number of intervals .......: 20
2018-06-15T08:44:45:  Energy range ..............: 100 GeV - 100 TeV
2018-06-15T08:44:45:  Bin 1 .....................: 1.221688e-10 +/- 5.373786e-12 [< 1.759066e-10] erg/cm2/s (TS = 1069.409)
2018-06-15T08:44:45:  Bin 2 .....................: 1.067982e-10 +/- 4.919193e-12 [< 1.559900e-10] erg/cm2/s (TS = 1064.433)
2018-06-15T08:44:45:  Bin 3 .....................: 1.011262e-10 +/- 4.647995e-12 [< 1.476061e-10] erg/cm2/s (TS = 1131.287)
2018-06-15T08:44:45:  Bin 4 .....................: 8.776835e-11 +/- 4.286944e-12 [< 1.306378e-10] erg/cm2/s (TS = 1046.271)
2018-06-15T08:44:45:  Bin 5 .....................: 7.005584e-11 +/- 3.853282e-12 [< 1.085886e-10] erg/cm2/s (TS = 853.251)
2018-06-15T08:44:45:  Bin 6 .....................: 5.759921e-11 +/- 3.576058e-12 [< 9.335975e-11] erg/cm2/s (TS = 684.670)
2018-06-15T08:44:45:  Bin 7 .....................: 4.871038e-11 +/- 3.294468e-12 [< 8.165503e-11] erg/cm2/s (TS = 598.169)
2018-06-15T08:44:45:  Bin 8 .....................: 4.619379e-11 +/- 3.306368e-12 [< 7.925744e-11] erg/cm2/s (TS = 530.022)
2018-06-15T08:44:45:  Bin 9 .....................: 3.258974e-11 +/- 2.893858e-12 [< 6.152829e-11] erg/cm2/s (TS = 355.286)
2018-06-15T08:44:45:  Bin 10 ....................: 3.437429e-11 +/- 3.165327e-12 [< 6.602754e-11] erg/cm2/s (TS = 334.213)
2018-06-15T08:44:45:  Bin 11 ....................: 1.745382e-11 +/- 2.478647e-12 [< 4.224027e-11] erg/cm2/s (TS = 129.701)
2018-06-15T08:44:45:  Bin 12 ....................: 1.818815e-11 +/- 2.835562e-12 [< 4.654376e-11] erg/cm2/s (TS = 106.882)
2018-06-15T08:44:45:  Bin 13 ....................: 2.153167e-11 +/- 3.448843e-12 [< 5.602009e-11] erg/cm2/s (TS = 125.532)
2018-06-15T08:44:45:  Bin 14 ....................: 1.819759e-11 +/- 3.681961e-12 [< 5.501719e-11] erg/cm2/s (TS = 72.438)
2018-06-15T08:44:45:  Bin 15 ....................: 1.406384e-11 +/- 3.759826e-12 [< 5.166209e-11] erg/cm2/s (TS = 45.063)
2018-06-15T08:44:45:  Bin 16 ....................: 8.405661e-12 +/- 3.432607e-12 [< 3.710673e-11] erg/cm2/s (TS = 19.313)
2018-06-15T08:44:45:  Bin 17 ....................: 1.760389e-11 +/- 5.869693e-12 [< 7.527111e-11] erg/cm2/s (TS = 28.969)
2018-06-15T08:44:45:  Bin 18 ....................: 1.095723e-11 +/- 5.480227e-12 [< 5.062560e-11] erg/cm2/s (TS = 12.875)
2018-06-15T08:44:45:  Bin 19 ....................: No event in this bin. Likelihood is zero. Bin is skipped.
2018-06-15T08:44:45:  Bin 20 ....................: No event in this bin. Likelihood is zero. Bin is skipped.


$ csspec debug=yes statistic=wstat
Input event list, counts cube, or observation definition XML file [crab_onoff_nobgd.xml] onoff_obs.xml
Input model definition XML file [crab_onoff.xml] crab_onoff_nobgd.xml
Source name [Crab] 
Spectrum generation method (SLICE|NODES|AUTO) [AUTO] 
Algorithm for defining energy bins (FILE|LIN|LOG) [LOG] 
Start value for first energy bin in TeV [0.1] 
Stop value for last energy bin in TeV [100.0] 
Number of energy bins (1-200) [20] 
Output spectrum file [spectrum.fits]
2018-06-15T08:45:39: +===================+
2018-06-15T08:45:39: | Generate spectrum |
2018-06-15T08:45:39: +===================+
2018-06-15T08:45:39: === GEbounds ===
2018-06-15T08:45:39:  Number of intervals .......: 20
2018-06-15T08:45:39:  Energy range ..............: 100 GeV - 100 TeV
2018-06-15T08:45:39:  Bin 1 .....................: 1.221688e-10 +/- 5.373786e-12 [< 1.329479e-10] erg/cm2/s
2018-06-15T08:45:39:  Bin 2 .....................: 1.067982e-10 +/- 4.919193e-12 [< 1.166917e-10] erg/cm2/s
2018-06-15T08:45:39:  Bin 3 .....................: 1.011262e-10 +/- 4.647995e-12 [< 1.104821e-10] erg/cm2/s
2018-06-15T08:45:39:  Bin 4 .....................: 8.776835e-11 +/- 4.286944e-12 [< 9.641762e-11] erg/cm2/s
2018-06-15T08:45:39:  Bin 5 .....................: 7.005584e-11 +/- 3.853282e-12 [< 7.786330e-11] erg/cm2/s
2018-06-15T08:45:39:  Bin 6 .....................: 5.759921e-11 +/- 3.576058e-12 [< 6.487794e-11] erg/cm2/s
2018-06-15T08:45:39:  Bin 7 .....................: 4.871038e-11 +/- 3.294468e-12 [< 5.544242e-11] erg/cm2/s
2018-06-15T08:45:39:  Bin 8 .....................: 4.619379e-11 +/- 3.306368e-12 [< 5.296566e-11] erg/cm2/s
2018-06-15T08:45:39:  Bin 9 .....................: 3.258974e-11 +/- 2.893858e-12 [< 3.858303e-11] erg/cm2/s
2018-06-15T08:45:39:  Bin 10 ....................: 3.437429e-11 +/- 3.165327e-12 [< 4.094432e-11] erg/cm2/s
2018-06-15T08:45:39:  Bin 11 ....................: 1.745382e-11 +/- 2.478647e-12 [< 2.274533e-11] erg/cm2/s
2018-06-15T08:45:39:  Bin 12 ....................: 1.818815e-11 +/- 2.835562e-12 [< 2.428989e-11] erg/cm2/s
2018-06-15T08:45:39:  Bin 13 ....................: 2.153167e-11 +/- 3.448843e-12 [< 2.901399e-11] erg/cm2/s
2018-06-15T08:45:39:  Bin 14 ....................: 1.819759e-11 +/- 3.681961e-12 [< 2.638040e-11] erg/cm2/s
2018-06-15T08:45:39:  Bin 15 ....................: 1.406384e-11 +/- 3.759826e-12 [< 2.277227e-11] erg/cm2/s
2018-06-15T08:45:39:  Bin 16 ....................: 8.405661e-12 +/- 3.432607e-12 [< 1.703334e-11] erg/cm2/s
2018-06-15T08:45:39:  Bin 17 ....................: 1.760389e-11 +/- 5.869693e-12 [< 3.173481e-11] erg/cm2/s
2018-06-15T08:45:39:  Bin 18 ....................: 1.095723e-11 +/- 5.480227e-12 [< 2.545617e-11] erg/cm2/s
2018-06-15T08:45:39:  Bin 19 ....................: 6.484044e-12 +/- 0.000000e+00 [< 1.386468e-11] erg/cm2/s
2018-06-15T08:45:39:  Bin 20 ....................: 5.493459e-12 +/- 0.000000e+00 [< 1.581612e-11] erg/cm2/s

There are differences (no TS evaluation in second case, no computations for bin 19 and 20 in first case), but in both cases csspec has made fits and computed upper limits.

Can you post a specific example where the code is not working in that way?

#4 Updated by Knödlseder Jürgen almost 6 years ago

  • % Done changed from 10 to 20

The issue was that in the absence of a background model, the model used for TS computation is in fact empty, and therefore has no parameters, which led to a likelihood value of 0. I corrected GObservations::likelihood::eval() and CTAOnOffObservation::likelihood_wstat() so that they also return now meaningful likelihood values in the absence of model parameters.

In addition, an additional check was needed in csspec to not return any results if there are no events in the energy bin. This now gives the following result for the XML model without any background, which is identical to the result obtained with a background component:

2018-06-15T12:36:57: +===================+
2018-06-15T12:36:57: | Generate spectrum |
2018-06-15T12:36:57: +===================+
2018-06-15T12:36:57: === GEbounds ===
2018-06-15T12:36:57:  Number of intervals .......: 20
2018-06-15T12:36:57:  Energy range ..............: 100 GeV - 100 TeV
2018-06-15T12:36:57:  Bin 1 .....................: 1.221688e-10 +/- 5.373786e-12 [< 1.759066e-10] erg/cm2/s (TS = 1069.409)
2018-06-15T12:36:57:  Bin 2 .....................: 1.067982e-10 +/- 4.919193e-12 [< 1.559900e-10] erg/cm2/s (TS = 1064.433)
2018-06-15T12:36:57:  Bin 3 .....................: 1.011262e-10 +/- 4.647995e-12 [< 1.476061e-10] erg/cm2/s (TS = 1131.287)
2018-06-15T12:36:57:  Bin 4 .....................: 8.776835e-11 +/- 4.286944e-12 [< 1.306378e-10] erg/cm2/s (TS = 1046.271)
2018-06-15T12:36:57:  Bin 5 .....................: 7.005584e-11 +/- 3.853282e-12 [< 1.085886e-10] erg/cm2/s (TS = 853.251)
2018-06-15T12:36:57:  Bin 6 .....................: 5.759921e-11 +/- 3.576058e-12 [< 9.335975e-11] erg/cm2/s (TS = 684.670)
2018-06-15T12:36:57:  Bin 7 .....................: 4.871038e-11 +/- 3.294468e-12 [< 8.165503e-11] erg/cm2/s (TS = 598.169)
2018-06-15T12:36:57:  Bin 8 .....................: 4.619379e-11 +/- 3.306368e-12 [< 7.925744e-11] erg/cm2/s (TS = 530.022)
2018-06-15T12:36:57:  Bin 9 .....................: 3.258974e-11 +/- 2.893858e-12 [< 6.152829e-11] erg/cm2/s (TS = 355.286)
2018-06-15T12:36:57:  Bin 10 ....................: 3.437429e-11 +/- 3.165327e-12 [< 6.602754e-11] erg/cm2/s (TS = 334.213)
2018-06-15T12:36:57:  Bin 11 ....................: 1.745382e-11 +/- 2.478647e-12 [< 4.224027e-11] erg/cm2/s (TS = 129.701)
2018-06-15T12:36:57:  Bin 12 ....................: 1.818815e-11 +/- 2.835562e-12 [< 4.654376e-11] erg/cm2/s (TS = 106.882)
2018-06-15T12:36:57:  Bin 13 ....................: 2.153167e-11 +/- 3.448843e-12 [< 5.602009e-11] erg/cm2/s (TS = 125.532)
2018-06-15T12:36:57:  Bin 14 ....................: 1.819759e-11 +/- 3.681961e-12 [< 5.501719e-11] erg/cm2/s (TS = 72.438)
2018-06-15T12:36:57:  Bin 15 ....................: 1.406384e-11 +/- 3.759826e-12 [< 5.166209e-11] erg/cm2/s (TS = 45.063)
2018-06-15T12:36:57:  Bin 16 ....................: 8.405661e-12 +/- 3.432607e-12 [< 3.710673e-11] erg/cm2/s (TS = 19.313)
2018-06-15T12:36:57:  Bin 17 ....................: 1.760389e-11 +/- 5.869693e-12 [< 7.527111e-11] erg/cm2/s (TS = 28.969)
2018-06-15T12:36:57:  Bin 18 ....................: 1.095723e-11 +/- 5.480227e-12 [< 5.062560e-11] erg/cm2/s (TS = 12.875)
2018-06-15T12:36:57:  Bin 19 ....................: No event in this bin. Likelihood is zero. Bin is skipped.
2018-06-15T12:36:57:  Bin 20 ....................: No event in this bin. Likelihood is zero. Bin is skipped.

#5 Updated by Specovius Andreas almost 6 years ago

Knödlseder Jürgen wrote:

The issue was that in the absence of a background model, the model used for TS computation is in fact empty, and therefore has no parameters, which led to a likelihood value of 0. I corrected GObservations::likelihood::eval() and CTAOnOffObservation::likelihood_wstat() so that they also return now meaningful likelihood values in the absence of model parameters.

In addition, an additional check was needed in csspec to not return any results if there are no events in the energy bin. This now gives the following result for the XML model without any background, which is identical to the result obtained with a background component:
[...]

Good to hear!
I just coded an analysis script that reproduces the upper limits issue using CasA 1dc data but when I wanted to upload it here I saw your new post.
Do you still need these lines?

#6 Updated by Knödlseder Jürgen almost 6 years ago

  • Status changed from In Progress to Feedback
  • % Done changed from 10 to 80

Well, it’s not clear to me whether your issue is fixed by my modifications. The code is now in the devel branch. Could you please check whether the actual code works for you?

#7 Updated by Specovius Andreas almost 6 years ago

I ran the same script, which behaved in the described way, again using the latest ctools&gammalib code from the devel branch.
As what I expected from your notes above, the script now produces exactly the same spectrum for a model including and for a model excluding a description of the background.

It seems you solved the issue! :)

#8 Updated by Knödlseder Jürgen almost 6 years ago

  • Status changed from Feedback to Closed
  • % Done changed from 80 to 100

Also available in: Atom PDF