Action #2868: Add radial CTA Gaussian background model with energy dependence of sigma - GammaLib - CTA IRAP Project Gateway

Action #2868

Add radial CTA Gaussian background model with energy dependence of sigma

Added by Knödlseder Jürgen about 5 years ago. Updated almost 5 years ago.

Status:

Closed

Start date:

04/13/2019

Priority:

Urgent

Due date:

Assigned To:

Knödlseder Jürgen

% Done:

100%

Category:

Target version:

1.6.0

Duration:

Description

The analysis of the H.E.S.S. observations of RX J1713 makes it necessary to have a radial Gaussian model where sigma depends on energy. This can be implemented by having a class with a spectral model member of type GModelSpectral, where the spectral law is used to compute sigma.

I suggest a class GCTAModelSpatialGaussSpectrum that implements a spatial model of type EnergyDependentGaussian which has a <sigma> tab in the XML file:

  <source name="My model" type="CTABackground" instrument="CTA">
    <spatialModel type="EnergyDependentGaussian">
      <sigma type="PowerLaw">
        <parameter name="Prefactor"   scale="1"   value="1" min="0.01" max="1000.0" free="1"/>
        <parameter name="Index"       scale="1"   value="0" min="-5.0" max="+5.0"   free="1"/>
        <parameter name="PivotEnergy" scale="1e6" value="1" min="0.01" max="1000.0" free="0"/>
      </sigma>
    </spatialModel>
    <spectrum type="...">
       ...
    </spectrum>
  </source>

Recurrence

No recurrence.

History

#1 Updated by Knödlseder Jürgen about 5 years ago

Status changed from New to In Progress
% Done changed from 0 to 50

I implemented a first version of GCTAModelSpatialGaussSpectrum and applied it to an Off observation using a constant spectral model for sigma. Here the beginning of the XML file:

<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<source_library title="source library">
  <source name="Background_020275" type="CTABackground" instrument="HESS" id="020275">
    <spatialModel type="Multiplicative">
      <spatialModel type="EnergyDependentGaussian">
        <sigma type="Constant">
          <parameter name="Normalization" scale="1" value="3.0"  min="1.5" max="10.0" free="1"/>
        </sigma>
      </spatialModel>
      <spatialModel type="Gradient">
        <parameter name="Grad_DETX" value="0.0449630797077924" error="0.0144756636305313" scale="1" free="1" />
        <parameter name="Grad_DETY" value="0.0555447282348084" error="0.0147448614132178" scale="1" free="1" />
      </spatialModel>
    </spatialModel>
    <spectrum type="NodeFunction">
      ...
    </spectrum>
  </source>
</source_library>

And here the resulting model fit:

2019-04-13T13:30:23: === GOptimizerLM ===
2019-04-13T13:30:23:  Optimized function value ..: 41650.817
2019-04-13T13:30:23:  Absolute precision ........: 0.005
2019-04-13T13:30:23:  Acceptable value decrease .: 2
2019-04-13T13:30:23:  Optimization status .......: converged
2019-04-13T13:30:23:  Number of parameters ......: 20
2019-04-13T13:30:23:  Number of free parameters .: 11
2019-04-13T13:30:23:  Number of iterations ......: 8
2019-04-13T13:30:23:  Lambda ....................: 1e-05
2019-04-13T13:30:23:  Maximum log likelihood ....: -41650.817
2019-04-13T13:30:23:  Observed events  (Nobs) ...: 5222.000
2019-04-13T13:30:23:  Predicted events (Npred) ..: 5221.998 (Nobs - Npred = 0.0019373436089154)
2019-04-13T13:30:23: === GModels ===
2019-04-13T13:30:23:  Number of models ..........: 1
2019-04-13T13:30:23:  Number of parameters ......: 20
2019-04-13T13:30:23: === GCTAModelBackground ===
2019-04-13T13:30:23:  Name ......................: Background_020275
2019-04-13T13:30:23:  Instruments ...............: HESS
2019-04-13T13:30:23:  Instrument scale factors ..: unity
2019-04-13T13:30:23:  Observation identifiers ...: 020275
2019-04-13T13:30:23:  Model type ................: CTABackground
2019-04-13T13:30:23:  Model components ..........: "Multiplicative" * "NodeFunction" * "Constant" 
2019-04-13T13:30:23:  Number of parameters ......: 20
2019-04-13T13:30:23:  Number of spatial par's ...: 3
2019-04-13T13:30:23:   1:Normalization ..........: 3.57531367292878 +/- 0.140346064408971 [1.5,10]  (free,scale=1,gradient)
2019-04-13T13:30:23:   2:Grad_DETX ..............: 0.0449122753455452 +/- 0.014465368333852 deg^-1 (free,scale=1,gradient)
2019-04-13T13:30:23:   2:Grad_DETY ..............: 0.0556605852663954 +/- 0.014729652264425 deg^-1 (free,scale=1,gradient)

For comparison, here the model fit with the usual radial Gaussian model. The results are basically the same, though they are not identical.

2019-04-13T13:20:28: === GOptimizerLM ===
2019-04-13T13:20:28:  Optimized function value ..: 41650.817
2019-04-13T13:20:28:  Absolute precision ........: 0.005
2019-04-13T13:20:28:  Acceptable value decrease .: 2
2019-04-13T13:20:28:  Optimization status .......: converged
2019-04-13T13:20:28:  Number of parameters ......: 20
2019-04-13T13:20:28:  Number of free parameters .: 11
2019-04-13T13:20:28:  Number of iterations ......: 7
2019-04-13T13:20:28:  Lambda ....................: 0.0001
2019-04-13T13:20:28:  Maximum log likelihood ....: -41650.817
2019-04-13T13:20:28:  Observed events  (Nobs) ...: 5222.000
2019-04-13T13:20:28:  Predicted events (Npred) ..: 5221.993 (Nobs - Npred = 0.00708173094517406)
2019-04-13T13:20:28: === GModels ===
2019-04-13T13:20:28:  Number of models ..........: 1
2019-04-13T13:20:28:  Number of parameters ......: 20
2019-04-13T13:20:28: === GCTAModelBackground ===
2019-04-13T13:20:28:  Name ......................: Background_020275
2019-04-13T13:20:28:  Instruments ...............: HESS
2019-04-13T13:20:28:  Instrument scale factors ..: unity
2019-04-13T13:20:28:  Observation identifiers ...: 020275
2019-04-13T13:20:28:  Model type ................: CTABackground
2019-04-13T13:20:28:  Model components ..........: "Multiplicative" * "NodeFunction" * "Constant" 
2019-04-13T13:20:28:  Number of parameters ......: 20
2019-04-13T13:20:28:  Number of spatial par's ...: 3
2019-04-13T13:20:28:   1:Sigma ..................: 3.57489131652591 +/- 0.140296325361186 [1.5,10] deg2 (free,scale=1,gradient)
2019-04-13T13:20:28:   2:Grad_DETX ..............: 0.0449117605696517 +/- 0.0144653779058488 deg^-1 (free,scale=1,gradient)
2019-04-13T13:20:28:   2:Grad_DETY ..............: 0.0556527129382587 +/- 0.0147296777463999 deg^-1 (free,scale=1,gradient)

#2 Updated by Knödlseder Jürgen about 5 years ago

% Done changed from 50 to 60

I now tried fitting with a two-node function for the sigma energy dependence. Here again the beginning of the XML file:

<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<source_library title="source library">
  <source name="Background_020275" type="CTABackground" instrument="HESS" id="020275">
    <spatialModel type="Multiplicative">
      <spatialModel type="EnergyDependentGaussian">
        <sigma type="NodeFunction">
          <node>
            <parameter name="Energy"    value="1"   scale="300000" free="0" />
            <parameter name="Intensity" value="3.0" scale="1.0"    min="1.5" max="100.0" free="1" />
          </node>
          <node>
            <parameter name="Energy"    value="1"   scale="10000000" free="0" />
            <parameter name="Intensity" value="4.0" scale="1.0"      min="1.5" max="100.0"  free="1" />
          </node>
        </sigma>
      </spatialModel>
      <spatialModel type="Gradient">
        <parameter name="Grad_DETX" value="0.0449630797077924" error="0.0144756636305313" scale="1" free="1" />
        <parameter name="Grad_DETY" value="0.0555447282348084" error="0.0147448614132178" scale="1" free="1" />
      </spatialModel>
    </spatialModel>
    <spectrum type="NodeFunction">

And here the result. The fit clearly found a sigma value that increases with energy.

2019-04-13T13:40:26: === GOptimizerLM ===
2019-04-13T13:40:26:  Optimized function value ..: 41637.746
2019-04-13T13:40:26:  Absolute precision ........: 0.005
2019-04-13T13:40:26:  Acceptable value decrease .: 2
2019-04-13T13:40:26:  Optimization status .......: converged
2019-04-13T13:40:26:  Number of parameters ......: 23
2019-04-13T13:40:26:  Number of free parameters .: 12
2019-04-13T13:40:26:  Number of iterations ......: 8
2019-04-13T13:40:26:  Lambda ....................: 1e-05
2019-04-13T13:40:26:  Maximum log likelihood ....: -41637.746
2019-04-13T13:40:26:  Observed events  (Nobs) ...: 5222.000
2019-04-13T13:40:26:  Predicted events (Npred) ..: 5221.993 (Nobs - Npred = 0.00719846033371141)
2019-04-13T13:40:26: === GModels ===
2019-04-13T13:40:26:  Number of models ..........: 1
2019-04-13T13:40:26:  Number of parameters ......: 23
2019-04-13T13:40:26: === GCTAModelBackground ===
2019-04-13T13:40:26:  Name ......................: Background_020275
2019-04-13T13:40:26:  Instruments ...............: HESS
2019-04-13T13:40:26:  Instrument scale factors ..: unity
2019-04-13T13:40:26:  Observation identifiers ...: 020275
2019-04-13T13:40:26:  Model type ................: CTABackground
2019-04-13T13:40:26:  Model components ..........: "Multiplicative" * "NodeFunction" * "Constant" 
2019-04-13T13:40:26:  Number of parameters ......: 23
2019-04-13T13:40:26:  Number of spatial par's ...: 6
2019-04-13T13:40:26:   Energy0 ..................: 300000 MeV (fixed,scale=300000)
2019-04-13T13:40:26:   Intensity0 ...............: 2.90914400029964 +/- 0.164805004035297 [1.5,100] ph/cm2/s/MeV (free,scale=1,gradient)
2019-04-13T13:40:26:   Energy1 ..................: 10000000 MeV (fixed,scale=10000000)
2019-04-13T13:40:26:   Intensity1 ...............: 9.76611185305207 +/- 3.50792380821302 [1.5,100] ph/cm2/s/MeV (free,scale=1,gradient)
2019-04-13T13:40:26:   2:Grad_DETX ..............: 0.0449196320577307 +/- 0.0144696377910105 deg^-1 (free,scale=1,gradient)
2019-04-13T13:40:26:   2:Grad_DETY ..............: 0.055651617302785 +/- 0.0147297692879211 deg^-1 (free,scale=1,gradient)

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

And here the results for 3 nodes. Things seem to work well.

2019-04-13T13:45:38: === GOptimizerLM ===
2019-04-13T13:45:38:  Optimized function value ..: 41636.950
2019-04-13T13:45:38:  Absolute precision ........: 0.005
2019-04-13T13:45:38:  Acceptable value decrease .: 2
2019-04-13T13:45:38:  Optimization status .......: converged
2019-04-13T13:45:38:  Number of parameters ......: 25
2019-04-13T13:45:38:  Number of free parameters .: 13
2019-04-13T13:45:38:  Number of iterations ......: 9
2019-04-13T13:45:38:  Lambda ....................: 1e-06
2019-04-13T13:45:38:  Maximum log likelihood ....: -41636.950
2019-04-13T13:45:38:  Observed events  (Nobs) ...: 5222.000
2019-04-13T13:45:38:  Predicted events (Npred) ..: 5221.923 (Nobs - Npred = 0.0772168028734086)
2019-04-13T13:45:38: === GModels ===
2019-04-13T13:45:38:  Number of models ..........: 1
2019-04-13T13:45:38:  Number of parameters ......: 25
2019-04-13T13:45:38: === GCTAModelBackground ===
2019-04-13T13:45:38:  Name ......................: Background_020275
2019-04-13T13:45:38:  Instruments ...............: HESS
2019-04-13T13:45:38:  Instrument scale factors ..: unity
2019-04-13T13:45:38:  Observation identifiers ...: 020275
2019-04-13T13:45:38:  Model type ................: CTABackground
2019-04-13T13:45:38:  Model components ..........: "Multiplicative" * "NodeFunction" * "Constant" 
2019-04-13T13:45:38:  Number of parameters ......: 25
2019-04-13T13:45:38:  Number of spatial par's ...: 8
2019-04-13T13:45:38:   Energy0 ..................: 300000 MeV (fixed,scale=300000)
2019-04-13T13:45:38:   Intensity0 ...............: 3.06016962093601 +/- 0.196471316341562 [1.5,100] ph/cm2/s/MeV (free,scale=1,gradient)
2019-04-13T13:45:38:   Energy1 ..................: 1000000 MeV (fixed,scale=1000000)
2019-04-13T13:45:38:   Intensity1 ...............: 3.87176208906406 +/- 0.453552927361097 [1.5,100] ph/cm2/s/MeV (free,scale=1,gradient)
2019-04-13T13:45:38:   Energy2 ..................: 10000000 MeV (fixed,scale=10000000)
2019-04-13T13:45:38:   Intensity2 ...............: 19.4095263235039 +/- 25.523403379984 [1.5,100] ph/cm2/s/MeV (free,scale=1,gradient)
2019-04-13T13:45:38:   2:Grad_DETX ..............: 0.0449309872630205 +/- 0.0144691948908213 deg^-1 (free,scale=1,gradient)
2019-04-13T13:45:38:   2:Grad_DETY ..............: 0.0556728207404466 +/- 0.0147337594535947 deg^-1 (free,scale=1,gradient)

#4 Updated by Knödlseder Jürgen about 5 years ago

% Done changed from 60 to 70

The new code was merged into devel.

Next step is to update csbkgmodel to take the new model into account. This is treated for simplicity in the same issue.

Here an excerpt of the proposed new parameter file:

#
# Script parameters
#==================
instrument, s, a, CTA,,, "Instrument name (for multi-instrument analysis)" 
spatial,    s, a, GAUSS,IRF|AEFF|GAUSS|GAUSS(E),, "Spatial model component" 
snumbins,   i, a, 3,1,200, "Number of energy nodes for sigma(E)" 
smin,       r, a, 0.1,,, "Lower energy limit for sigma(E) energy nodes (TeV)" 
smax,       r, a, 10.0,,, "Upper energy limit for sigma(E) energy nodes (TeV)" 
gradient,   b, a, yes,,, "Allow for spatial gradient?" 
spectral,   s, a, NODES,PLAW|NODES,, "Spectral model component" 
ebinalg,    s, a, LOG,FILE|LIN|LOG,, "Algorithm for defining spectral energy nodes" 
emin,       r, a, 0.1,,, "Lower energy limit (TeV)" 
emax,       r, a, 100.0,,, "Upper energy limit (TeV)" 
enumbins,   i, a, 8,1,200, "Number of spectral energy nodes" 
ebinfile,   f, a, NONE,,, "Name of the file containing the energy node definition" 
runwise,    b, a, yes,,, "Generate runwise background model?" 
rad,        r, h, 2.0,,, "Radius for event selection (degrees)"

A new spatial model GAUSS(E) was added, and new parameters snumbins, smin and smax were added that control the number and energy range of the the nodes for sigma. If GAUSS(E) is choses a EnergyDependentGaussian model will be used as spatial model with a NodeFunction as spectral law. This seems to work well for RX J1713. Before adopting the model should also be tested on other sources.

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

Status changed from In Progress to Closed
% Done changed from 70 to 100

Code was merged into devel.

Also available in: Atom PDF

GammaLib

Issues

Tags