Bug #1357
ctobssim does not simulate background for high upper energy limit.
| Status: | Closed | Start date: | 11/05/2014 | |
|---|---|---|---|---|
| Priority: | Urgent | Due date: | ||
| Assigned To: |  Knödlseder Jürgen | % Done: | 100% | |
| Category: | - | |||
| Target version: | 00-08-00 | |||
| Duration: |
Description
Running ctobssim for the 0.01-200 TeV energy range gave zero background events:
2014-11-05T13:33:52: Simulation area ...........: 1.9635e+11 cm2 2014-11-05T13:33:52: Simulation cone ...........: RA=83.63 deg, Dec=22.01 deg, r=10.5 deg 2014-11-05T13:33:52: Time interval .............: 0 - 3600 s 2014-11-05T13:33:52: Photon energy range .......: 10 GeV - 200 TeV 2014-11-05T13:33:52: Event energy range ........: 10 GeV - 200 TeV 2014-11-05T13:34:03: MC source photons .........: 12538094 [Crab] 2014-11-05T13:34:03: MC source events ..........: 5072 [Crab] 2014-11-05T13:34:03: MC source events ..........: 5072 (all source models) 2014-11-05T13:34:04: MC background events ......: 0 2014-11-05T13:34:04: MC events .................: 5072 (all models)
while 0.01-100 TeV gives background events:
2014-11-05T13:32:53: Simulation area ...........: 1.9635e+11 cm2 2014-11-05T13:32:53: Simulation cone ...........: RA=83.63 deg, Dec=22.01 deg, r=10.5 deg 2014-11-05T13:32:53: Time interval .............: 0 - 3600 s 2014-11-05T13:32:53: Photon energy range .......: 10 GeV - 100 TeV 2014-11-05T13:32:53: Event energy range ........: 10 GeV - 100 TeV 2014-11-05T13:33:04: MC source photons .........: 12538086 [Crab] 2014-11-05T13:33:04: MC source events ..........: 5083 [Crab] 2014-11-05T13:33:04: MC source events ..........: 5083 (all source models) 2014-11-05T13:33:07: MC background events ......: 179583 2014-11-05T13:33:07: MC events .................: 184666 (all models)
Recurrence
No recurrence.
History
#1
Updated by Knödlseder Jürgen about 10 years ago
- Status changed from New to In Progress
- % Done changed from 0 to 10
The problem seems to be related to some discontinuity in the background spectrum:
set_flux_cache: fmin=3.82746e-11 fmax=2.7339e-11 emin=1.34896e+08 emax=1.41254e+08 epivot=1.38038e+08 gamma=-7.3064 /...=1.18322 prefactor=3.23479e-11 flux=0.000206356 set_flux_cache: fmin=2.7339e-11 fmax=1.64034e-11 emin=1.41254e+08 emax=1.47911e+08 epivot=1.44544e+08 gamma=-11.0924 /...=1.29099 prefactor=2.11767e-11 flux=0.000142235 set_flux_cache: fmin=1.64034e-11 fmax=5.4678e-12 emin=1.47911e+08 emax=1.54882e+08 epivot=1.51356e+08 gamma=-23.8561 /...=1.73205 prefactor=9.47051e-12 flux=6.91011e-05 set_flux_cache: fmin=5.4678e-12 fmax=1.93588e-14 emin=1.54882e+08 emax=1.62181e+08 epivot=1.58489e+08 gamma=-122.547 /...=16.8061 prefactor=3.25346e-13 flux=nan set_flux_cache: fmin=1.93588e-14 fmax=5.80765e-14 emin=1.62181e+08 emax=1.69824e+08 epivot=1.65959e+08 gamma=23.8561 /...=0.57735 prefactor=3.35305e-14 flux=2.70484e-07 set_flux_cache: fmin=5.80765e-14 fmax=9.67942e-14 emin=1.69824e+08 emax=1.77828e+08 epivot=1.7378e+08 gamma=11.0924 /...=0.774597 prefactor=7.49765e-14 flux=6.0781e-07 set_flux_cache: fmin=9.67942e-14 fmax=1.35512e-13 emin=1.77828e+08 emax=1.86209e+08 epivot=1.8197e+08 gamma=7.3064 /...=0.845154 prefactor=1.14528e-13 flux=9.65615e-07 set_flux_cache: fmin=1.35512e-13 fmax=1.7423e-13 emin=1.86209e+08 emax=1.94984e+08 epivot=1.90546e+08 gamma=5.45722 /...=0.881917 prefactor=1.53656e-13 flux=1.3533e-06
Note that the background spectrum covers 21 bins running from 0.0126 TeV to 199.53 TeV. Not clear so far why a discontinuity exists at 155-162 TeV. Energy bins at the end are:
50.12- 79.43 79.43-125.89 125.89-199.53
#2
Updated by Knödlseder Jürgen about 10 years ago
- % Done changed from 10 to 100
It turned out that the problem was due to a bad formulation of the power law integration in gammalib::plaw_photon_flux. The formula made use of pow(energy, -gamma) which for a large negative gamma led to an overflow. Reformulating the integration by substituting x=E/E0 avoids the additional pow call, speeding up at the same time the computations.
The same problem was inherent to gammalib::plaw_energy_flux and has been corrected.
Unit tests have been added to check from now on the computations.
#3
Updated by Knödlseder Jürgen about 10 years ago
- Status changed from In Progress to Closed