Bug #3624

bug on csspec with Composite model

Added by Burtovoi Aleksandr over 3 years ago. Updated over 3 years ago.

Status:ClosedStart date:04/27/2021
Priority:NormalDue date:
Assigned To:Knödlseder Jürgen% Done:

100%

Category:-
Target version:1.7.4
Duration:

Description

csspec task, applied to the source with a Composite spectral model, returns a 0.0 value of the flux upper limit in all energy bins.
ctools version: 1.7.3
Attached files represent the test-analysis

sim_model.xml Magnifier - input xml-file for ctobssim (1.95 KB) Burtovoi Aleksandr, 04/27/2021 03:42 PM

ctobssim.log - log of the test-simulations (13.1 KB) Burtovoi Aleksandr, 04/27/2021 03:42 PM

ctlike_outmodel.xml Magnifier - xml-file provided as an output of ctlike (1.88 KB) Burtovoi Aleksandr, 04/27/2021 03:42 PM

spec.fits - output of csspec (11.3 KB) Burtovoi Aleksandr, 04/27/2021 03:42 PM

csspec.log - log of csspec (8.64 KB) Burtovoi Aleksandr, 04/27/2021 03:42 PM

spec_aleksandr.png (29.8 KB) Knödlseder Jürgen, 04/28/2021 10:42 AM

spec_old.png (30.6 KB) Knödlseder Jürgen, 04/28/2021 10:42 AM

spec_kludge.png (30.7 KB) Knödlseder Jürgen, 04/28/2021 11:28 AM

Spec_aleksandr Spec_old Spec_kludge

Recurrence

No recurrence.

Related issues

Related to GammaLib - Action #3625: Implement a more efficient Monte-Carlo sampler for expone... Closed 04/27/2021
Related to GammaLib - Feature #3626: Add GModelSpectralFunc model constructor Closed 04/28/2021

History

#1 Updated by Knödlseder Jürgen over 3 years ago

  • Status changed from New to In Progress
  • Assigned To set to Knödlseder Jürgen

I started investigating the problem. As a first side note I note the long time spent by ctobssim in simulating the first two energy ranges while all subsequent energy ranges basically need no time:

2021-04-23T14:32:50: === CTA observation ===
2021-04-23T14:32:50:  Simulation cone ...........: RA=128.836 deg, Dec=-45.1764 deg, radius=1.5 deg
2021-04-23T14:32:50:  Time interval .............: 0 - 180000 s
2021-04-23T14:32:50:  Photon energy range .......: 30 GeV - 72.3259983091833 GeV
2021-04-23T14:32:50:  Event energy range ........: 30 GeV - 72.3259983091833 GeV
2021-04-23T14:32:50:   Simulation area ..........: 1.73683e+09 cm2
2021-04-23T14:32:50:   Use model ................: PSR_J0835-4510
2021-04-23T14:32:50:   Normalization ............: 1 [PSR_J0835-4510]
2021-04-23T14:32:50:   Flux .....................: 1.40407e-10 [PSR_J0835-4510] photons/cm2/s
2021-04-23T14:32:50:   Normalized flux ..........: 1.40407e-10 [PSR_J0835-4510] photons/cm2/s
2021-04-23T14:32:50:   Photon rate ..............: 0.243863 photons/s [PSR_J0835-4510]
2021-04-23T14:51:33:   MC source photons ........: 43656 [PSR_J0835-4510]
2021-04-23T14:51:33:   MC source events .........: 8343 [PSR_J0835-4510]
2021-04-23T14:51:33:   MC source events .........: 8343 (all source models)
2021-04-23T14:51:33:  Photon energy range .......: 72.3259983091833 GeV - 174.368334380666 GeV
2021-04-23T14:51:33:  Event energy range ........: 72.3259983091833 GeV - 174.368334380666 GeV
2021-04-23T14:51:33:   Simulation area ..........: 4.85741e+09 cm2
2021-04-23T14:51:33:   Use model ................: PSR_J0835-4510
2021-04-23T14:51:33:   Normalization ............: 1 [PSR_J0835-4510]
2021-04-23T14:51:33:   Flux .....................: 1.95471e-12 [PSR_J0835-4510] photons/cm2/s
2021-04-23T14:51:33:   Normalized flux ..........: 1.95471e-12 [PSR_J0835-4510] photons/cm2/s
2021-04-23T14:51:33:   Photon rate ..............: 0.00949484 photons/s [PSR_J0835-4510]
2021-04-23T15:45:59:   MC source photons ........: 1678 [PSR_J0835-4510]
2021-04-23T15:45:59:   MC source events .........: 425 [PSR_J0835-4510]
2021-04-23T15:45:59:   MC source events .........: 425 (all source models)
2021-04-23T15:45:59:  Photon energy range .......: 174.368334380666 GeV - 420.37879525304 GeV
2021-04-23T15:45:59:  Event energy range ........: 174.368334380666 GeV - 420.37879525304 GeV
2021-04-23T15:45:59:   Simulation area ..........: 9.34853e+09 cm2
2021-04-23T15:45:59:   Use model ................: PSR_J0835-4510
2021-04-23T15:45:59:   Normalization ............: 1 [PSR_J0835-4510]
2021-04-23T15:45:59:   Flux .....................: 7.28198e-13 [PSR_J0835-4510] photons/cm2/s
2021-04-23T15:45:59:   Normalized flux ..........: 7.28198e-13 [PSR_J0835-4510] photons/cm2/s
2021-04-23T15:45:59:   Photon rate ..............: 0.00680758 photons/s [PSR_J0835-4510]
2021-04-23T15:45:59:   MC source photons ........: 1195 [PSR_J0835-4510]
2021-04-23T15:45:59:   MC source events .........: 417 [PSR_J0835-4510]
2021-04-23T15:45:59:   MC source events .........: 417 (all source models)

#2 Updated by Knödlseder Jürgen over 3 years ago

  • Related to Action #3625: Implement a more efficient Monte-Carlo sampler for exponentially and super-exponentially cut-off power laws added

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

  • % Done changed from 0 to 10

The reason for the long simulation times is the first super-exponential cut-off model in the composite spectrum

<spectrum type="SuperExponentialCutoffPowerLaw">
    <parameter name="Prefactor" value="17.64" scale="1e-09" min="1e-07" max="1000" free="0" />
    <parameter name="Index1" value="0.913" scale="-1" min="0" max="5" free="0" />
    <parameter name="CutoffEnergy" value="0.143" scale="1000" min="0.01" max="1000" free="0" />
    <parameter name="Index2" value="0.439" scale="1" min="0.1" max="5" free="0" />
    <parameter name="PivotEnergy" value="1" scale="1000" min="0.01" max="1000" free="0" />
</spectrum>
that has a cut-off energy of 143 MeV. The Monte Carlo sampling in GModelSpectralSuperExpPlaw::mc is done using a rejection method that draws random energies for a power law and compares then the super-exponential cut-off model to the power law to decide whether to keep or to reject the energy. Since for the first energy bins the sampling is done deeply in the cut-off of the model, the energy is very often rejected before being accepted. Below some debugging of the Monte Carlo method, where samples indicates the number of random energies that are drawn before one is accepted. This number can go up to millions!
GModelSpectralSuperExpPlaw::mc(30 GeV,72.3259983091833 GeV,-315576065.450442 s (TT)): energy=31.8619801747531 GeV samples=274597
GModelSpectralSuperExpPlaw::mc(30 GeV,72.3259983091833 GeV,-315576049.460237 s (TT)): energy=31.8702496902379 GeV samples=453435
GModelSpectralSuperExpPlaw::mc(30 GeV,72.3259983091833 GeV,-315576046.466759 s (TT)): energy=32.8847377698187 GeV samples=17424
GModelSpectralSuperExpPlaw::mc(30 GeV,72.3259983091833 GeV,-315576043.428479 s (TT)): energy=31.0090430106833 GeV samples=49213
GModelSpectralSuperExpPlaw::mc(30 GeV,72.3259983091833 GeV,-315576041.090847 s (TT)): energy=47.7142965496299 GeV samples=80922
GModelSpectralSuperExpPlaw::mc(30 GeV,72.3259983091833 GeV,-315576039.582654 s (TT)): energy=32.0036966875382 GeV samples=429141
GModelSpectralSuperExpPlaw::mc(30 GeV,72.3259983091833 GeV,-315576039.367193 s (TT)): energy=30.836690685018 GeV samples=244683
GModelSpectralSuperExpPlaw::mc(30 GeV,72.3259983091833 GeV,-315576039.138902 s (TT)): energy=37.6197471157148 GeV samples=1025250
I added issue #3625 for the implementation of a more efficient Monte-Carlo sampler that takes care of this problem.

#4 Updated by Knödlseder Jürgen over 3 years ago

Fixing issue #3625 reduces the event simulation duration from 2254.24 seconds to 42.49 seconds on my Mac OS power book.

#5 Updated by Knödlseder Jürgen over 3 years ago

  • % Done changed from 10 to 20

Investigating csspec here a first remark of something that probably does not work as expected for composite spectral models. When running the script with chatter=4 the following is written in the log file:

2021-04-28T06:41:26: +=========================+
2021-04-28T06:41:26: | Adjust model parameters |
2021-04-28T06:41:26: +=========================+
2021-04-28T06:41:26: === PSR_J0835-4510 ===
2021-04-28T06:41:26:  Fixing "1:Prefactor" 
2021-04-28T06:41:26:  Fixing "1:Index" 
2021-04-28T06:41:26:  Fixing "2:Prefactor" 
2021-04-28T06:41:26:  Fixing "2:Index" 
2021-04-28T06:41:26:  Fixing "2:CutoffEnergy" 
2021-04-28T06:41:26:  Freeing "1:Prefactor" 
2021-04-28T06:41:26: === CTABackgroundModel ===
Obviously, only the Prefactor of the first spectral model is free’d. This allows generation of a spectrum, but is rather arbitrary. Of course, freeing both Prefactor parameters would probably lead to some degeneracy, but it’s not clear whether freeing always the first one is the right choice. Ideally, both pre factors should be combined in a single one, yet this functionality is not yet available in GammaLib. Here the code that fixes the model parameters:
# Fix all model parameters
for par in model:
    if par.is_free():
        self._log_string(gammalib.EXPLICIT, ' Fixing "'+par.name()+'"')
        par.fix()

# Free the normalisation parameter which is assumed to be the first spectral parameter
normpar = model.spectral()[0]
if normpar.is_fixed():
    self._log_string(gammalib.EXPLICIT, ' Freeing "'+normpar.name()+'"')
    normpar.free()
The strong assumption here is that the relevant parameter to free is the first parameter of the spectral model.

Since all parameters of the source of interest are fixed in any case, expect of the normalisation parameter, the best thing to do would be to convert the spectral model of the source into a file function, and leave only the normalisation of the file function as a free parameter. This would cope with all kinds of models, at probably fix also the upper limit computations.

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

  • % Done changed from 20 to 30

For the record, here is what the log file with chatter=4 states about the upper limit calculation. It obviously failed:

2021-04-28T07:04:07: === Computing upper limit for energy bin ===
2021-04-28T07:04:07: Upper limit calculation failed.
2021-04-28T07:04:07:  Bin 1 .....................: 1.427690e-11 +/- 5.416689e-13 erg/cm2/s (TS = 772.289)

#7 Updated by Knödlseder Jürgen over 3 years ago

  • Target version set to 1.7.4

Unfortunately, the GModelSpectralFunc methods that allow for an in-memory construction of a file function only become available in version 2.0.0 of the code, hence I have to implement a kludge for the bugfix-1.7.4 release to do the conversion.

#8 Updated by Knödlseder Jürgen over 3 years ago

  • Related to Feature #3626: Add GModelSpectralFunc model constructor added

#9 Updated by Knödlseder Jürgen over 3 years ago

For reference, here are the spectral results found by Aleksandr

2021-04-26T08:11:02: +===================+
2021-04-26T08:11:02: | Generate spectrum |
2021-04-26T08:11:02: +===================+
2021-04-26T08:11:02: === GEbounds ===
2021-04-26T08:11:02:  Number of intervals .......: 20
2021-04-26T08:11:02:  Energy range ..............: 30 GeV - 199 TeV
2021-04-26T08:11:23:  Bin 1 .....................: 1.301655e-11 +/- 5.383611e-13 erg/cm2/s (TS = 645.591)
2021-04-26T08:11:36:  Bin 2 .....................: 2.462624e-12 +/- 2.973857e-13 erg/cm2/s (TS = 79.929)
2021-04-26T08:11:45:  Bin 3 .....................: 2.969283e-13 +/- 4.089534e-13 erg/cm2/s (TS = 4.707)
2021-04-26T08:11:53:  Bin 4 .....................: 2.379248e-13 +/- 0.000000e+00 erg/cm2/s (TS = 7.175)
2021-04-26T08:11:59:  Bin 5 .....................: 3.165973e-13 +/- 7.452551e-13 erg/cm2/s (TS = 24.758)
2021-04-26T08:12:03:  Bin 6 .....................: 3.663549e-13 +/- 3.277524e-12 erg/cm2/s (TS = 102.696)
2021-04-26T08:12:08:  Bin 7 .....................: 4.842922e-13 +/- 3.771110e-13 erg/cm2/s (TS = 201.624)
2021-04-26T08:12:11:  Bin 8 .....................: 5.477164e-13 +/- 8.313567e-13 erg/cm2/s (TS = 410.724)
2021-04-26T08:12:15:  Bin 9 .....................: 6.135478e-13 +/- 0.000000e+00 erg/cm2/s (TS = 718.972)
2021-04-26T08:12:19:  Bin 10 ....................: 6.862635e-13 +/- 0.000000e+00 erg/cm2/s (TS = 1105.524)
2021-04-26T08:12:24:  Bin 11 ....................: 7.206437e-13 +/- 0.000000e+00 erg/cm2/s (TS = 1529.553)
2021-04-26T08:12:28:  Bin 12 ....................: 7.169484e-13 +/- 8.534060e-13 erg/cm2/s (TS = 1256.404)
2021-04-26T08:12:31:  Bin 13 ....................: 5.713853e-13 +/- 2.095115e-12 erg/cm2/s (TS = 779.375)
2021-04-26T08:12:36:  Bin 14 ....................: 3.621812e-13 +/- 0.000000e+00 erg/cm2/s (TS = 401.523)
2021-04-26T08:12:39:  Bin 15 ....................: 1.631374e-13 +/- 4.508672e-12 erg/cm2/s (TS = 125.031)
2021-04-26T08:12:43:  Bin 16 ....................: 4.491535e-14 +/- 2.430208e-13 erg/cm2/s
2021-04-26T08:12:48:  Bin 17 ....................: 4.284596e-15 +/- 0.000000e+00 erg/cm2/s
2021-04-26T08:12:52:  Bin 18 ....................: 1.133950e-16 +/- 0.000000e+00 erg/cm2/s
2021-04-26T08:12:56:  Bin 19 ....................: 3.557983e-19 +/- 0.000000e+00 erg/cm2/s
2021-04-26T08:13:00:  Bin 20 ....................: 4.073476e-23 +/- 0.000000e+00 erg/cm2/s
and here are the one found in my analysis before implementing the kludge. Note that my results differ slightly from the those found by Aleksandr because I implemented in the mean time the more efficient Monte Carlo sampling of the exponential power laws (see #3625) which will produce a different event sampling for the test case.
2021-04-28T07:51:44: +===================+
2021-04-28T07:51:44: | Generate spectrum |
2021-04-28T07:51:44: +===================+
2021-04-28T07:51:44: === GEbounds ===
2021-04-28T07:51:44:  Number of intervals .......: 20
2021-04-28T07:51:44:  Energy range ..............: 30 GeV - 199 TeV
2021-04-28T07:04:07:  Bin 1 .....................: 1.427690e-11 +/- 5.416689e-13 erg/cm2/s (TS = 772.289)
2021-04-28T07:04:07:  Bin 2 .....................: 1.826697e-12 +/- 3.019027e-13 erg/cm2/s (TS = 43.077)
2021-04-28T07:04:07:  Bin 3 .....................: 3.697824e-13 +/- 3.024720e-13 erg/cm2/s (TS = 4.162)
2021-04-28T07:04:07:  Bin 4 .....................: 3.523538e-13 +/- 2.630165e-13 erg/cm2/s (TS = 15.176)
2021-04-28T07:04:07:  Bin 5 .....................: 2.945674e-13 +/- 0.000000e+00 erg/cm2/s (TS = 27.663)
2021-04-28T07:04:07:  Bin 6 .....................: 4.478029e-13 +/- 2.061951e-13 erg/cm2/s (TS = 129.545)
2021-04-28T07:04:07:  Bin 7 .....................: 4.654380e-13 +/- 6.369592e-13 erg/cm2/s (TS = 224.080)
2021-04-28T07:04:07:  Bin 8 .....................: 5.697476e-13 +/- 4.411305e-13 erg/cm2/s (TS = 501.889)
2021-04-28T07:04:07:  Bin 9 .....................: 6.441232e-13 +/- 7.019721e-13 erg/cm2/s (TS = 826.930)
2021-04-28T07:04:07:  Bin 10 ....................: 7.380578e-13 +/- 4.754821e-13 erg/cm2/s (TS = 1312.308)
2021-04-28T07:04:07:  Bin 11 ....................: 7.971278e-13 +/- 3.610715e-13 erg/cm2/s (TS = 1785.076)
2021-04-28T07:04:07:  Bin 12 ....................: 6.860986e-13 +/- 0.000000e+00 erg/cm2/s (TS = 1199.027)
2021-04-28T07:04:07:  Bin 13 ....................: 5.727868e-13 +/- 1.897827e-12 erg/cm2/s (TS = 845.929)
2021-04-28T07:04:07:  Bin 14 ....................: 4.012355e-13 +/- 4.052026e-13 erg/cm2/s (TS = 469.790)
2021-04-28T07:04:07:  Bin 15 ....................: 1.962754e-13 +/- 2.047394e-13 erg/cm2/s (TS = 175.953)
2021-04-28T07:04:07:  Bin 16 ....................: 7.272739e-14 +/- 6.218244e-14 erg/cm2/s
2021-04-28T07:04:07:  Bin 17 ....................: 1.331823e-14 +/- 2.182493e-14 erg/cm2/s
2021-04-28T07:04:07:  Bin 18 ....................: 1.133950e-16 +/- 0.000000e+00 erg/cm2/s
2021-04-28T07:04:07:  Bin 19 ....................: 3.557983e-19 +/- 0.000000e+00 erg/cm2/s
2021-04-28T07:04:07:  Bin 20 ....................: 4.073476e-23 +/- 0.000000e+00 erg/cm2/s

#10 Updated by Knödlseder Jürgen over 3 years ago

Below the corresponding spectra, left panel is for Aleksandr, right from my analysis. Note that there are fewer spectral points in the left panel, which is due to differences in the error estimates between the two analyses (the script only plots spectral points if the statistical flux error is smaller than the flux; in the other cases upper limits are plotted, yet since the upper limits are all set to zero, nothing is shown).

#11 Updated by Knödlseder Jürgen over 3 years ago

Here the results after adding the kludge, replacing the spectral model by a file function with 20 nodes within each energy bin. Things looks now much better, and also the spectrum that comes out looks fine.

2021-04-28T08:31:44: +===================+
2021-04-28T08:31:44: | Generate spectrum |
2021-04-28T08:31:44: +===================+
2021-04-28T08:31:44: === GEbounds ===
2021-04-28T08:31:44:  Number of intervals .......: 20
2021-04-28T08:31:44:  Energy range ..............: 30 GeV - 199 TeV
2021-04-28T08:54:49:  Bin 1 .....................: 1.428555e-11 +/- 5.380575e-13 [< 2.994514e-11] erg/cm2/s (TS = 772.341)
2021-04-28T08:54:49:  Bin 2 .....................: 1.809213e-12 +/- 2.801406e-13 [< 3.765293e-12] erg/cm2/s (TS = 43.195)
2021-04-28T08:54:49:  Bin 3 .....................: 3.312245e-13 +/- 1.569851e-13 [< 7.906755e-13] erg/cm2/s (TS = 4.528)
2021-04-28T08:54:49:  Bin 4 .....................: 3.522863e-13 +/- 9.332650e-14 [< 7.708354e-13] erg/cm2/s (TS = 15.080)
2021-04-28T08:54:49:  Bin 5 .....................: 2.864654e-13 +/- 5.824141e-14 [< 6.340299e-13] erg/cm2/s (TS = 27.521)
2021-04-28T08:54:49:  Bin 6 .....................: 4.664197e-13 +/- 4.885011e-14 [< 1.122977e-12] erg/cm2/s (TS = 129.358)
2021-04-28T08:54:49:  Bin 7 .....................: 4.693020e-13 +/- 4.224375e-14 [< 1.281683e-12] erg/cm2/s (TS = 224.009)
2021-04-28T08:54:49:  Bin 8 .....................: 5.916918e-13 +/- 4.261850e-14 [< 1.978291e-12] erg/cm2/s (TS = 502.580)
2021-04-28T08:54:49:  Bin 9 .....................: 6.579367e-13 +/- 4.261173e-14 [< 2.616455e-12] erg/cm2/s (TS = 827.177)
2021-04-28T08:54:49:  Bin 10 ....................: 7.431215e-13 +/- 4.259600e-14 [< 3.386862e-12] erg/cm2/s (TS = 1311.580)
2021-04-28T08:54:49:  Bin 11 ....................: 8.233050e-13 +/- 4.524378e-14 [< 4.342176e-12] erg/cm2/s (TS = 1785.263)
2021-04-28T08:54:49:  Bin 12 ....................: 6.813613e-13 +/- 4.668986e-14 [< 3.703221e-12] erg/cm2/s (TS = 1199.037)
2021-04-28T08:54:49:  Bin 13 ....................: 5.956886e-13 +/- 5.104040e-14 [< 3.454869e-12] erg/cm2/s (TS = 846.316)
2021-04-28T08:54:49:  Bin 14 ....................: 4.082578e-13 +/- 4.947661e-14 [< 2.550423e-12] erg/cm2/s (TS = 469.634)
2021-04-28T08:54:49:  Bin 15 ....................: 2.080721e-13 +/- 4.108257e-14 [< 1.309689e-12] erg/cm2/s (TS = 176.456)
2021-04-28T08:54:49:  Bin 16 ....................: 7.422105e-14 +/- 2.810298e-14 erg/cm2/s
2021-04-28T08:54:49:  Bin 17 ....................: 1.304906e-14 +/- 1.304906e-14 erg/cm2/s
2021-04-28T08:54:49:  Bin 18 ....................: 1.133622e-16 +/- 0.000000e+00 erg/cm2/s
2021-04-28T08:54:49:  Bin 19 ....................: 3.556966e-19 +/- 0.000000e+00 erg/cm2/s
2021-04-28T08:54:49:  Bin 20 ....................: 4.072843e-23 +/- 0.000000e+00 erg/cm2/s

#12 Updated by Knödlseder Jürgen over 3 years ago

  • Status changed from In Progress to Closed
  • % Done changed from 50 to 100

The code is not in the new release 1.7.4 that you can download as usual from our web site (http://cta.irap.omp.eu/ctools/index.html). The code is also merged in the devel branch. I close the issue now.

Also available in: Atom PDF