Adapt to new spectral model names (#1821)

Jürgen Knödlseder
1 parent d51de83b
Showing 11 changed files with 178 additions and 70 deletions
ChangeLog
NEWS
README.md
cscripts/cspull.py
doc/source/users/user_manual/getting_started/models.rst
test/data/obs_stacked_two.xml
test/data/sciver/crab_eplaw.xml
test/data/sciver/crab_plaw2.xml
test/data/sciver/crab_supeplaw.xml
test/science_verification.py
test/test_cspull.py
-2016-07-21  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
+2016-07-24  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
         * Version 1.1.0 released
         ========================
+        Adapt to new spatial (#1818) and spectral model names (#1821)
         Fix stacked analysis bias (#1819)
         Add background model science verification (#1815)
         Add cstsmapsplit script (#1797)
 New Features and Important Changes in ctools 1.1.0
-21 July 2016
+24 July 2016
 Introduction
 ctools information
 ==================
-* Version:             1.1.0.dev1 (21 July 2016)
+* Version:             1.1.0.dev1 (24 July 2016)
 * Author:              Juergen Knoedlseder (jurgen.knodlseder@irap.omp.eu)
 * GammaLib dependency: 1.1.0.dev1
@@ -173,8 +173,8 @@ class cspull(ctools.cscript):
         # Compute spatial binning for point spread function and
         # energy dispersion cubes
         binsz = 10.0 * self['binsz'].real()
-        nxpix = self['nxpix'].integer() // 10  # Make sure result is int
-        nypix = self['nypix'].integer() // 10  # Make sure result is int
+        nxpix = self['npix'].integer() // 10  # Make sure result is int
+        nypix = self['npix'].integer() // 10  # Make sure result is int
         if nxpix < 2:
             nxpix = 2
         if nypix < 2:
@@ -118,16 +118,14 @@ term will be appended to each ``&lt;parameter&gt;`` tag.
    parameters so that the effective parameters to be optimised (the ``value`` terms) 
    are all of about unity.
-.. warning::
+.. note::
-   The syntax of the model definition XML file has been adopted from the 
-   syntax used by the Fermi/LAT ScienceTools.
-   This should allow for exchange of model definition files between ctools 
-   and the Fermi/LAT ScienceTools.
-   Note, however, that ctools implements a number of extensions with 
-   respect to the Fermi/LAT ScienceTools syntax, hence eventually minor 
-   manual adjustments may be needed to use ctools XML files for 
-   Fermi/LAT analyses.
+   The syntax of the model definition XML file has been inspired from the
+   syntax used by the Fermi/LAT ScienceTools, but for reasons of clarity and
+   homogenity of the various model and parameter names we have made some
+   modifications.
+   Nevertheless, the exact format used by the Fermi/LAT ScienceTools is also
+   supported.
 .. _sec_spatial_models:
@@ -155,8 +153,8 @@ Point source
   .. note::
-    For compatibility with the Fermi/LAT ScienceTools the model type ``PointSource``
-    can be replaced by ``SkyDirFunction``.
+    For compatibility with the Fermi/LAT ScienceTools the model type
+    ``PointSource`` can be replaced by ``SkyDirFunction``.
 Radial source
@@ -373,8 +371,6 @@ in ctools.
 .. warning::
-   The units of the spectral model depend on the spatial model component.
-   Units quoted below apply to spatial models for a source component.
    Source intensities are generally given in units of
    :math:`{\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm MeV}^{-1}`.
@@ -405,11 +401,11 @@ Constant
   .. code-block:: xml
-   <spectrum type="ConstantValue">
+   <spectrum type="Constant">
      <parameter name="Normalization" scale="1e-16" value="5.7" min="1e-07" max="1000.0" free="1"/>
    </spectrum>
-  The ``ConstantValue`` model component implements the constant function
+  This spectral model component implements the constant function
   .. math::
     \frac{dN}{dE} = N_0
@@ -419,6 +415,11 @@ Constant
   * :math:`N_0` = ``Normalization``
     :math:`({\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm MeV}^{-1})`
+  .. note::
+
+    For compatibility with the Fermi/LAT ScienceTools the model type
+    ``Constant`` can be replaced by ``ConstantValue`` and the parameter
+    ``Normalization`` by ``Value``.
 Power law
@@ -432,7 +433,7 @@ Power law
      <parameter name="PivotEnergy" scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="0"/>
    </spectrum>
-  The ``PowerLaw`` model component implements the power law function
+  This spectral model component implements the power law function
   .. math::
     \frac{dN}{dE} = k_0 \left( \frac{E}{E_0} \right)^{\gamma}
@@ -454,19 +455,19 @@ Power law
     For compatibility with the Fermi/LAT ScienceTools the parameter
     ``PivotEnergy`` can be replaced by ``Scale``.
-  An alternative power law function that uses the integral photon flux as parameter
-  rather than the Prefactor is specified by
+  An alternative power law function that uses the integral photon flux as
+  parameter rather than the Prefactor is specified by
   .. code-block:: xml
-   <spectrum type="PowerLawPhotonFlux">
+   <spectrum type="PowerLaw">
      <parameter scale="1e-07" name="PhotonFlux" min="1e-07" max="1000.0"    value="1.0" free="1"/>
      <parameter scale="1.0"   name="Index"      min="-5.0"  max="+5.0"      value="-2.0" free="1"/>
      <parameter scale="1.0"   name="LowerLimit" min="10.0"  max="1000000.0" value="100.0" free="0"/>
      <parameter scale="1.0"   name="UpperLimit" min="10.0"  max="1000000.0" value="500000.0" free="0"/>
    </spectrum>
-  The ``PowerLawPhotonFlux`` model component implements the power law function
+  This spectral model component implements the power law function
   .. math::
     \frac{dN}{dE} = \frac{N(\gamma+1)E^{\gamma}}
@@ -484,8 +485,8 @@ Power law
   .. warning::
-    The ``LowerLimit`` and ``UpperLimit`` parameters are always treated as fixed and,
-    the flux given by the ``PhotonFlux`` parameter is computed over the
+    The ``LowerLimit`` and ``UpperLimit`` parameters are always treated as fixed
+    and the flux given by the ``PhotonFlux`` parameter is computed over the
     range set by these two parameters.
     Use of this model allows the errors on the integral flux to be evaluated directly
     by :ref:`ctlike`.
@@ -493,8 +494,8 @@ Power law
   .. note::
     For compatibility with the Fermi/LAT ScienceTools the model type
-    ``PowerLawPhotonFlux`` can be replaced by ``PowerLaw2`` and the parameter
-    ``PhotonFlux`` can be replaced by ``Integral``.
+    ``PowerLaw`` can be replaced by ``PowerLaw2`` and the parameter
+    ``PhotonFlux`` by ``Integral``.
 Exponentially cut-off power law
@@ -502,18 +503,19 @@ Exponentially cut-off power law
   .. code-block:: xml
-   <spectrum type="ExpCutoff">
-     <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="Cutoff"      scale="1e6"   value="1.0"  min="0.01"  max="1000.0" free="1"/>
-     <parameter name="PivotEnergy" scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="0"/>
+   <spectrum type="ExponentialCutoffPowerLaw">
+     <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="CutoffEnergy" scale="1e6"   value="1.0"  min="0.01"  max="1000.0" free="1"/>
+     <parameter name="PivotEnergy"  scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="0"/>
    </spectrum>
-  The ``ExpCutoff`` model component implements the exponentially cut-off power law function
+  This spectral model component implements the exponentially cut-off power law
+  function
   .. math::
     \frac{dN}{dE} = k_0 \left( \frac{E}{E_0} \right)^{\gamma}
-                    \exp \left( \frac{-E}{\tt E_{\rm cut}} \right)
+                    \exp \left( \frac{-E}{E_{\rm cut}} \right)
   where
@@ -522,7 +524,7 @@ Exponentially cut-off power law
   * :math:`\gamma` = ``Index``
   * :math:`E_0` = ``PivotEnergy``
     :math:`({\rm MeV})`
-  * :math:`E_{\rm cut}` = ``Cutoff``
+  * :math:`E_{\rm cut}` = ``CutoffEnergy``
     :math:`({\rm MeV})`
   .. warning::
@@ -531,8 +533,61 @@ Exponentially cut-off power law
   .. note::
-    For compatibility with the Fermi/LAT ScienceTools the parameter
-    ``PivotEnergy`` can be replaced by ``Scale``.
+    For compatibility with the Fermi/LAT ScienceTools the model type
+    ``ExponentialCutoffPowerLaw`` can be replaced by ``ExpCutoff`` and
+    the parameters ``CutoffEnergy`` by ``Cutoff`` and ``PivotEnergy``
+    by ``Scale``.
+
+
+Super exponentially cut-off power law
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+  .. code-block:: xml
+
+   <spectrum type="SuperExponentialCutoffPowerLaw">
+    <parameter name="Prefactor"    scale="1e-16" value="1.0" min="1e-07" max="1000.0" free="1"/>
+    <parameter name="Index1"       scale="-1"    value="2.0" min="0.0"   max="+5.0"   free="1"/>
+    <parameter name="CutoffEnergy" scale="1e6"   value="1.0" min="0.01"  max="1000.0" free="1"/>
+    <parameter name="Index2"       scale="1.0"   value="1.5" min="0.1"   max="5.0"    free="1"/>
+    <parameter name="PivotEnergy"  scale="1e6"   value="1.0" min="0.01"  max="1000.0" free="0"/>
+   </spectrum>
+
+  This spectral model component implements the super exponentially cut-off power
+  law function
+
+  .. math::
+    \frac{dN}{dE} = k_0 \left( \frac{E}{E_0} \right)^{\gamma}
+                    \exp \left( 
+                      -\left( \frac{E}{E_{\rm cut}} \right)^{\alpha}
+                    \right)
+
+  where
+
+  * :math:`k_0` = ``Prefactor``
+    :math:`({\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm MeV}^{-1})`
+  * :math:`\gamma` = ``Index1``
+  * :math:`\alpha` = ``Index2``
+  * :math:`E_0` = ``PivotEnergy``
+    :math:`({\rm MeV})`
+  * :math:`E_{\rm cut}` = ``CutoffEnergy``
+    :math:`({\rm MeV})`
+
+  .. warning::
+
+    The ``PivotEnergy`` parameter is not intended to be fitted.
+
+  An alternative XML format is supported for compatibility with the Fermi/LAT
+  XML format:
+
+  .. code-block:: xml
+
+   <spectrum type="PLSuperExpCutoff">
+    <parameter name="Prefactor"   scale="1e-16" value="1.0" min="1e-07" max="1000.0" free="1"/>
+    <parameter name="Index1"      scale="-1"    value="2.0" min="0.0"   max="+5.0"   free="1"/>
+    <parameter name="Cutoff"      scale="1e6"   value="1.0" min="0.01"  max="1000.0" free="1"/>
+    <parameter name="Index2"      scale="1.0"   value="1.5" min="0.1"   max="5.0"    free="1"/>
+    <parameter name="Scale"       scale="1e6"   value="1.0" min="0.01"  max="1000.0" free="0"/>
+   </spectrum>
 Broken power law
@@ -541,13 +596,13 @@ Broken power law
   .. code-block:: xml
    <spectrum type="BrokenPowerLaw">
-     <parameter name="Prefactor"  scale="1e-16" value="5.7"  min="1e-07" max="1000.0" free="1"/>
-     <parameter name="Index1"     scale="-1"    value="2.48" min="0.0"   max="+5.0"   free="1"/>
-     <parameter name="BreakValue" scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="1"/>
-     <parameter name="Index2"     scale="-1"    value="2.70" min="0.01"  max="1000.0" free="1"/>
+     <parameter name="Prefactor"   scale="1e-16" value="5.7"  min="1e-07" max="1000.0" free="1"/>
+     <parameter name="Index1"      scale="-1"    value="2.48" min="0.0"   max="+5.0"   free="1"/>
+     <parameter name="BreakEnergy" scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="1"/>
+     <parameter name="Index2"      scale="-1"    value="2.70" min="0.01"  max="1000.0" free="1"/>
    </spectrum>
-  The ``BrokenPowerLaw`` model component implements the broken power law function
+  This spectral model component implements the broken power law function
   .. math::
@@ -564,15 +619,20 @@ Broken power law
     :math:`({\rm ph}\,\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm MeV}^{-1})`
   * :math:`\gamma_1` = ``Index1``
   * :math:`\gamma_2` = ``Index2``
-  * :math:`E_b` = ``BreakValue``
+  * :math:`E_b` = ``BreakEnergy``
     :math:`({\rm MeV})`
   .. warning::
-    Note that the ``BreakValue`` parameter may be poorly constrained if 
+    Note that the ``BreakEnergy`` parameter may be poorly constrained if
     there is no clear spectral cut-off in the spectrum.
     This model may lead to complications in the maximum likelihood fitting.
+  .. note::
+
+    For compatibility with the Fermi/LAT ScienceTools the parameters
+    ``BreakEnergy`` can be replaced by ``BreakValue``.
+
 Log parabola
 ^^^^^^^^^^^^
@@ -586,7 +646,7 @@ Log parabola
      <parameter name="PivotEnergy" scale="1e6"   value="1.0"     min="0.01"  max="1000.0" free="0"/>
    </spectrum>
-  The ``LogParabola`` model component implements the log parabola function
+  This spectral model component implements the log parabola function
   .. math::
     \frac{dN}{dE} = k_0 \left( \frac{E}{E_0} \right)^{\gamma+\eta \ln(E/E_0)}
@@ -633,7 +693,7 @@ Gaussian
      <parameter name="Sigma"         scale="1e6"   value="1.0"  min="0.01"  max="100.0"  free="1"/>
    </spectrum>
-  The ``Gaussian`` model component implements the gaussian function
+  This spectral model component implements the gaussian function
   .. math::
     \frac{dN}{dE} = \frac{N_0}{\sqrt{2\pi}\sigma}
@@ -658,7 +718,7 @@ File function
      <parameter scale="1.0" name="Normalization" min="0.0" max="1000.0" value="1.0" free="1"/>
    </spectrum>
-  The ``FileFunction`` model component implements an arbitrary function 
+  This spectral model component implements an arbitrary function
   that is defined by intensity values at specific energies.
   The energy and intensity values are defined using an ASCII file with
   columns of energy and differential flux values.
@@ -677,9 +737,10 @@ File function
   .. warning::
-    If the file name is given as relative path, the path is relative to the
-    working directory where the ctool has been launched.
-    Alternatively, an absolute path may be specified.
+    If the file name is given without a path it is expected that the file
+    resides in the same directory than the XML file.
+    If the file resides in a different directory, an absolute path name should
+    be specified.
     Any environment variable present in the path name will be expanded.
@@ -699,7 +760,7 @@ Node function
      </node>
    </spectrum>
-  The ``NodeFunction`` model component implements a generalised broken 
+  This spectral model component implements a generalised broken 
   power law which is defined by a set of energy and intensity values
   (the so called nodes) that are piecewise connected by power laws.
   Energies are given in units of
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<observation_list title="observation list">
+  <observation name="Crab" id="000001" instrument="CTA">
+    <parameter name="Pointing" ra="83.6331" dec="21.0145" />
+    <parameter name="EnergyBoundaries" emin="50000" emax="1e+08" />
+    <parameter name="GoodTimeIntervals" tmin="0" tmax="1800" />
+    <parameter name="TimeReference" mjdrefi="51544" mjdreff="0.5" timeunit="s" timesys="TT" timeref="LOCAL" />
+    <parameter name="RegionOfInterest" ra="83.6331" dec="21.0145" rad="1.5" />
+    <parameter name="Deadtime" deadc="0.95" />
+    <parameter name="Calibration" database="prod2" response="South_0.5h" />
+  </observation>
+  <observation name="Crab" id="000002" instrument="CTA">
+    <parameter name="Pointing" ra="83.6331" dec="23.0145" />
+    <parameter name="EnergyBoundaries" emin="70000" emax="1e+08" />
+    <parameter name="GoodTimeIntervals" tmin="1800" tmax="3600" />
+    <parameter name="TimeReference" mjdrefi="51544" mjdreff="0.5" timeunit="s" timesys="TT" timeref="LOCAL" />
+    <parameter name="RegionOfInterest" ra="83.6331" dec="23.0145" rad="1.5" />
+    <parameter name="Deadtime" deadc="0.95" />
+    <parameter name="Calibration" database="prod2" response="South_0.5h" />
+  </observation>
+</observation_list>
 <?xml version="1.0" standalone="no"?>
 <source_library title="source library">
   <source name="Crab" type="PointSource">
-    <spectrum type="ExpCutoff">
-      <parameter name="Prefactor"   scale="1e-15" value="2.0"  min="1e-07" max="1000.0" free="1"/>
-      <parameter name="Index"       scale="-1"    value="2.1"  min="0.0"   max="+5.0"   free="1"/>
-      <parameter name="Cutoff"      scale="1e6"   value="1.0"  min="0.01"  max="1000.0" free="1"/>
-      <parameter name="PivotEnergy" scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="0"/>
+    <spectrum type="ExponentialCutoffPowerLaw">
+      <parameter name="Prefactor"    scale="1e-15" value="2.0"  min="1e-07" max="1000.0" free="1"/>
+      <parameter name="Index"        scale="-1"    value="2.1"  min="0.0"   max="+5.0"   free="1"/>
+      <parameter name="CutoffEnergy" scale="1e6"   value="1.0"  min="0.01"  max="1000.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"/>
 <?xml version="1.0" standalone="no"?>
 <source_library title="source library">
   <source name="Crab" type="PointSource">
-    <spectrum type="PowerLawPhotonFlux">
+    <spectrum type="PowerLaw">
       <parameter name="PhotonFlux" scale="1e-5" value="1.617"    min="1e-07" max="1000.0"    free="1"/>
       <parameter name="Index"      scale="-1.0" value="2.48"     min="-5.0"  max="+5.0"      free="1"/>
       <parameter name="LowerLimit" scale="1.0"  value="100.0"    min="10.0"  max="1000000.0" free="0"/>
 <?xml version="1.0" standalone="no"?>
 <source_library title="source library">
   <source name="Crab" type="PointSource">
-    <spectrum type="PLSuperExpCutoff">
-      <parameter name="Prefactor"   scale="1e-16" value="5.7"  min="1e-07" max="1000.0" free="1"/>
-      <parameter name="Index1"      scale="-1"    value="2.48" min="0.0"   max="+5.0"   free="1"/>
-      <parameter name="Cutoff"      scale="1e6"   value="1.0"  min="0.01"  max="1000.0" free="1"/>
-      <parameter name="PivotEnergy" scale="1e6"   value="1.0"  min="0.01"  max="1000.0" free="0"/>
-      <parameter name="Index2"      scale="1.0"   value="1.5"  min="0.1"   max="5.0"    free="1"/>
+    <spectrum type="SuperExponentialCutoffPowerLaw">
+      <parameter name="Prefactor"    scale="1e-16" value="5.7"  min="1e-07" max="1000.0" free="1"/>
+      <parameter name="Index1"       scale="-1"    value="2.48" min="0.0"   max="+5.0"   free="1"/>
+      <parameter name="CutoffEnergy" scale="1e6"   value="1.0"  min="0.01"  max="1000.0" free="1"/>
+      <parameter name="PivotEnergy"  scale="1e6"   value="1.0"  min="0.01"  max="1000.0" free="0"/>
+      <parameter name="Index2"       scale="1.0"   value="1.5"  min="0.1"   max="5.0"    free="1"/>
     </spectrum>
     <spatialModel type="PointSource">
       <parameter name="RA"  scale="1.0" value="83.6331" min="-360" max="360" free="0"/>
@@ -424,7 +424,7 @@ class sciver(gammalib.GPythonTestSuite):
         self.pull('data/sciver/crab_eplaw')
         self.test('Pull_Crab_Prefactor')
         self.test('Pull_Crab_Index')
-        self.test('Pull_Crab_Cutoff')
+        self.test('Pull_Crab_CutoffEnergy')
         self.test('Pull_Background_Prefactor')
         self.test('Pull_Background_Index')
         return
@@ -438,7 +438,7 @@ class sciver(gammalib.GPythonTestSuite):
         self.test('Pull_Crab_Prefactor')
         self.test('Pull_Crab_Index1')
         self.test('Pull_Crab_Index2')
-        self.test('Pull_Crab_Cutoff')
+        self.test('Pull_Crab_CutoffEnergy')
         self.test('Pull_Background_Prefactor')
         self.test('Pull_Background_Index')
         return
@@ -44,6 +44,7 @@ class Test(test):
         # Set members
         self._inobs         = self._datadir + '/obs_stacked.xml'
+        self._inobs_two     = self._datadir + '/obs_stacked_two.xml'
         self._stacked_model = self._datadir + '/crab_bkgcube.xml'
         # Return
@@ -209,14 +210,13 @@ class Test(test):
         # Check pull distribution file
         self._check_pull_file('cspull_py4.dat')
-        # Set-up stacked cspull
+        # Set-up stacked cspull with one observation where response cubes
+        # are specified in the input observation
         pull = cscripts.cspull()
         pull['inobs']    = self._inobs
         pull['inmodel']  = self._stacked_model
         pull['outfile']  = 'cspull_py5.dat'
         pull['ntrials']  = 3
-        pull['caldb']    = self._caldb
-        pull['irf']      = self._irf
         pull['enumbins'] = 0
         pull['logfile']  = 'cspull_py5.log'
         pull['chatter']  = 4
@@ -227,6 +227,31 @@ class Test(test):
         # Check pull distribution file
         self._check_pull_file('cspull_py5.dat')
+        # Set-up stacked cspull with two observations from which response
+        # cubes will be computed internally. The IRFs are also specified
+        # in the input observation, hence we do not need to query the IRF
+        # parameters.
+        pull = cscripts.cspull()
+        pull['inobs']    = self._inobs_two
+        pull['inmodel']  = self._model
+        pull['outfile']  = 'cspull_py6.dat'
+        pull['ntrials']  = 3
+        pull['emin']     = 0.02
+        pull['emax']     = 100.0
+        pull['enumbins'] = 10
+        pull['npix']     = 20
+        pull['binsz']    = 0.2
+        pull['coordsys'] = 'CEL'
+        pull['proj']     = 'TAN'
+        pull['logfile']  = 'cspull_py6.log'
+        pull['chatter']  = 4
+
+        # Execute cspull script
+        pull.execute()
+
+        # Check pull distribution file
+        self._check_pull_file('cspull_py6.dat')
+
         # Return
         return