Jürgen Knödlseder / gammalib

Commit cc742ac2c9badf7a736861f3813fc1acb87d5a9e

Authored by Jürgen Knödlseder
1 parent e2622d33

Rework interface of eval, eval_gradients and mc methods. 

The eval() and eval_gradients() methods now take a GPhoton object as argument to allow the implementation of energy and time dependent source morphologies.

GEnergy and GTime objects have been added to the mc() method interface to support equally the implementation of energy and time dependent source morphologies.

So far, no class really implements energy and time dependent source morphologies, but now we have the infrastructure to do so.
Inline Side-by-side
Showing 37 changed files with 1007 additions and 364 deletions
include/GModelSky.hpp
  View file @cc742ac
... ... @@ -59,26 +59,36 @@ class GObservation;
59 59 * a spectral and a temporal component. The factorization is given by
60 60 *
61 61 * \f[
62   - * S(\vec{p}, E, t) = S_{\rm p}(\vec{p}, E, t) \,
63   - * S_{\rm E}(\vec{p}, E, t) \,
64   - * S_{\rm t}(\vec{p}, E, t)
  62 + * S(\vec{p}, E, t) = S_{\rm p}(\vec{p} | E, t) \,
  63 + * S_{\rm E}(E | t) \,
  64 + * S_{\rm t}(t)
65 65 * \f]
66 66 *
67 67 * where
68   - * - \f$S_{\rm p}(\vec{p}, E, t)\f$ is the spatial,
69   - * - \f$S_{\rm E}(\vec{p}, E, t)\f$ is the spectral, and
70   - * - \f$S_{\rm t}(\vec{p}, E, t)\f$ is the temporal component of the
71   - * model.
72   - *
73   - * Note that the spatial, spectral and temporal components depend on the
74   - * three photon properties (true sky direction \f$\vec{p}\f$, true energy
75   - * \f$E\f$ and true arrival time \f$t\f$) so that any correlation between
76   - * these properties can be implemented. The factorization is in the sens
77   - * that the spatial component shall return a normalization that depends
78   - * on the spatial morphology distribution, the spectral component shall
79   - * return an intensity that depends on the spectral energy distribution,
80   - * and the temporal component shall return a normalization that depends
81   - * on the light curve of the source.
  68 + * - \f$S_{\rm p}(\vec{p} | E, t)\f$ is the spatial component,
  69 + * - \f$S_{\rm E}(E | t)\f$ is the spectral component, and
  70 + * - \f$S_{\rm t}(t)\f$ is the temporal component.
  71 + *
  72 + * The spatial component describes the energy and time dependent morphology
  73 + * of the source. It satisfies
  74 + * \f[
  75 + * \int_{\Omega} S_{\rm p}(\vec{p} | E, t) d\Omega = 1
  76 + * \f]
  77 + * for all \f$E\f$ and \f$t\f$, hence the spatial component does not
  78 + * impact the spatially integrated spectral and temporal properties of the
  79 + * source.
  80 + *
  81 + * The spectral component describes the spatially integrated time dependent
  82 + * spectral distribution of the source. It satisfies
  83 + * \f[
  84 + * \int_{E} S_{\rm E}(E | t) dE = \Phi
  85 + * \f]
  86 + * for all \f$t\f$, where \f$\Phi\f$ is the spatially and spectrally
  87 + * integrated total source flux. The spectral component does not impact
  88 + * the temporal properties of the integrated flux \f$\Phi\f$.
  89 + *
  90 + * The temporal component describes the temporal variation of the total
  91 + * integrated flux \f$\Phi\f$ of the source.
82 92 *
83 93 * The class has two methods for model evaluation that evaluate the model
84 94 * for a specific event, given an observation. The eval() method returns
... ...
include/GModelSpatial.hpp
  View file @cc742ac
... ... @@ -32,7 +32,10 @@
32 32 #include <vector>
33 33 #include "GBase.hpp"
34 34 #include "GModelPar.hpp"
  35 +#include "GPhoton.hpp"
35 36 #include "GSkyDir.hpp"
  37 +#include "GEnergy.hpp"
  38 +#include "GTime.hpp"
36 39 #include "GXmlElement.hpp"
37 40 #include "GRan.hpp"
38 41  
... ... @@ -64,9 +67,10 @@ public:
64 67 virtual void clear(void) = 0;
65 68 virtual GModelSpatial* clone(void) const = 0;
66 69 virtual std::string type(void) const = 0;
67   - virtual double eval(const GSkyDir& srcDir) const = 0;
68   - virtual double eval_gradients(const GSkyDir& srcDir) const = 0;
69   - virtual GSkyDir mc(GRan& ran) const = 0;
  70 + virtual double eval(const GPhoton& photon) const = 0;
  71 + virtual double eval_gradients(const GPhoton& photon) const = 0;
  72 + virtual GSkyDir mc(const GEnergy& energy, const GTime& time,
  73 + GRan& ran) const = 0;
70 74 virtual void read(const GXmlElement& xml) = 0;
71 75 virtual void write(GXmlElement& xml) const = 0;
72 76 virtual std::string print(void) const = 0;
... ... @@ -85,6 +89,7 @@ protected:
85 89 std::vector<GModelPar*> m_pars; //!< Parameter pointers
86 90 };
87 91  
  92 +
88 93 /***********************************************************************//**
89 94 * @brief Return number of parameters
90 95 *
... ...
include/GModelSpatialDiffuse.hpp
  View file @cc742ac
... ... @@ -30,7 +30,10 @@
30 30 /* __ Includes ___________________________________________________________ */
31 31 #include <string>
32 32 #include "GModelSpatial.hpp"
  33 +#include "GPhoton.hpp"
33 34 #include "GSkyDir.hpp"
  35 +#include "GEnergy.hpp"
  36 +#include "GTime.hpp"
34 37 #include "GXmlElement.hpp"
35 38 #include "GRan.hpp"
36 39  
... ... @@ -58,9 +61,10 @@ public:
58 61 virtual void clear(void) = 0;
59 62 virtual GModelSpatialDiffuse* clone(void) const = 0;
60 63 virtual std::string type(void) const = 0;
61   - virtual double eval(const GSkyDir& srcDir) const = 0;
62   - virtual double eval_gradients(const GSkyDir& srcDir) const = 0;
63   - virtual GSkyDir mc(GRan& ran) const = 0;
  64 + virtual double eval(const GPhoton& photon) const = 0;
  65 + virtual double eval_gradients(const GPhoton& photon) const = 0;
  66 + virtual GSkyDir mc(const GEnergy& energy, const GTime& time,
  67 + GRan& ran) const = 0;
64 68 virtual void read(const GXmlElement& xml) = 0;
65 69 virtual void write(GXmlElement& xml) const = 0;
66 70 virtual std::string print(void) const = 0;
... ...
include/GModelSpatialDiffuseConst.hpp
  View file @cc742ac
... ... @@ -60,9 +60,11 @@ public:
60 60 virtual void clear(void);
61 61 virtual GModelSpatialDiffuseConst* clone(void) const;
62 62 virtual std::string type(void) const;
63   - virtual double eval(const GSkyDir& srcDir) const;
64   - virtual double eval_gradients(const GSkyDir& srcDir) const;
65   - virtual GSkyDir mc(GRan& ran) const;
  63 + virtual double eval(const GPhoton& photon) const;
  64 + virtual double eval_gradients(const GPhoton& photon) const;
  65 + virtual GSkyDir mc(const GEnergy& energy,
  66 + const GTime& time,
  67 + GRan& ran) const;
66 68 virtual void read(const GXmlElement& xml);
67 69 virtual void write(GXmlElement& xml) const;
68 70 virtual std::string print(void) const;
... ...
include/GModelSpatialDiffuseCube.hpp
  View file @cc742ac
... ... @@ -69,9 +69,11 @@ public:
69 69 virtual void clear(void);
70 70 virtual GModelSpatialDiffuseCube* clone(void) const;
71 71 virtual std::string type(void) const;
72   - virtual double eval(const GSkyDir& srcDir) const;
73   - virtual double eval_gradients(const GSkyDir& srcDir) const;
74   - virtual GSkyDir mc(GRan& ran) const;
  72 + virtual double eval(const GPhoton& photon) const;
  73 + virtual double eval_gradients(const GPhoton& photon) const;
  74 + virtual GSkyDir mc(const GEnergy& energy,
  75 + const GTime& time,
  76 + GRan& ran) const;
75 77 virtual void read(const GXmlElement& xml);
76 78 virtual void write(GXmlElement& xml) const;
77 79 virtual std::string print(void) const;
... ...
include/GModelSpatialDiffuseMap.hpp
  View file @cc742ac
... ... @@ -61,10 +61,12 @@ public:
61 61 // Implemented pure virtual methods
62 62 virtual void clear(void);
63 63 virtual GModelSpatialDiffuseMap* clone(void) const;
64   - virtual std::string type(void) const { return "SpatialMap"; }
65   - virtual double eval(const GSkyDir& srcDir) const;
66   - virtual double eval_gradients(const GSkyDir& srcDir) const;
67   - virtual GSkyDir mc(GRan& ran) const;
  64 + virtual std::string type(void) const;
  65 + virtual double eval(const GPhoton& photon) const;
  66 + virtual double eval_gradients(const GPhoton& photon) const;
  67 + virtual GSkyDir mc(const GEnergy& energy,
  68 + const GTime& time,
  69 + GRan& ran) const;
68 70 virtual void read(const GXmlElement& xml);
69 71 virtual void write(GXmlElement& xml) const;
70 72 virtual std::string print(void) const;
... ... @@ -83,4 +85,18 @@ protected:
83 85 std::vector<double> m_mc_cache; //!< Monte Carlo cache
84 86 };
85 87  
  88 +/***********************************************************************//**
  89 + * @brief Return spatial model type
  90 + *
  91 + * @return "SpatialMap".
  92 + *
  93 + * Returns the type of the spatial map model.
  94 + ***************************************************************************/
  95 +inline
  96 +std::string GModelSpatialDiffuseMap::type(void) const
  97 +{
  98 + return "SpatialMap";
  99 +}
  100 +
  101 +
86 102 #endif /* GMODELSPATIALDIFFUSEMAP_HPP */
... ...
include/GModelSpatialElliptical.hpp
  View file @cc742ac
... ... @@ -33,6 +33,9 @@
33 33 #include "GModelSpatial.hpp"
34 34 #include "GModelPar.hpp"
35 35 #include "GSkyDir.hpp"
  36 +#include "GEnergy.hpp"
  37 +#include "GTime.hpp"
  38 +#include "GPhoton.hpp"
36 39 #include "GXmlElement.hpp"
37 40 #include "GRan.hpp"
38 41  
... ... @@ -62,27 +65,33 @@ public:
62 65 virtual void clear(void) = 0;
63 66 virtual GModelSpatialElliptical* clone(void) const = 0;
64 67 virtual std::string type(void) const = 0;
65   - virtual double eval(const double& theta,
66   - const double& posangle) const = 0;
67   - virtual double eval_gradients(const double& theta,
68   - const double& posangle) const = 0;
69   - virtual GSkyDir mc(GRan& ran) const = 0;
  68 + virtual double eval(const double& theta,
  69 + const double& posangle,
  70 + const GEnergy& energy,
  71 + const GTime& time) const = 0;
  72 + virtual double eval_gradients(const double& theta,
  73 + const double& posangle,
  74 + const GEnergy& energy,
  75 + const GTime& time) const = 0;
  76 + virtual GSkyDir mc(const GEnergy& energy,
  77 + const GTime& time,
  78 + GRan& ran) const = 0;
70 79 virtual double theta_max(void) const = 0;
71 80 virtual std::string print(void) const = 0;
72 81  
73   - // Implemented virtual methods
74   - virtual double eval(const GSkyDir& srcDir) const;
75   - virtual double eval_gradients(const GSkyDir& srcDir) const;
  82 + // Implemented virtual base class methods
  83 + virtual double eval(const GPhoton& photon) const;
  84 + virtual double eval_gradients(const GPhoton& photon) const;
76 85 virtual void read(const GXmlElement& xml);
77 86 virtual void write(GXmlElement& xml) const;
78 87  
79 88 // Other methods
80   - double ra(void) const { return m_ra.value(); }
81   - double dec(void) const { return m_dec.value(); }
82   - double posangle(void) const { return m_posangle.value(); }
  89 + double ra(void) const;
  90 + double dec(void) const;
  91 + double posangle(void) const;
83 92 GSkyDir dir(void) const;
84 93 void dir(const GSkyDir& dir);
85   - void posangle(const double& posangle) { m_posangle.value(posangle); }
  94 + void posangle(const double& posangle);
86 95  
87 96 protected:
88 97 // Protected methods
... ... @@ -96,4 +105,63 @@ protected:
96 105 GModelPar m_posangle; //!< Position angle from North, counterclockwise (deg)
97 106 };
98 107  
  108 +
  109 +/***********************************************************************//**
  110 + * @brief Return Right Ascencion of model centre
  111 + *
  112 + * @return Right Ascencion of model centre (degrees).
  113 + *
  114 + * Returns the Right Ascension of the model centre in degrees.
  115 + ***************************************************************************/
  116 +inline
  117 +double GModelSpatialElliptical::ra(void) const
  118 +{
  119 + return (m_ra.value());
  120 +}
  121 +
  122 +
  123 +/***********************************************************************//**
  124 + * @brief Return Declination of model centre
  125 + *
  126 + * @return Declination of model centre (degrees).
  127 + *
  128 + * Returns the Declination of the model centre in degrees.
  129 + ***************************************************************************/
  130 +inline
  131 +double GModelSpatialElliptical::dec(void) const
  132 +{
  133 + return (m_dec.value());
  134 +}
  135 +
  136 +
  137 +/***********************************************************************//**
  138 + * @brief Return Position Angle of model
  139 + *
  140 + * @return Position Angle of model (degrees).
  141 + *
  142 + * Returns the Position Angle of model in degrees, measured counterclockwise
  143 + * from celestial North.
  144 + ***************************************************************************/
  145 +inline
  146 +double GModelSpatialElliptical::posangle(void) const
  147 +{
  148 + return (m_posangle.value());
  149 +}
  150 +
  151 +
  152 +/***********************************************************************//**
  153 + * @brief Set Position Angle of model
  154 + *
  155 + * @param[in] posangle Position Angle of model (degrees).
  156 + *
  157 + * Sets the Position Angle of model in degrees, measured counterclockwise
  158 + * from celestial North.
  159 + ***************************************************************************/
  160 +inline
  161 +void GModelSpatialElliptical::posangle(const double& posangle)
  162 +{
  163 + m_posangle.value(posangle);
  164 + return;
  165 +}
  166 +
99 167 #endif /* GMODELSPATIALELLIPTICAL_HPP */
... ...
include/GModelSpatialEllipticalDisk.hpp
  View file @cc742ac
... ... @@ -63,22 +63,28 @@ public:
63 63 // Implemented pure virtual base class methods
64 64 virtual void clear(void);
65 65 virtual GModelSpatialEllipticalDisk* clone(void) const;
66   - virtual std::string type(void) const { return "EllipticalDisk"; }
67   - virtual double eval(const double& theta,
68   - const double& posangle) const;
69   - virtual double eval_gradients(const double& theta,
70   - const double& posangle) const;
71   - virtual GSkyDir mc(GRan& ran) const;
  66 + virtual std::string type(void) const;
  67 + virtual double eval(const double& theta,
  68 + const double& posangle,
  69 + const GEnergy& energy,
  70 + const GTime& time) const;
  71 + virtual double eval_gradients(const double& theta,
  72 + const double& posangle,
  73 + const GEnergy& energy,
  74 + const GTime& time) const;
  75 + virtual GSkyDir mc(const GEnergy& energy,
  76 + const GTime& time,
  77 + GRan& ran) const;
72 78 virtual double theta_max(void) const;
73 79 virtual void read(const GXmlElement& xml);
74 80 virtual void write(GXmlElement& xml) const;
75 81 virtual std::string print(void) const;
76 82  
77 83 // Other methods
78   - double semiminor(void) const { return m_semiminor.value(); }
79   - double semimajor(void) const { return m_semimajor.value(); }
80   - void semiminor(const double& semiminor) { m_semiminor.value(semiminor); }
81   - void semimajor(const double& semimajor) { m_semimajor.value(semimajor); }
  84 + double semiminor(void) const;
  85 + double semimajor(void) const;
  86 + void semiminor(const double& semiminor);
  87 + void semimajor(const double& semimajor);
82 88  
83 89 protected:
84 90 // Protected methods
... ... @@ -99,4 +105,76 @@ protected:
99 105 mutable double m_norm; //!< Normalization
100 106 };
101 107  
  108 +
  109 +/***********************************************************************//**
  110 + * @brief Return model type
  111 + *
  112 + * @return "EllipticalDisk".
  113 + *
  114 + * Returns the type of the elliptical disk model.
  115 + ***************************************************************************/
  116 +inline
  117 +std::string GModelSpatialEllipticalDisk::type(void) const
  118 +{
  119 + return "EllipticalDisk";
  120 +}
  121 +
  122 +
  123 +/***********************************************************************//**
  124 + * @brief Return semi-minor axis of ellipse
  125 + *
  126 + * @return Semi-minor axis of ellipse (degrees).
  127 + *
  128 + * Returns the semi-minor axis of the ellipse in degrees.
  129 + ***************************************************************************/
  130 +inline
  131 +double GModelSpatialEllipticalDisk::semiminor(void) const
  132 +{
  133 + return (m_semiminor.value());
  134 +}
  135 +
  136 +
  137 +/***********************************************************************//**
  138 + * @brief Set semi-minor axis of ellipse
  139 + *
  140 + * @param[in] semiminor Semi-minor axis of ellipse (degrees)
  141 + *
  142 + * Sets the semi-minor axis of the ellipse in degrees.
  143 + ***************************************************************************/
  144 +inline
  145 +void GModelSpatialEllipticalDisk::semiminor(const double& semiminor)
  146 +{
  147 + m_semiminor.value(semiminor);
  148 + return;
  149 +}
  150 +
  151 +
  152 +/***********************************************************************//**
  153 + * @brief Return semi-major axis of ellipse
  154 + *
  155 + * @return Semi-major axis of ellipse (degrees).
  156 + *
  157 + * Returns the semi-major axis of the ellipse in degrees.
  158 + ***************************************************************************/
  159 +inline
  160 +double GModelSpatialEllipticalDisk::semimajor(void) const
  161 +{
  162 + return (m_semimajor.value());
  163 +}
  164 +
  165 +
  166 +/***********************************************************************//**
  167 + * @brief Set semi-major axis of ellipse
  168 + *
  169 + * @param[in] semimajor Semi-major axis of ellipse (degrees)
  170 + *
  171 + * Sets the semi-major axis of the ellipse in degrees.
  172 + ***************************************************************************/
  173 +inline
  174 +void GModelSpatialEllipticalDisk::semimajor(const double& semimajor)
  175 +{
  176 + m_semimajor.value(semimajor);
  177 + return;
  178 +}
  179 +
102 180 #endif /* GMODELSPATIALELLIPTICALDISK_HPP */
... ...
include/GModelSpatialPointSource.hpp
  View file @cc742ac
... ... @@ -62,17 +62,19 @@ public:
62 62 // Implemented pure virtual methods
63 63 virtual void clear(void);
64 64 virtual GModelSpatialPointSource* clone(void) const;
65   - virtual std::string type(void) const { return "SkyDirFunction"; }
66   - virtual double eval(const GSkyDir& srcDir) const;
67   - virtual double eval_gradients(const GSkyDir& srcDir) const;
68   - virtual GSkyDir mc(GRan& ran) const;
  65 + virtual std::string type(void) const;
  66 + virtual double eval(const GPhoton& photon) const;
  67 + virtual double eval_gradients(const GPhoton& photon) const;
  68 + virtual GSkyDir mc(const GEnergy& energy,
  69 + const GTime& time,
  70 + GRan& ran) const;
69 71 virtual void read(const GXmlElement& xml);
70 72 virtual void write(GXmlElement& xml) const;
71 73 virtual std::string print(void) const;
72 74  
73 75 // Other methods
74   - double ra(void) const { return m_ra.value(); }
75   - double dec(void) const { return m_dec.value(); }
  76 + double ra(void) const;
  77 + double dec(void) const;
76 78 GSkyDir dir(void) const;
77 79 void dir(const GSkyDir& dir);
78 80  
... ... @@ -87,4 +89,47 @@ protected:
87 89 GModelPar m_dec; //!< Declination (deg)
88 90 };
89 91  
  92 +
  93 +/***********************************************************************//**
  94 + * @brief Return model type
  95 + *
  96 + * @return "SkyDirFunction".
  97 + *
  98 + * Returns the type of the spatial model.
  99 + ***************************************************************************/
  100 +inline
  101 +std::string GModelSpatialPointSource::type(void) const
  102 +{
  103 + return "SkyDirFunction";
  104 +}
  105 +
  106 +
  107 +/***********************************************************************//**
  108 + * @brief Return Right Ascencion of model centre
  109 + *
  110 + * @return Right Ascencion of model centre (degrees).
  111 + *
  112 + * Returns the Right Ascension of the model centre in degrees.
  113 + ***************************************************************************/
  114 +inline
  115 +double GModelSpatialPointSource::ra(void) const
  116 +{
  117 + return (m_ra.value());
  118 +}
  119 +
  120 +
  121 +/***********************************************************************//**
  122 + * @brief Return Declination of model centre
  123 + *
  124 + * @return Declination of model centre (degrees).
  125 + *
  126 + * Returns the Declination of the model centre in degrees.
  127 + ***************************************************************************/
  128 +inline
  129 +double GModelSpatialPointSource::dec(void) const
  130 +{
  131 + return (m_dec.value());
  132 +}
  133 +
  134 +
90 135 #endif /* GMODELSPATIALPOINTSOURCE_HPP */
... ...
include/GModelSpatialRadial.hpp
  View file @cc742ac
... ... @@ -33,6 +33,9 @@
33 33 #include "GModelSpatial.hpp"
34 34 #include "GModelPar.hpp"
35 35 #include "GSkyDir.hpp"
  36 +#include "GEnergy.hpp"
  37 +#include "GTime.hpp"
  38 +#include "GPhoton.hpp"
36 39 #include "GXmlElement.hpp"
37 40 #include "GRan.hpp"
38 41  
... ... @@ -62,21 +65,27 @@ public:
62 65 virtual void clear(void) = 0;
63 66 virtual GModelSpatialRadial* clone(void) const = 0;
64 67 virtual std::string type(void) const = 0;
65   - virtual double eval(const double& theta) const = 0;
66   - virtual double eval_gradients(const double& theta) const = 0;
67   - virtual GSkyDir mc(GRan& ran) const = 0;
  68 + virtual double eval(const double& theta,
  69 + const GEnergy& energy,
  70 + const GTime& time) const = 0;
  71 + virtual double eval_gradients(const double& theta,
  72 + const GEnergy& energy,
  73 + const GTime& time) const = 0;
  74 + virtual GSkyDir mc(const GEnergy& energy,
  75 + const GTime& time,
  76 + GRan& ran) const = 0;
68 77 virtual double theta_max(void) const = 0;
69 78 virtual std::string print(void) const = 0;
70 79  
71   - // Implemented virtual methods
72   - virtual double eval(const GSkyDir& srcDir) const;
73   - virtual double eval_gradients(const GSkyDir& srcDir) const;
  80 + // Implemented virtual base class methods
  81 + virtual double eval(const GPhoton& photon) const;
  82 + virtual double eval_gradients(const GPhoton& photon) const;
74 83 virtual void read(const GXmlElement& xml);
75 84 virtual void write(GXmlElement& xml) const;
76 85  
77 86 // Other methods
78   - double ra(void) const { return m_ra.value(); }
79   - double dec(void) const { return m_dec.value(); }
  87 + double ra(void) const;
  88 + double dec(void) const;
80 89 GSkyDir dir(void) const;
81 90 void dir(const GSkyDir& dir);
82 91  
... ... @@ -91,4 +100,32 @@ protected:
91 100 GModelPar m_dec; //!< Declination (deg)
92 101 };
93 102  
  103 +
  104 +/***********************************************************************//**
  105 + * @brief Return Right Ascencion of model centre
  106 + *
  107 + * @return Right Ascencion of model centre (degrees).
  108 + *
  109 + * Returns the Right Ascension of the model centre in degrees.
  110 + ***************************************************************************/
  111 +inline
  112 +double GModelSpatialRadial::ra(void) const
  113 +{
  114 + return (m_ra.value());
  115 +}
  116 +
  117 +
  118 +/***********************************************************************//**
  119 + * @brief Return Declination of model centre
  120 + *
  121 + * @return Declination of model centre (degrees).
  122 + *
  123 + * Returns the Declination of the model centre in degrees.
  124 + ***************************************************************************/
  125 +inline
  126 +double GModelSpatialRadial::dec(void) const
  127 +{
  128 + return (m_dec.value());
  129 +}
  130 +
94 131 #endif /* GMODELSPATIALRADIAL_HPP */
... ...
include/GModelSpatialRadialDisk.hpp
  View file @cc742ac
... ... @@ -57,21 +57,27 @@ public:
57 57 // Operators
58 58 virtual GModelSpatialRadialDisk& operator=(const GModelSpatialRadialDisk& model);
59 59  
60   - // Implemented pure virtual methods
  60 + // Implemented pure virtual base class methods
61 61 virtual void clear(void);
62 62 virtual GModelSpatialRadialDisk* clone(void) const;
63   - virtual std::string type(void) const { return "DiskFunction"; }
64   - virtual double eval(const double& theta) const;
65   - virtual double eval_gradients(const double& theta) const;
66   - virtual GSkyDir mc(GRan& ran) const;
  63 + virtual std::string type(void) const;
  64 + virtual double eval(const double& theta,
  65 + const GEnergy& energy,
  66 + const GTime& time) const;
  67 + virtual double eval_gradients(const double& theta,
  68 + const GEnergy& energy,
  69 + const GTime& time) const;
  70 + virtual GSkyDir mc(const GEnergy& energy,
  71 + const GTime& time,
  72 + GRan& ran) const;
67 73 virtual double theta_max(void) const;
68 74 virtual void read(const GXmlElement& xml);
69 75 virtual void write(GXmlElement& xml) const;
70 76 virtual std::string print(void) const;
71 77  
72 78 // Other methods
73   - double radius(void) const { return m_radius.value(); }
74   - void radius(const double& radius) { m_radius.value(radius); }
  79 + double radius(void) const;
  80 + void radius(const double& radius);
75 81  
76 82 protected:
77 83 // Protected methods
... ... @@ -89,4 +95,47 @@ protected:
89 95 mutable double m_norm; //!< Normalization
90 96 };
91 97  
  98 +
  99 +/***********************************************************************//**
  100 + * @brief Return model type
  101 + *
  102 + * @return "DiskFunction".
  103 + *
  104 + * Returns the type of the radial disk model.
  105 + ***************************************************************************/
  106 +inline
  107 +std::string GModelSpatialRadialDisk::type(void) const
  108 +{
  109 + return "DiskFunction";
  110 +}
  111 +
  112 +
  113 +/***********************************************************************//**
  114 + * @brief Return disk radius
  115 + *
  116 + * @return Disk radius (degrees).
  117 + *
  118 + * Returns the radius of the disk in degrees.
  119 + ***************************************************************************/
  120 +inline
  121 +double GModelSpatialRadialDisk::radius(void) const
  122 +{
  123 + return (m_radius.value());
  124 +}
  125 +
  126 +
  127 +/***********************************************************************//**
  128 + * @brief Set disk radius
  129 + *
  130 + * @param[in] radius Disk radius (degrees).
  131 + *
  132 + * Sets the radius of the disk in degrees.
  133 + ***************************************************************************/
  134 +inline
  135 +void GModelSpatialRadialDisk::radius(const double& radius)
  136 +{
  137 + m_radius.value(radius);
  138 + return;
  139 +}
  140 +
92 141 #endif /* GMODELSPATIALRADIALDISK_HPP */
... ...
include/GModelSpatialRadialGauss.hpp
  View file @cc742ac
... ... @@ -32,6 +32,8 @@
32 32 #include "GModelSpatialRadial.hpp"
33 33 #include "GModelPar.hpp"
34 34 #include "GSkyDir.hpp"
  35 +#include "GEnergy.hpp"
  36 +#include "GTime.hpp"
35 37 #include "GXmlElement.hpp"
36 38  
37 39  
... ... @@ -59,18 +61,24 @@ public:
59 61 // Implemented pure virtual methods
60 62 virtual void clear(void);
61 63 virtual GModelSpatialRadialGauss* clone(void) const;
62   - virtual std::string type(void) const { return "GaussFunction"; }
63   - virtual double eval(const double& theta) const;
64   - virtual double eval_gradients(const double& theta) const;
65   - virtual GSkyDir mc(GRan& ran) const;
  64 + virtual std::string type(void) const;
  65 + virtual double eval(const double& theta,
  66 + const GEnergy& energy,
  67 + const GTime& time) const;
  68 + virtual double eval_gradients(const double& theta,
  69 + const GEnergy& energy,
  70 + const GTime& time) const;
  71 + virtual GSkyDir mc(const GEnergy& energy,
  72 + const GTime& time,
  73 + GRan& ran) const;
66 74 virtual double theta_max(void) const;
67 75 virtual void read(const GXmlElement& xml);
68 76 virtual void write(GXmlElement& xml) const;
69 77 virtual std::string print(void) const;
70 78  
71 79 // Other methods
72   - double sigma(void) const { return m_sigma.value(); }
73   - void sigma(const double& sigma) { m_sigma.value(sigma); }
  80 + double sigma(void) const;
  81 + void sigma(const double& sigma);
74 82  
75 83 protected:
76 84 // Protected methods
... ... @@ -82,4 +90,47 @@ protected:
82 90 GModelPar m_sigma; //!< Gaussian width (deg)
83 91 };
84 92  
  93 +
  94 +/***********************************************************************//**
  95 + * @brief Return model type
  96 + *
  97 + * @return "GaussFunction".
  98 + *
  99 + * Returns the type of the radial Gauss model.
  100 + ***************************************************************************/
  101 +inline
  102 +std::string GModelSpatialRadialGauss::type(void) const
  103 +{
  104 + return "GaussFunction";
  105 +}
  106 +
  107 +
  108 +/***********************************************************************//**
  109 + * @brief Return Gaussian sigma
  110 + *
  111 + * @return Gaussian sigma (degrees).
  112 + *
  113 + * Returns the Gaussian sigma in degrees.
  114 + ***************************************************************************/
  115 +inline
  116 +double GModelSpatialRadialGauss::sigma(void) const
  117 +{
  118 + return (m_sigma.value());
  119 +}
  120 +
  121 +
  122 +/***********************************************************************//**
  123 + * @brief Set Gaussian sigma
  124 + *
  125 + * @param[in] sigma Gaussian sigma (degrees).
  126 + *
  127 + * Sets the Gaussian sigma in degrees.
  128 + ***************************************************************************/
  129 +inline
  130 +void GModelSpatialRadialGauss::sigma(const double& sigma)
  131 +{
  132 + m_sigma.value(sigma);
  133 + return;
  134 +}
  135 +
85 136 #endif /* GMODELSPATIALRADIALGAUSS_HPP */
... ...
include/GModelSpatialRadialShell.hpp
  View file @cc742ac
... ... @@ -72,22 +72,28 @@ public:
72 72 // Implemented pure virtual methods
73 73 virtual void clear(void);
74 74 virtual GModelSpatialRadialShell* clone(void) const;
75   - virtual std::string type(void) const { return "ShellFunction"; }
76   - virtual double eval(const double& theta) const;
77   - virtual double eval_gradients(const double& theta) const;
78   - virtual GSkyDir mc(GRan& ran) const;
  75 + virtual std::string type(void) const;
  76 + virtual double eval(const double& theta,
  77 + const GEnergy& energy,
  78 + const GTime& time) const;
  79 + virtual double eval_gradients(const double& theta,
  80 + const GEnergy& energy,
  81 + const GTime& time) const;
  82 + virtual GSkyDir mc(const GEnergy& energy,
  83 + const GTime& time,
  84 + GRan& ran) const;
79 85 virtual double theta_max(void) const;
80 86 virtual void read(const GXmlElement& xml);
81 87 virtual void write(GXmlElement& xml) const;
82 88 virtual std::string print(void) const;
83 89  
84 90 // Other methods
85   - double radius(void) const { return m_radius.value(); }
86   - double width(void) const { return m_width.value(); }
87   - bool small_angle(void) const { return m_small_angle; }
88   - void radius(const double& radius) { m_radius.value(radius); }
89   - void width(const double& width) { m_width.value(width); }
90   - void small_angle(const bool& small_angle) { m_small_angle = small_angle; }
  91 + double radius(void) const;
  92 + double width(void) const;
  93 + const bool& small_angle(void) const;
  94 + void radius(const double& radius);
  95 + void width(const double& width);
  96 + void small_angle(const bool& small_angle);
91 97  
92 98 protected:
93 99 // Protected methods
... ... @@ -113,4 +119,107 @@ protected:
113 119 mutable double m_norm; //!< Shell normalization
114 120 };
115 121  
  122 +
  123 +/***********************************************************************//**
  124 + * @brief Return model type
  125 + *
  126 + * @return "ShellFunction".
  127 + *
  128 + * Returns the type of the radial shell model.
  129 + ***************************************************************************/
  130 +inline
  131 +std::string GModelSpatialRadialShell::type(void) const
  132 +{
  133 + return "ShellFunction";
  134 +}
  135 +
  136 +
  137 +/***********************************************************************//**
  138 + * @brief Return shell radius
  139 + *
  140 + * @return Shell radius (degrees).
  141 + *
  142 + * Returns the shell radius in degrees.
  143 + ***************************************************************************/
  144 +inline
  145 +double GModelSpatialRadialShell::radius(void) const
  146 +{
  147 + return (m_radius.value());
  148 +}
  149 +
  150 +
  151 +/***********************************************************************//**
  152 + * @brief Set shell radius
  153 + *
  154 + * @param[in] radius Shell radius (degrees).
  155 + *
  156 + * Sets the shell radius in degrees.
  157 + ***************************************************************************/
  158 +inline
  159 +void GModelSpatialRadialShell::radius(const double& radius)
  160 +{
  161 + m_radius.value(radius);
  162 + return;
  163 +}
  164 +
  165 +
  166 +/***********************************************************************//**
  167 + * @brief Return shell width
  168 + *
  169 + * @return Shell width (degrees).
  170 + *
  171 + * Returns the shell width in degrees.
  172 + ***************************************************************************/
  173 +inline
  174 +double GModelSpatialRadialShell::width(void) const
  175 +{
  176 + return (m_width.value());
  177 +}
  178 +
  179 +
  180 +/***********************************************************************//**
  181 + * @brief Set width radius
  182 + *
  183 + * @param[in] width Shell width (degrees).
  184 + *
  185 + * Sets the shell width in degrees.
  186 + ***************************************************************************/
  187 +inline
  188 +void GModelSpatialRadialShell::width(const double& width)
  189 +{
  190 + m_width.value(width);
  191 + return;
  192 +}
  193 +
  194 +
  195 +/***********************************************************************//**
  196 + * @brief Return shell approximation flag
  197 + *
  198 + * @return True if small angle approximation is taken.
  199 + *
  200 + * Signals if the small angle approximation is taken.
  201 + ***************************************************************************/
  202 +inline
  203 +const bool& GModelSpatialRadialShell::small_angle(void) const
  204 +{
  205 + return (m_small_angle);
  206 +}
  207 +
  208 +
  209 +/***********************************************************************//**
  210 + * @brief Set shell approximation flag
  211 + *
  212 + * @param[in] small_angle Shell approximation flag.
  213 + *
  214 + * Sets the shell approximation flag. True specifies that the small angle
  215 + * approximation shall be taken, false specifies that the correct spherical
  216 + * computation shall be used.
  217 + ***************************************************************************/
  218 +inline
  219 +void GModelSpatialRadialShell::small_angle(const bool& small_angle)
  220 +{
  221 + m_small_angle = small_angle;
  222 + return;
  223 +}
  224 +
116 225 #endif /* GMODELSPATIALRADIALSHELL_HPP */
... ...
include/GResponse.hpp
  View file @cc742ac
... ... @@ -117,12 +117,12 @@ protected:
117 117 // Npred theta integration kernel for radial model
118 118 class npred_radial_kern_theta : public GFunction {
119 119 public:
120   - npred_radial_kern_theta(const GResponse* rsp,
121   - const GModelSpatialRadial* spatial,
122   - const GEnergy* srcEng,
123   - const GTime* srcTime,
124   - const GObservation* obs,
125   - const GMatrix* rot) :
  120 + npred_radial_kern_theta(const GResponse& rsp,
  121 + const GModelSpatialRadial& spatial,
  122 + const GEnergy& srcEng,
  123 + const GTime& srcTime,
  124 + const GObservation& obs,
  125 + const GMatrix& rot) :
126 126 m_rsp(rsp),
127 127 m_spatial(spatial),
128 128 m_srcEng(srcEng),
... ... @@ -131,24 +131,24 @@ protected:
131 131 m_rot(rot) { }
132 132 double eval(double theta);
133 133 protected:
134   - const GResponse* m_rsp; //!< Pointer to response
135   - const GModelSpatialRadial* m_spatial; //!< Pointer to spatial model
136   - const GEnergy* m_srcEng; //!< Pointer to true photon energy
137   - const GTime* m_srcTime; //!< Pointer to true photon arrival time
138   - const GObservation* m_obs; //!< Pointer to observation
139   - const GMatrix* m_rot; //!< Pointer to rotation matrix
  134 + const GResponse& m_rsp; //!< Response
  135 + const GModelSpatialRadial& m_spatial; //!< Spatial model
  136 + const GEnergy& m_srcEng; //!< True photon energy
  137 + const GTime& m_srcTime; //!< True photon arrival time
  138 + const GObservation& m_obs; //!< Observation
  139 + const GMatrix& m_rot; //!< Rotation matrix
140 140 };
141 141  
142 142 // Npred phi integration kernel for radial model
143 143 class npred_radial_kern_phi : public GFunction {
144 144 public:
145   - npred_radial_kern_phi(const GResponse* rsp,
146   - const GEnergy* srcEng,
147   - const GTime* srcTime,
148   - const GObservation* obs,
149   - const GMatrix* rot,
150   - double theta,
151   - double sin_theta) :
  145 + npred_radial_kern_phi(const GResponse& rsp,
  146 + const GEnergy& srcEng,
  147 + const GTime& srcTime,
  148 + const GObservation& obs,
  149 + const GMatrix& rot,
  150 + const double& theta,
  151 + const double& sin_theta) :
152 152 m_rsp(rsp),
153 153 m_srcEng(srcEng),
154 154 m_srcTime(srcTime),
... ... @@ -159,25 +159,25 @@ protected:
159 159 m_sin_theta(sin_theta) { }
160 160 double eval(double phi);
161 161 protected:
162   - const GResponse* m_rsp; //!< Pointer to response
163   - const GEnergy* m_srcEng; //!< Pointer to true photon energy
164   - const GTime* m_srcTime; //!< Pointer to true photon arrival time
165   - const GObservation* m_obs; //!< Pointer to observation
166   - const GMatrix* m_rot; //!< Pointer to rotation matrix
167   - double m_theta; //!< Offset angle (radians)
168   - double m_cos_theta; //!< cosine of offset angle
169   - double m_sin_theta; //!< Sine of offset angle
  162 + const GResponse& m_rsp; //!< Response
  163 + const GEnergy& m_srcEng; //!< True photon energy
  164 + const GTime& m_srcTime; //!< True photon arrival time
  165 + const GObservation& m_obs; //!< Observation
  166 + const GMatrix& m_rot; //!< Rotation matrix
  167 + const double& m_theta; //!< Offset angle (radians)
  168 + double m_cos_theta; //!< cosine of offset angle
  169 + const double& m_sin_theta; //!< Sine of offset angle
170 170 };
171 171  
172 172 // Npred theta integration kernel for elliptical model
173 173 class npred_elliptical_kern_theta : public GFunction {
174 174 public:
175   - npred_elliptical_kern_theta(const GResponse* rsp,
176   - const GModelSpatialElliptical* spatial,
177   - const GEnergy* srcEng,
178   - const GTime* srcTime,
179   - const GObservation* obs,
180   - const GMatrix* rot) :
  175 + npred_elliptical_kern_theta(const GResponse& rsp,
  176 + const GModelSpatialElliptical& spatial,
  177 + const GEnergy& srcEng,
  178 + const GTime& srcTime,
  179 + const GObservation& obs,
  180 + const GMatrix& rot) :
181 181 m_rsp(rsp),
182 182 m_spatial(spatial),
183 183 m_srcEng(srcEng),
... ... @@ -186,25 +186,25 @@ protected:
186 186 m_rot(rot) { }
187 187 double eval(double theta);
188 188 protected:
189   - const GResponse* m_rsp; //!< Pointer to response
190   - const GModelSpatialElliptical* m_spatial; //!< Pointer to spatial model
191   - const GEnergy* m_srcEng; //!< Pointer to true photon energy
192   - const GTime* m_srcTime; //!< Pointer to true photon arrival time
193   - const GObservation* m_obs; //!< Pointer to observation
194   - const GMatrix* m_rot; //!< Pointer to rotation matrix
  189 + const GResponse& m_rsp; //!< Response
  190 + const GModelSpatialElliptical& m_spatial; //!< Spatial model
  191 + const GEnergy& m_srcEng; //!< True photon energy
  192 + const GTime& m_srcTime; //!< True photon arrival time
  193 + const GObservation& m_obs; //!< Observation
  194 + const GMatrix& m_rot; //!< Rotation matrix
195 195 };
196 196  
197 197 // Npred phi integration kernel for elliptical model
198 198 class npred_elliptical_kern_phi : public GFunction {
199 199 public:
200   - npred_elliptical_kern_phi(const GResponse* rsp,
201   - const GModelSpatialElliptical* spatial,
202   - const GEnergy* srcEng,
203   - const GTime* srcTime,
204   - const GObservation* obs,
205   - const GMatrix* rot,
206   - double theta,
207   - double sin_theta) :
  200 + npred_elliptical_kern_phi(const GResponse& rsp,
  201 + const GModelSpatialElliptical& spatial,
  202 + const GEnergy& srcEng,
  203 + const GTime& srcTime,
  204 + const GObservation& obs,
  205 + const GMatrix& rot,
  206 + const double& theta,
  207 + const double& sin_theta) :
208 208 m_rsp(rsp),
209 209 m_spatial(spatial),
210 210 m_srcEng(srcEng),
... ... @@ -216,15 +216,15 @@ protected:
216 216 m_sin_theta(sin_theta) { }
217 217 double eval(double phi);
218 218 protected:
219   - const GResponse* m_rsp; //!< Pointer to response
220   - const GModelSpatialElliptical* m_spatial; //!< Pointer to spatial model
221   - const GEnergy* m_srcEng; //!< Pointer to true photon energy
222   - const GTime* m_srcTime; //!< Pointer to true photon arrival time
223   - const GObservation* m_obs; //!< Pointer to observation
224   - const GMatrix* m_rot; //!< Pointer to rotation matrix
225   - double m_theta; //!< Offset angle (radians)
  219 + const GResponse& m_rsp; //!< Response
  220 + const GModelSpatialElliptical& m_spatial; //!< Spatial model
  221 + const GEnergy& m_srcEng; //!< True photon energy
  222 + const GTime& m_srcTime; //!< True photon arrival time
  223 + const GObservation& m_obs; //!< Observation
  224 + const GMatrix& m_rot; //!< Rotation matrix
  225 + const double& m_theta; //!< Offset angle (radians)
226 226 double m_cos_theta; //!< cosine of offset angle
227   - double m_sin_theta; //!< Sine of offset angle
  227 + const double& m_sin_theta; //!< Sine of offset angle
228 228 };
229 229  
230 230 };
... ...
pyext/GModelSpatial.i
  View file @cc742ac
... ... @@ -26,6 +26,7 @@
26 26 %{
27 27 /* Put headers and other declarations here that are needed for compilation */
28 28 #include "GModelSpatial.hpp"
  29 +#include "GTools.hpp"
29 30 %}
30 31  
31 32  
... ... @@ -34,7 +35,7 @@
34 35 *
35 36 * @brief Abstract spatial model base class
36 37 ***************************************************************************/
37   -class GModelSpatial {
  38 +class GModelSpatial : public GBase {
38 39  
39 40 public:
40 41 // Constructors and destructors
... ... @@ -46,9 +47,10 @@ public:
46 47 virtual void clear(void) = 0;
47 48 virtual GModelSpatial* clone(void) const = 0;
48 49 virtual std::string type(void) const = 0;
49   - virtual double eval(const GSkyDir& srcDir) const = 0;
50   - virtual double eval_gradients(const GSkyDir& srcDir) const = 0;
51   - virtual GSkyDir mc(GRan& ran) const = 0;
  50 + virtual double eval(const GPhoton& photon) const = 0;
  51 + virtual double eval_gradients(const GPhoton& photon) const = 0;
  52 + virtual GSkyDir mc(const GEnergy& energy, const GTime& time,
  53 + GRan& ran) const = 0;
52 54 virtual void read(const GXmlElement& xml) = 0;
53 55 virtual void write(GXmlElement& xml) const = 0;
54 56  
... ...
pyext/GModelSpatialDiffuse.i
  View file @cc742ac
... ... @@ -45,9 +45,10 @@ public:
45 45 virtual void clear(void) = 0;
46 46 virtual GModelSpatialDiffuse* clone(void) const = 0;
47 47 virtual std::string type(void) const = 0;
48   - virtual double eval(const GSkyDir& srcDir) const = 0;
49   - virtual double eval_gradients(const GSkyDir& srcDir) const = 0;
50   - virtual GSkyDir mc(GRan& ran) const = 0;
  48 + virtual double eval(const GPhoton& photon) const = 0;
  49 + virtual double eval_gradients(const GPhoton& photon) const = 0;
  50 + virtual GSkyDir mc(const GEnergy& energy, const GTime& time,
  51 + GRan& ran) const = 0;
51 52 virtual void read(const GXmlElement& xml) = 0;
52 53 virtual void write(GXmlElement& xml) const = 0;
53 54 };
... ...
pyext/GModelSpatialDiffuseConst.i
  View file @cc742ac
... ... @@ -48,9 +48,11 @@ public:
48 48 virtual void clear(void);
49 49 virtual GModelSpatialDiffuseConst* clone(void) const;
50 50 virtual std::string type(void) const;
51   - virtual double eval(const GSkyDir& srcDir) const;
52   - virtual double eval_gradients(const GSkyDir& srcDir) const;
53   - virtual GSkyDir mc(GRan& ran) const;
  51 + virtual double eval(const GPhoton& photon) const;
  52 + virtual double eval_gradients(const GPhoton& photon) const;
  53 + virtual GSkyDir mc(const GEnergy& energy,
  54 + const GTime& time,
  55 + GRan& ran) const;
54 56 virtual void read(const GXmlElement& xml);
55 57 virtual void write(GXmlElement& xml) const;
56 58  
... ...
pyext/GModelSpatialDiffuseCube.i
  View file @cc742ac
... ... @@ -51,9 +51,11 @@ public:
51 51 virtual void clear(void);
52 52 virtual GModelSpatialDiffuseCube* clone(void) const;
53 53 virtual std::string type(void) const;
54   - virtual double eval(const GSkyDir& srcDir) const;
55   - virtual double eval_gradients(const GSkyDir& srcDir) const;
56   - virtual GSkyDir mc(GRan& ran) const;
  54 + virtual double eval(const GPhoton& photon) const;
  55 + virtual double eval_gradients(const GPhoton& photon) const;
  56 + virtual GSkyDir mc(const GEnergy& energy,
  57 + const GTime& time,
  58 + GRan& ran) const;
57 59 virtual void read(const GXmlElement& xml);
58 60 virtual void write(GXmlElement& xml) const;
59 61  
... ...
pyext/GModelSpatialDiffuseMap.i
  View file @cc742ac
... ... @@ -48,9 +48,11 @@ public:
48 48 virtual void clear(void);
49 49 virtual GModelSpatialDiffuseMap* clone(void) const;
50 50 virtual std::string type(void) const;
51   - virtual double eval(const GSkyDir& srcDir) const;
52   - virtual double eval_gradients(const GSkyDir& srcDir) const;
53   - virtual GSkyDir mc(GRan& ran) const;
  51 + virtual double eval(const GPhoton& photon) const;
  52 + virtual double eval_gradients(const GPhoton& photon) const;
  53 + virtual GSkyDir mc(const GEnergy& energy,
  54 + const GTime& time,
  55 + GRan& ran) const;
54 56 virtual void read(const GXmlElement& xml);
55 57 virtual void write(GXmlElement& xml) const;
56 58 };
... ...
pyext/GModelSpatialElliptical.i
  View file @cc742ac
... ... @@ -46,16 +46,22 @@ public:
46 46 virtual void clear(void) = 0;
47 47 virtual GModelSpatialElliptical* clone(void) const = 0;
48 48 virtual std::string type(void) const = 0;
49   - virtual double eval(const double& theta,
50   - const double& posangle) const = 0;
51   - virtual double eval_gradients(const double& theta,
52   - const double& posangle) const = 0;
53   - virtual GSkyDir mc(GRan& ran) const = 0;
  49 + virtual double eval(const double& theta,
  50 + const double& posangle,
  51 + const GEnergy& energy,
  52 + const GTime& time) const = 0;
  53 + virtual double eval_gradients(const double& theta,
  54 + const double& posangle,
  55 + const GEnergy& energy,
  56 + const GTime& time) const = 0;
  57 + virtual GSkyDir mc(const GEnergy& energy,
  58 + const GTime& time,
  59 + GRan& ran) const = 0;
54 60 virtual double theta_max(void) const = 0;
55 61  
56 62 // Implemented virtual methods
57   - virtual double eval(const GSkyDir& srcDir) const;
58   - virtual double eval_gradients(const GSkyDir& srcDir) const;
  63 + virtual double eval(const GPhoton& photon) const;
  64 + virtual double eval_gradients(const GPhoton& photon) const;
59 65 virtual void read(const GXmlElement& xml);
60 66 virtual void write(GXmlElement& xml) const;
61 67  
... ...
pyext/GModelSpatialEllipticalDisk.i
  View file @cc742ac
... ... @@ -51,11 +51,17 @@ public:
51 51 virtual void clear(void);
52 52 virtual GModelSpatialEllipticalDisk* clone(void) const;
53 53 virtual std::string type(void) const;
54   - virtual double eval(const double& theta,
55   - const double& posangle) const;
56   - virtual double eval_gradients(const double& theta,
57   - const double& posangle) const;
58   - virtual GSkyDir mc(GRan& ran) const;
  54 + virtual double eval(const double& theta,
  55 + const double& posangle,
  56 + const GEnergy& energy,
  57 + const GTime& time) const;
  58 + virtual double eval_gradients(const double& theta,
  59 + const double& posangle,
  60 + const GEnergy& energy,
  61 + const GTime& time) const;
  62 + virtual GSkyDir mc(const GEnergy& energy,
  63 + const GTime& time,
  64 + GRan& ran) const;
59 65 virtual double theta_max(void) const;
60 66 virtual void read(const GXmlElement& xml);
61 67 virtual void write(GXmlElement& xml) const;
... ... @@ -75,10 +81,10 @@ public:
75 81 GModelSpatialEllipticalDisk copy() {
76 82 return (*self);
77 83 }
78   - double eval(const GSkyDir& srcDir) const {
79   - return self->GModelSpatialElliptical::eval(srcDir);
  84 + double eval(const GPhoton& photon) const {
  85 + return self->GModelSpatialElliptical::eval(photon);
80 86 }
81   - double eval_gradients(const GSkyDir& srcDir) const {
82   - return self->GModelSpatialElliptical::eval_gradients(srcDir);
  87 + double eval_gradients(const GPhoton& photon) const {
  88 + return self->GModelSpatialElliptical::eval_gradients(photon);
83 89 }
84 90 };
... ...
pyext/GModelSpatialPointSource.i
  View file @cc742ac
... ... @@ -48,9 +48,11 @@ public:
48 48 virtual void clear(void);
49 49 virtual GModelSpatialPointSource* clone(void) const;
50 50 virtual std::string type(void) const;
51   - virtual double eval(const GSkyDir& srcDir) const;
52   - virtual double eval_gradients(const GSkyDir& srcDir) const;
53   - virtual GSkyDir mc(GRan& ran) const;
  51 + virtual double eval(const GPhoton& photon) const;
  52 + virtual double eval_gradients(const GPhoton& photon) const;
  53 + virtual GSkyDir mc(const GEnergy& energy,
  54 + const GTime& time,
  55 + GRan& ran) const;
54 56 virtual void read(const GXmlElement& xml);
55 57 virtual void write(GXmlElement& xml) const;
56 58  
... ...
pyext/GModelSpatialRadial.i
  View file @cc742ac
... ... @@ -47,14 +47,20 @@ public:
47 47 virtual void clear(void) = 0;
48 48 virtual GModelSpatialRadial* clone(void) const = 0;
49 49 virtual std::string type(void) const = 0;
50   - virtual double eval(const double& theta) const = 0;
51   - virtual double eval_gradients(const double& theta) const = 0;
52   - virtual GSkyDir mc(GRan& ran) const = 0;
  50 + virtual double eval(const double& theta,
  51 + const GEnergy& energy,
  52 + const GTime& time) const = 0;
  53 + virtual double eval_gradients(const double& theta,
  54 + const GEnergy& energy,
  55 + const GTime& time) const = 0;
  56 + virtual GSkyDir mc(const GEnergy& energy,
  57 + const GTime& time,
  58 + GRan& ran) const = 0;
53 59 virtual double theta_max(void) const = 0;
54 60  
55   - // Implemented virtual methods
56   - virtual double eval(const GSkyDir& srcDir) const;
57   - virtual double eval_gradients(const GSkyDir& srcDir) const;
  61 + // Implemented virtual base class methods
  62 + virtual double eval(const GPhoton& photon) const;
  63 + virtual double eval_gradients(const GPhoton& photon) const;
58 64 virtual void read(const GXmlElement& xml);
59 65 virtual void write(GXmlElement& xml) const;
60 66  
... ...
pyext/GModelSpatialRadialDisk.i
  View file @cc742ac
... ... @@ -51,9 +51,15 @@ public:
51 51 virtual void clear(void);
52 52 virtual GModelSpatialRadialDisk* clone(void) const;
53 53 virtual std::string type(void) const;
54   - virtual double eval(const double& theta) const;
55   - virtual double eval_gradients(const double& theta) const;
56   - virtual GSkyDir mc(GRan& ran) const;
  54 + virtual double eval(const double& theta,
  55 + const GEnergy& energy,
  56 + const GTime& time) const;
  57 + virtual double eval_gradients(const double& theta,
  58 + const GEnergy& energy,
  59 + const GTime& time) const;
  60 + virtual GSkyDir mc(const GEnergy& energy,
  61 + const GTime& time,
  62 + GRan& ran) const;
57 63 virtual double theta_max(void) const;
58 64 virtual void read(const GXmlElement& xml);
59 65 virtual void write(GXmlElement& xml) const;
... ... @@ -77,10 +83,10 @@ public:
77 83 GModelSpatialRadialDisk copy() {
78 84 return (*self);
79 85 }
80   - double eval(const GSkyDir& srcDir) const {
81   - return self->GModelSpatialRadial::eval(srcDir);
  86 + double eval(const GPhoton& photon) const {
  87 + return self->GModelSpatialRadial::eval(photon);
82 88 }
83   - double eval_gradients(const GSkyDir& srcDir) const {
84   - return self->GModelSpatialRadial::eval_gradients(srcDir);
  89 + double eval_gradients(const GPhoton& photon) const {
  90 + return self->GModelSpatialRadial::eval_gradients(photon);
85 91 }
86 92 };
... ...
pyext/GModelSpatialRadialGauss.i
  View file @cc742ac
... ... @@ -48,9 +48,15 @@ public:
48 48 virtual void clear(void);
49 49 virtual GModelSpatialRadialGauss* clone(void) const;
50 50 virtual std::string type(void) const;
51   - virtual double eval(const double& theta) const;
52   - virtual double eval_gradients(const double& theta) const;
53   - virtual GSkyDir mc(GRan& ran) const;
  51 + virtual double eval(const double& theta,
  52 + const GEnergy& energy,
  53 + const GTime& time) const;
  54 + virtual double eval_gradients(const double& theta,
  55 + const GEnergy& energy,
  56 + const GTime& time) const;
  57 + virtual GSkyDir mc(const GEnergy& energy,
  58 + const GTime& time,
  59 + GRan& ran) const;
54 60 virtual double theta_max(void) const;
55 61 virtual void read(const GXmlElement& xml);
56 62 virtual void write(GXmlElement& xml) const;
... ... @@ -74,10 +80,10 @@ public:
74 80 GModelSpatialRadialGauss copy() {
75 81 return (*self);
76 82 }
77   - double eval(const GSkyDir& srcDir) const {
78   - return self->GModelSpatialRadial::eval(srcDir);
  83 + double eval(const GPhoton& photon) const {
  84 + return self->GModelSpatialRadial::eval(photon);
79 85 }
80   - double eval_gradients(const GSkyDir& srcDir) const {
81   - return self->GModelSpatialRadial::eval_gradients(srcDir);
  86 + double eval_gradients(const GPhoton& photon) const {
  87 + return self->GModelSpatialRadial::eval_gradients(photon);
82 88 }
83 89 };
... ...
pyext/GModelSpatialRadialShell.i
  View file @cc742ac
... ... @@ -49,9 +49,15 @@ public:
49 49 virtual void clear(void);
50 50 virtual GModelSpatialRadialShell* clone(void) const;
51 51 virtual std::string type(void) const;
52   - virtual double eval(const double& theta) const;
53   - virtual double eval_gradients(const double& theta) const;
54   - virtual GSkyDir mc(GRan& ran) const;
  52 + virtual double eval(const double& theta,
  53 + const GEnergy& energy,
  54 + const GTime& time) const;
  55 + virtual double eval_gradients(const double& theta,
  56 + const GEnergy& energy,
  57 + const GTime& time) const;
  58 + virtual GSkyDir mc(const GEnergy& energy,
  59 + const GTime& time,
  60 + GRan& ran) const;
55 61 virtual double theta_max(void) const;
56 62 virtual void read(const GXmlElement& xml);
57 63 virtual void write(GXmlElement& xml) const;
... ... @@ -79,10 +85,10 @@ public:
79 85 GModelSpatialRadialShell copy() {
80 86 return (*self);
81 87 }
82   - double eval(const GSkyDir& srcDir) const {
83   - return self->GModelSpatialRadial::eval(srcDir);
  88 + double eval(const GPhoton& photon) const {
  89 + return self->GModelSpatialRadial::eval(photon);
84 90 }
85   - double eval_gradients(const GSkyDir& srcDir) const {
86   - return self->GModelSpatialRadial::eval_gradients(srcDir);
  91 + double eval_gradients(const GPhoton& photon) const {
  92 + return self->GModelSpatialRadial::eval_gradients(photon);
87 93 }
88 94 };
... ...
src/model/GModelSky.cpp
  View file @cc742ac
... ... @@ -414,7 +414,7 @@ double GModelSky::value(const GPhoton&amp; photon)
414 414 double source = 1.0;
415 415  
416 416 // Evaluate source model
417   - if (m_spatial != NULL) source *= m_spatial->eval(photon.dir());
  417 + if (m_spatial != NULL) source *= m_spatial->eval(photon);
418 418 if (m_spectral != NULL) source *= m_spectral->eval(photon.energy());
419 419 if (m_temporal != NULL) source *= m_temporal->eval(photon.time());
420 420  
... ... @@ -436,7 +436,7 @@ double GModelSky::value(const GPhoton&amp; photon)
436 436 GVector GModelSky::gradients(const GPhoton& photon)
437 437 {
438 438 // Evaluate source model gradients
439   - if (m_spatial != NULL) m_spatial->eval_gradients(photon.dir());
  439 + if (m_spatial != NULL) m_spatial->eval_gradients(photon);
440 440 if (m_spectral != NULL) m_spectral->eval_gradients(photon.energy());
441 441 if (m_temporal != NULL) m_temporal->eval_gradients(photon.time());
442 442  
... ... @@ -787,6 +787,8 @@ void GModelSky::write(GXmlElement&amp; xml) const
787 787 * computations.
788 788 * @todo Implement photon arrival direction simulation for diffuse models
789 789 * @todo Implement unique model ID to assign as Monte Carlo ID
  790 + *
  791 + * @todo THIS METHOD SO FAR ONLY WORKS FOR FACTORIZED SOURCE MODELS!!!!!
790 792 ***************************************************************************/
791 793 GPhotons GModelSky::mc(const double& area,
792 794 const GSkyDir& dir, const double& radius,
... ... @@ -853,11 +855,12 @@ GPhotons GModelSky::mc(const double&amp; area,
853 855 // Set photon arrival time
854 856 photon.time(times[i]);
855 857  
856   - // Set incident photon direction
857   - photon.dir(m_spatial->mc(ran));
858   -
859 858 // Set photon energy
860 859 photon.energy(m_spectral->mc(emin, emax, ran));
  860 + //photon.energy(m_spectral->mc(emin, emax, photon.time(), ran));
  861 +
  862 + // Set incident photon direction
  863 + photon.dir(m_spatial->mc(photon.energy(), photon.time(), ran));
861 864  
862 865 // Append photon
863 866 photons.append(photon);
... ...
src/model/GModelSpatialDiffuseConst.cpp
  View file @cc742ac
... ... @@ -229,12 +229,13 @@ GModelSpatialDiffuseConst* GModelSpatialDiffuseConst::clone(void) const
229 229 /***********************************************************************//**
230 230 * @brief Evaluate function
231 231 *
232   - * @param[in] srcDir True photon arrival direction.
  232 + * @param[in] photon Incident photon (ignored).
  233 + * @return Model value.
233 234 *
234 235 * Evaluates the spatial part for an isotropic source model. By definition
235 236 * this value is independent from the sky direction.
236 237 ***************************************************************************/
237   -double GModelSpatialDiffuseConst::eval(const GSkyDir& srcDir) const
  238 +double GModelSpatialDiffuseConst::eval(const GPhoton& photon) const
238 239 {
239 240 // Return value
240 241 return (m_value.value());
... ... @@ -244,13 +245,14 @@ double GModelSpatialDiffuseConst::eval(const GSkyDir&amp; srcDir) const
244 245 /***********************************************************************//**
245 246 * @brief Evaluate function and gradients
246 247 *
247   - * @param[in] srcDir True photon arrival direction.
  248 + * @param[in] photon Incident photon (ignored).
  249 + * @return Model value.
248 250 *
249 251 * Evaluates the spatial part for an isotropic source model and set the
250 252 * parameter gradient. By definition, the value and gradient is independent
251 253 * from the sky direction. The value is 1, the parameter gradient is 0.
252 254 ***************************************************************************/
253   -double GModelSpatialDiffuseConst::eval_gradients(const GSkyDir& srcDir) const
  255 +double GModelSpatialDiffuseConst::eval_gradients(const GPhoton& photon) const
254 256 {
255 257 // Compute function value
256 258 double value = m_value.value();
... ... @@ -269,11 +271,16 @@ double GModelSpatialDiffuseConst::eval_gradients(const GSkyDir&amp; srcDir) const
269 271 /***********************************************************************//**
270 272 * @brief Return MC sky direction
271 273 *
272   - * @param[in] ran Random number generator.
  274 + * @param[in] energy Photon energy (ignored).
  275 + * @param[in] time Photon arrival time (ignored).
  276 + * @param[in,out] ran Random number generator.
  277 + * @return Sky direction.
273 278 *
274 279 * Returns an abitrary position on the celestial sphere.
275 280 ***************************************************************************/
276   -GSkyDir GModelSpatialDiffuseConst::mc(GRan& ran) const
  281 +GSkyDir GModelSpatialDiffuseConst::mc(const GEnergy& energy,
  282 + const GTime& time,
  283 + GRan& ran) const
277 284 {
278 285 // Simulate Right Ascension and Declination
279 286 double ra = twopi * ran.uniform();
... ...
src/model/GModelSpatialDiffuseCube.cpp
  View file @cc742ac
... ... @@ -268,14 +268,15 @@ GModelSpatialDiffuseCube* GModelSpatialDiffuseCube::clone(void) const
268 268 /***********************************************************************//**
269 269 * @brief Evaluate function
270 270 *
271   - * @param[in] srcDir True photon arrival direction.
  271 + * @param[in] photon Incident photon.
  272 + * @return Model value.
272 273 *
273 274 * @exception GException::feature_not_implemented
274 275 * Method not yet implemented
275 276 *
276 277 * @todo Implement method.
277 278 ***************************************************************************/
278   -double GModelSpatialDiffuseCube::eval(const GSkyDir& srcDir) const
  279 +double GModelSpatialDiffuseCube::eval(const GPhoton& photon) const
279 280 {
280 281 // Dump warning that method is not yet implemented
281 282 throw GException::feature_not_implemented(G_EVAL);
... ... @@ -288,11 +289,12 @@ double GModelSpatialDiffuseCube::eval(const GSkyDir&amp; srcDir) const
288 289 /***********************************************************************//**
289 290 * @brief Evaluate function and gradients
290 291 *
291   - * @param[in] srcDir True photon arrival direction.
  292 + * @param[in] photon Incident photon.
  293 + * @return Model value.
292 294 *
293 295 * @todo Implement method.
294 296 ***************************************************************************/
295   -double GModelSpatialDiffuseCube::eval_gradients(const GSkyDir& srcDir) const
  297 +double GModelSpatialDiffuseCube::eval_gradients(const GPhoton& photon) const
296 298 {
297 299 // Dump warning that method is not yet implemented
298 300 throw GException::feature_not_implemented(G_EVAL_GRADIENTS);
... ... @@ -308,14 +310,19 @@ double GModelSpatialDiffuseCube::eval_gradients(const GSkyDir&amp; srcDir) const
308 310 /***********************************************************************//**
309 311 * @brief Returns MC sky direction
310 312 *
311   - * @param[in] ran Random number generator.
  313 + * @param[in] energy Photon energy.
  314 + * @param[in] time Photon arrival time.
  315 + * @param[in,out] ran Random number generator.
  316 + * @return Sky direction.
312 317 *
313 318 * @exception GException::feature_not_implemented
314 319 * Method not yet implemented
315 320 *
316 321 * @todo Implement method.
317 322 ***************************************************************************/
318   -GSkyDir GModelSpatialDiffuseCube::mc(GRan& ran) const
  323 +GSkyDir GModelSpatialDiffuseCube::mc(const GEnergy& energy,
  324 + const GTime& time,
  325 + GRan& ran) const
319 326 {
320 327 // Allocate sky direction
321 328 GSkyDir dir;
... ...
src/model/GModelSpatialDiffuseMap.cpp
  View file @cc742ac
... ... @@ -224,17 +224,17 @@ GModelSpatialDiffuseMap* GModelSpatialDiffuseMap::clone(void) const
224 224 /***********************************************************************//**
225 225 * @brief Return intensity of skymap
226 226 *
227   - * @param[in] srcDir True photon arrival direction.
  227 + * @param[in] photon Incident photon.
228 228 * @return Sky map intensity.
229 229 *
230 230 * Returns the intensity of the skymap at the specified sky direction
231 231 * multiplied by the normalization factor. If the sky direction falls outside
232 232 * the skymap, an intensity of 0 is returned.
233 233 ***************************************************************************/
234   -double GModelSpatialDiffuseMap::eval(const GSkyDir& srcDir) const
  234 +double GModelSpatialDiffuseMap::eval(const GPhoton& photon) const
235 235 {
236 236 // Get skymap intensity
237   - double intensity = m_map(srcDir);
  237 + double intensity = m_map(photon.dir());
238 238  
239 239 // Return intensity times normalization factor
240 240 return (intensity * m_value.value());
... ... @@ -244,7 +244,7 @@ double GModelSpatialDiffuseMap::eval(const GSkyDir&amp; srcDir) const
244 244 /***********************************************************************//**
245 245 * @brief Return intensity of skymap and gradient
246 246 *
247   - * @param[in] srcDir True photon arrival direction.
  247 + * @param[in] photon Incident photon.
248 248 * @return Sky map intensity.
249 249 *
250 250 * Returns the intensity of the skymap at the specified sky direction
... ... @@ -252,10 +252,10 @@ double GModelSpatialDiffuseMap::eval(const GSkyDir&amp; srcDir) const
252 252 * with respect to the normalization factor. If the sky direction falls
253 253 * outside the skymap, an intensity of 0 is returned.
254 254 ***************************************************************************/
255   -double GModelSpatialDiffuseMap::eval_gradients(const GSkyDir& srcDir) const
  255 +double GModelSpatialDiffuseMap::eval_gradients(const GPhoton& photon) const
256 256 {
257 257 // Get skymap intensity
258   - double intensity = m_map(srcDir);
  258 + double intensity = m_map(photon.dir());
259 259  
260 260 // Compute partial derivatives of the parameter values
261 261 double g_value = (m_value.isfree()) ? intensity * m_value.scale() : 0.0;
... ... @@ -271,23 +271,23 @@ double GModelSpatialDiffuseMap::eval_gradients(const GSkyDir&amp; srcDir) const
271 271 /***********************************************************************//**
272 272 * @brief Returns MC sky direction
273 273 *
274   - * @param[in] ran Random number generator.
  274 + * @param[in] energy Photon energy (ignored).
  275 + * @param[in] time Photon arrival time (ignored).
  276 + * @param[in,out] ran Random number generator.
275 277 * @return Sky direction.
276 278 *
277   - * @exception GException::feature_not_implemented
278   - * Method not yet implemented
279   - *
280   - * This method returns a random sky direction according to the intensity
281   - * distribution of the model sky map. It makes use of a cache array that
282   - * contains the normalized cumulative flux values of the skymap. Using a
283   - * uniform random number, this cache array is scanned using a bi-section
284   - * method to determine the skymap pixel for which the position should be
285   - * returned. To avoid binning problems, the exact position within the pixel
286   - * is set by a uniform random number generator (neglecting thus pixel
287   - * distortions). The fractional skymap pixel is then converted into a sky
288   - * direction.
  279 + * Returns a random sky direction according to the intensity distribution of
  280 + * the model sky map. It makes use of a cache array that contains the
  281 + * normalized cumulative flux values of the skymap. Using a uniform random
  282 + * number, this cache array is scanned using a bi-section method to determine
  283 + * the skymap pixel for which the position should be returned. To avoid
  284 + * binning problems, the exact position within the pixel is set by a uniform
  285 + * random number generator (neglecting thus pixel distortions). The
  286 + * fractional skymap pixel is then converted into a sky direction.
289 287 ***************************************************************************/
290   -GSkyDir GModelSpatialDiffuseMap::mc(GRan& ran) const
  288 +GSkyDir GModelSpatialDiffuseMap::mc(const GEnergy& energy,
  289 + const GTime& time,
  290 + GRan& ran) const
291 291 {
292 292 // Allocate sky direction
293 293 GSkyDir dir;
... ...
src/model/GModelSpatialElliptical.cpp
  View file @cc742ac
... ... @@ -67,8 +67,8 @@ GModelSpatialElliptical::GModelSpatialElliptical(void) : GModelSpatial()
67 67 * @param[in] xml XML element.
68 68 *
69 69 * Constructs an elliptical spatial model component by extracting information
70   - * from an XML element. See read() for more information about the expected
71   - * structure of the XML element.
  70 + * from an XML element. See the read() method for more information about the
  71 + * expected structure of the XML element.
72 72 ***************************************************************************/
73 73 GModelSpatialElliptical::GModelSpatialElliptical(const GXmlElement& xml) :
74 74 GModelSpatial()
... ... @@ -160,19 +160,20 @@ GModelSpatialElliptical&amp; GModelSpatialElliptical::operator=(const GModelSpatialE
160 160 /***********************************************************************//**
161 161 * @brief Return model value
162 162 *
163   - * @param[in] srcDir True photon arrival direction.
  163 + * @param[in] photon Incident Photon.
164 164 *
165   - * Evaluates the elliptical spatial model for a given true photon arrival
166   - * direction.
  165 + * Evaluates the elliptical spatial model value for a specific incident
  166 + * @p photon.
167 167 ***************************************************************************/
168   -double GModelSpatialElliptical::eval(const GSkyDir& srcDir) const
  168 +double GModelSpatialElliptical::eval(const GPhoton& photon) const
169 169 {
170 170 // Compute distance from source and position angle (in radians)
171   - double theta = dir().dist(srcDir);
172   - double posang = dir().posang(srcDir);
  171 + const GSkyDir& srcDir = photon.dir();
  172 + double theta = dir().dist(srcDir);
  173 + double posang = dir().posang(srcDir);
173 174  
174 175 // Evaluate model
175   - double value = eval(theta, posang);
  176 + double value = eval(theta, posang, photon.energy(), photon.time());
176 177  
177 178 // Return result
178 179 return value;
... ... @@ -182,19 +183,21 @@ double GModelSpatialElliptical::eval(const GSkyDir&amp; srcDir) const
182 183 /***********************************************************************//**
183 184 * @brief Return model value and set analytical gradients
184 185 *
185   - * @param[in] srcDir True photon arrival direction.
  186 + * @param[in] photon Incident Photon.
186 187 *
187   - * Evaluates the elliptical spatial model for a given true photon arrival
188   - * direction.
  188 + * Evaluates the elliptical spatial model value and analytical model
  189 + * parameter gradients for a specific incident @p photon.
189 190 ***************************************************************************/
190   -double GModelSpatialElliptical::eval_gradients(const GSkyDir& srcDir) const
  191 +double GModelSpatialElliptical::eval_gradients(const GPhoton& photon) const
191 192 {
192 193 // Compute distance from source and position angle (in radians)
193   - double theta = dir().dist(srcDir);
194   - double posang = dir().posang(srcDir);
  194 + const GSkyDir& srcDir = photon.dir();
  195 + double theta = dir().dist(srcDir);
  196 + double posang = dir().posang(srcDir);
195 197  
196 198 // Evaluate model and set gradients
197   - double value = eval_gradients(theta, posang);
  199 + double value = eval_gradients(theta, posang, photon.energy(),
  200 + photon.time());
198 201  
199 202 // Return result
200 203 return value;
... ... @@ -211,7 +214,7 @@ double GModelSpatialElliptical::eval_gradients(const GSkyDir&amp; srcDir) const
211 214 * @exception GException::model_invalid_parnames
212 215 * Invalid model parameter names found in XML element.
213 216 *
214   - * Read the elliptical source location and position angle information from
  217 + * Reads the elliptical source location and position angle information from
215 218 * an XML element in the following format
216 219 *
217 220 * <spatialModel type="...">
... ... @@ -316,7 +319,7 @@ void GModelSpatialElliptical::read(const GXmlElement&amp; xml)
316 319 * @exception GException::model_invalid_parnames
317 320 * Invalid model parameter names found in XML element.
318 321 *
319   - * Write the elliptical source location and position angle information into
  322 + * Writes the elliptical source location and position angle information into
320 323 * an XML element in the following format
321 324 *
322 325 * <spatialModel type="...">
... ...
src/model/GModelSpatialEllipticalDisk.cpp
  View file @cc742ac
... ... @@ -234,6 +234,9 @@ GModelSpatialEllipticalDisk* GModelSpatialEllipticalDisk::clone(void) const
234 234 *
235 235 * @param[in] theta Angular distance from disk centre (radians).
236 236 * @param[in] posangle Position angle (counterclockwise from North) (radians).
  237 + * @param[in] energy Photon energy.
  238 + * @param[in] time Photon arrival time.
  239 + * @return Model value.
237 240 *
238 241 * Evaluates the spatial component for an elliptical disk source model. The
239 242 * disk source model is an elliptical function \f$f(\theta, \phi)\f$, where
... ... @@ -273,8 +276,10 @@ GModelSpatialEllipticalDisk* GModelSpatialEllipticalDisk::clone(void) const
273 276 *
274 277 * (see the update() method).
275 278 ***************************************************************************/
276   -double GModelSpatialEllipticalDisk::eval(const double& theta,
277   - const double& posangle) const
  279 +double GModelSpatialEllipticalDisk::eval(const double& theta,
  280 + const double& posangle,
  281 + const GEnergy& energy,
  282 + const GTime& time) const
278 283 {
279 284 // Initialise value
280 285 double value = 0.0;
... ... @@ -323,6 +328,9 @@ double GModelSpatialEllipticalDisk::eval(const double&amp; theta,
323 328 *
324 329 * @param[in] theta Angular distance from disk centre (radians).
325 330 * @param[in] posangle Position angle (counterclockwise from North) (radians).
  331 + * @param[in] energy Photon energy.
  332 + * @param[in] time Photon arrival time.
  333 + * @return Model value.
326 334 *
327 335 * Evaluates the function value. No gradient computation is implemented as
328 336 * Elliptical models will be convolved with the instrument response and thus
... ... @@ -330,44 +338,60 @@ double GModelSpatialEllipticalDisk::eval(const double&amp; theta,
330 338 *
331 339 * See the eval() method for more information.
332 340 ***************************************************************************/
333   -double GModelSpatialEllipticalDisk::eval_gradients(const double& theta,
334   - const double& posangle) const
  341 +double GModelSpatialEllipticalDisk::eval_gradients(const double& theta,
  342 + const double& posangle,
  343 + const GEnergy& energy,
  344 + const GTime& time) const
335 345 {
336 346 // Return value
337   - return (eval(theta, posangle));
  347 + return (eval(theta, posangle, energy, time));
338 348 }
339 349  
340 350  
341 351 /***********************************************************************//**
342 352 * @brief Returns MC sky direction
343 353 *
344   - * @param[in] ran Random number generator.
  354 + * @param[in] energy Photon energy.
  355 + * @param[in] time Photon arrival time.
  356 + * @param[in,out] ran Random number generator.
  357 + * @return Sky direction.
345 358 *
346 359 * Draws an arbitrary sky position from the 2D disk distribution.
347   - * @todo test function
  360 + *
  361 + * @todo Test function
348 362 ***************************************************************************/
349   -GSkyDir GModelSpatialEllipticalDisk::mc(GRan& ran) const
  363 +GSkyDir GModelSpatialEllipticalDisk::mc(const GEnergy& energy,
  364 + const GTime& time,
  365 + GRan& ran) const
350 366 {
351 367 // Update precomputation cache
352 368 update();
353 369  
354   - // Initialise a radial disk model with radius of the semimajor axis
355   - GModelSpatialRadialDisk disk = GModelSpatialRadialDisk(dir(), semimajor());
356   -
357   - // Initialise SkyDir
358   - GSkyDir dir = GSkyDir();
  370 + // Initialise photon
  371 + GPhoton photon;
  372 + photon.energy(energy);
  373 + photon.time(time);
359 374  
360 375 // Draw randomly from the radial disk
361 376 // and reject the value if its outside the ellipse
362 377 do {
363 378  
364   - // Set SkyDir to random position
365   - dir = disk.mc(ran);
  379 + // Simulate offset from photon arrival direction
  380 + double cosrad = std::cos(semimajor()*deg2rad);
  381 + double theta = std::acos(1.0 - ran.uniform() * (1.0 - cosrad)) * rad2deg;
  382 + double phi = 360.0 * ran.uniform();
  383 +
  384 + // Rotate sky direction by offset
  385 + GSkyDir sky_dir = dir();
  386 + sky_dir.rotate_deg(phi, theta);
  387 +
  388 + // Set photon sky direction
  389 + photon.dir(sky_dir);
366 390  
367   - } while(GModelSpatialElliptical::eval(dir) <= 0.0);
  391 + } while(GModelSpatialElliptical::eval(photon) <= 0.0);
368 392  
369   - // Return SkyDir
370   - return dir;
  393 + // Return photon direction
  394 + return (photon.dir());
371 395  
372 396 }
373 397  
... ...
src/model/GModelSpatialPointSource.cpp
  View file @cc742ac
... ... @@ -71,7 +71,7 @@ GModelSpatialPointSource::GModelSpatialPointSource(void) : GModelSpatial()
71 71  
72 72  
73 73 /***********************************************************************//**
74   - * @brief Constructor
  74 + * @brief Sky direction constructor
75 75 *
76 76 * @param[in] dir Sky direction.
77 77 *
... ... @@ -97,8 +97,8 @@ GModelSpatialPointSource::GModelSpatialPointSource(const GSkyDir&amp; dir) :
97 97 * @param[in] xml XML element.
98 98 *
99 99 * Construct a point source spatial model by extracting information from an
100   - * XML element. See GModelSpatialPointSource::read() for more information
101   - * about the expected structure of the XML element.
  100 + * XML element. See the read() method for more information about the expected
  101 + * structure of the XML element.
102 102 ***************************************************************************/
103 103 GModelSpatialPointSource::GModelSpatialPointSource(const GXmlElement& xml) :
104 104 GModelSpatial()
... ... @@ -189,7 +189,7 @@ GModelSpatialPointSource&amp; GModelSpatialPointSource::operator=(const GModelSpatia
189 189 ==========================================================================*/
190 190  
191 191 /***********************************************************************//**
192   - * @brief Clear instance
  192 + * @brief Clear point source model
193 193 ***************************************************************************/
194 194 void GModelSpatialPointSource::clear(void)
195 195 {
... ... @@ -207,12 +207,13 @@ void GModelSpatialPointSource::clear(void)
207 207  
208 208  
209 209 /***********************************************************************//**
210   - * @brief Clone instance
  210 + * @brief Clone point source model
211 211 *
212   - * @return Pointer to deep copy of point source spatial model.
  212 + * @return Pointer to deep copy of point source model.
213 213 ***************************************************************************/
214 214 GModelSpatialPointSource* GModelSpatialPointSource::clone(void) const
215 215 {
  216 + // Clone point source model
216 217 return new GModelSpatialPointSource(*this);
217 218 }
218 219  
... ... @@ -220,18 +221,18 @@ GModelSpatialPointSource* GModelSpatialPointSource::clone(void) const
220 221 /***********************************************************************//**
221 222 * @brief Evaluate function
222 223 *
223   - * @param[in] srcDir True photon arrival direction.
224   - * @return Model value (1 if srcDir corresponds to source direction, 0 else).
  224 + * @param[in] photon Incident photon.
  225 + * @return Model value.
225 226 *
226 227 * Evaluates the spatial part for a point source model. It implements a delta
227 228 * function with respect to the coordinates of the source. For numerical
228   - * reasons a certain tolerance is accepted (typically 0.1 arcsec, i.e. well
  229 + * reasons, a certain tolerance is accepted (typically 0.1 arcsec, i.e. well
229 230 * below the angular resolution of gamma-ray telescopes).
230 231 ***************************************************************************/
231   -double GModelSpatialPointSource::eval(const GSkyDir& srcDir) const
  232 +double GModelSpatialPointSource::eval(const GPhoton& photon) const
232 233 {
233 234 // Set value dependent on source distance
234   - double value = (srcDir.dist_deg(dir()) < tolerance) ? 1.0 : 0.0;
  235 + double value = (photon.dir().dist_deg(dir()) < tolerance) ? 1.0 : 0.0;
235 236  
236 237 // Return value
237 238 return value;
... ... @@ -241,8 +242,8 @@ double GModelSpatialPointSource::eval(const GSkyDir&amp; srcDir) const
241 242 /***********************************************************************//**
242 243 * @brief Evaluate function and gradients
243 244 *
244   - * @param[in] srcDir True photon arrival direction.
245   - * @return Model value (1 if srcDir corresponds to source direction, 0 else).
  245 + * @param[in] photon Incident photon.
  246 + * @return Model value.
246 247 *
247 248 * Evaluates the spatial part for a point source model and the gradient.
248 249 * It implements a delta function with respect to the coordinates of the
... ... @@ -252,28 +253,31 @@ double GModelSpatialPointSource::eval(const GSkyDir&amp; srcDir) const
252 253 *
253 254 * This method does not provide valid parameter gradients.
254 255 ***************************************************************************/
255   -double GModelSpatialPointSource::eval_gradients(const GSkyDir& srcDir) const
  256 +double GModelSpatialPointSource::eval_gradients(const GPhoton& photon) const
256 257 {
257   - // Set value dependent on source distance
258   - double value = (srcDir.dist_deg(dir()) < tolerance) ? 1.0 : 0.0;
259   -
260 258 // Return value
261   - return value;
  259 + return (eval(photon));
262 260 }
263 261  
264 262  
265 263 /***********************************************************************//**
266 264 * @brief Returns MC sky direction
267 265 *
268   - * @param[in] ran Random number generator.
  266 + * @param[in] energy Photon energy.
  267 + * @param[in] time Photon arrival time.
  268 + * @param[in,out] ran Random number generator.
269 269 * @return Sky direction.
270 270 *
271   - * This method always returns the directon of the point source.
  271 + * Draws an arbitrary sky direction for the point source model. As the point
  272 + * source is a point in the sky, the methods always returns the directon of
  273 + * the point source.
272 274 ***************************************************************************/
273   -GSkyDir GModelSpatialPointSource::mc(GRan& ran) const
  275 +GSkyDir GModelSpatialPointSource::mc(const GEnergy& energy,
  276 + const GTime& time,
  277 + GRan& ran) const
274 278 {
275   - // Return point source direction
276   - return dir();
  279 + // Return sky direction
  280 + return (dir());
277 281 }
278 282  
279 283  
... ... @@ -460,6 +464,8 @@ std::string GModelSpatialPointSource::print(void) const
460 464 * @brief Return position of point source
461 465 *
462 466 * @return Point source sky direction.
  467 + *
  468 + * Returns the sky direction of the point source.
463 469 ***************************************************************************/
464 470 GSkyDir GModelSpatialPointSource::dir(void) const
465 471 {
... ... @@ -476,6 +482,8 @@ GSkyDir GModelSpatialPointSource::dir(void) const
476 482  
477 483 /***********************************************************************//**
478 484 * @brief Set position of point source
  485 + *
  486 + * Sets the sky direction of the point source.
479 487 ***************************************************************************/
480 488 void GModelSpatialPointSource::dir(const GSkyDir& dir)
481 489 {
... ...
src/model/GModelSpatialRadial.cpp
  View file @cc742ac
... ... @@ -66,9 +66,9 @@ GModelSpatialRadial::GModelSpatialRadial(void) : GModelSpatial()
66 66 *
67 67 * @param[in] xml XML element.
68 68 *
69   - * Creates instance of a radial spatial model by extracting information
70   - * from an XML element. See GModelSpatialRadial::read() for more information about
71   - * the expected structure of the XML element.
  69 + * Constructs a radial spatial model by extracting information from an XML
  70 + * element. See the read() method for more information about the expected
  71 + * structure of the XML element.
72 72 ***************************************************************************/
73 73 GModelSpatialRadial::GModelSpatialRadial(const GXmlElement& xml) :
74 74 GModelSpatial()
... ... @@ -83,6 +83,7 @@ GModelSpatialRadial::GModelSpatialRadial(const GXmlElement&amp; xml) :
83 83 return;
84 84 }
85 85  
  86 +
86 87 /***********************************************************************//**
87 88 * @brief Copy constructor
88 89 *
... ... @@ -125,6 +126,7 @@ GModelSpatialRadial::~GModelSpatialRadial(void)
125 126 * @brief Assignment operator
126 127 *
127 128 * @param[in] model Radial spatial model.
  129 + * @return Radial spatial model.
128 130 ***************************************************************************/
129 131 GModelSpatialRadial& GModelSpatialRadial::operator=(const GModelSpatialRadial& model)
130 132 {
... ... @@ -159,18 +161,18 @@ GModelSpatialRadial&amp; GModelSpatialRadial::operator=(const GModelSpatialRadial&amp; m
159 161 /***********************************************************************//**
160 162 * @brief Return model value
161 163 *
162   - * @param[in] srcDir True photon arrival direction.
  164 + * @param[in] photon Incident Photon.
163 165 *
164   - * Evaluates the radial spatial model for a given true photon arrival
165   - * direction.
  166 + * Evaluates the radial spatial model value for a specific incident
  167 + * @p photon.
166 168 ***************************************************************************/
167   -double GModelSpatialRadial::eval(const GSkyDir& srcDir) const
  169 +double GModelSpatialRadial::eval(const GPhoton& photon) const
168 170 {
169 171 // Compute distance from source (in radians)
170   - double theta = srcDir.dist(dir());
  172 + double theta = photon.dir().dist(dir());
171 173  
172 174 // Evaluate model
173   - double value = eval(theta);
  175 + double value = eval(theta, photon.energy(), photon.time());
174 176  
175 177 // Return result
176 178 return value;
... ... @@ -180,18 +182,18 @@ double GModelSpatialRadial::eval(const GSkyDir&amp; srcDir) const
180 182 /***********************************************************************//**
181 183 * @brief Return model value and set analytical gradients
182 184 *
183   - * @param[in] srcDir True photon arrival direction.
  185 + * @param[in] photon Incident Photon.
184 186 *
185   - * Evaluates the radial spatial model for a given true photon arrival
186   - * direction.
  187 + * Evaluates the radial spatial model value and analytical model parameter
  188 + * gradients for a specific incident @p photon.
187 189 ***************************************************************************/
188   -double GModelSpatialRadial::eval_gradients(const GSkyDir& srcDir) const
  190 +double GModelSpatialRadial::eval_gradients(const GPhoton& photon) const
189 191 {
190 192 // Compute distance from source (in radians)
191   - double theta = srcDir.dist(dir());
  193 + double theta = photon.dir().dist(dir());
192 194  
193 195 // Evaluate model and set gradients
194   - double value = eval_gradients(theta);
  196 + double value = eval_gradients(theta, photon.energy(), photon.time());
195 197  
196 198 // Return result
197 199 return value;
... ... @@ -208,9 +210,23 @@ double GModelSpatialRadial::eval_gradients(const GSkyDir&amp; srcDir) const
208 210 * @exception GException::model_invalid_parnames
209 211 * Invalid model parameter names found in XML element.
210 212 *
211   - * Read the radial source location information from an XML element. The XML
212   - * element is required to have at least 2 parameters.
213   - * The location is named either "RA" and "DEC" or "GLON" and "GLAT".
  213 + * Reads the radial source location and position angle information from an
  214 + * XML element in the following format
  215 + *
  216 + * <spatialModel type="...">
  217 + * <parameter name="RA" scale="1" value="83.63" min="-360" max="360" free="1"/>
  218 + * <parameter name="DEC" scale="1" value="22.01" min="-90" max="90" free="1"/>
  219 + * ...
  220 + * </spatialModel>
  221 + *
  222 + * or
  223 + *
  224 + * <spatialModel type="...">
  225 + * <parameter name="GLON" scale="1" value="83.63" min="-360" max="360" free="1"/>
  226 + * <parameter name="GLAT" scale="1" value="22.01" min="-90" max="90" free="1"/>
  227 + * ...
  228 + * </spatialModel>
  229 + *
214 230 ***************************************************************************/
215 231 void GModelSpatialRadial::read(const GXmlElement& xml)
216 232 {
... ... @@ -291,8 +307,14 @@ void GModelSpatialRadial::read(const GXmlElement&amp; xml)
291 307 * @exception GException::model_invalid_parnames
292 308 * Invalid model parameter names found in XML element.
293 309 *
294   - * Write the radial source location information into an XML element. The XML
295   - * element will have at least 2 parameter leafs named "RA" and "DEC".
  310 + * Writes the radial source location and position angle information into an
  311 + * XML element in the following format
  312 + *
  313 + * <spatialModel type="...">
  314 + * <parameter name="RA" scale="1" value="83.63" min="-360" max="360" free="1"/>
  315 + * <parameter name="DEC" scale="1" value="22.01" min="-90" max="90" free="1"/>
  316 + * ...
  317 + * </spatialModel>
296 318 *
297 319 * @todo The case that an existing spatial XML element with "GLON" and "GLAT"
298 320 * as coordinates is not supported.
... ...
src/model/GModelSpatialRadialDisk.cpp
  View file @cc742ac
... ... @@ -97,8 +97,8 @@ GModelSpatialRadialDisk::GModelSpatialRadialDisk(const GSkyDir&amp; dir,
97 97 * @param[in] xml XML element.
98 98 *
99 99 * Creates instance of radial disk model by extracting information from an
100   - * XML element. See GModelSpatialRadialDisk::read() for more information about the
101   - * expected structure of the XML element.
  100 + * XML element. See the read() method for more information about the expected
  101 + * structure of the XML element.
102 102 ***************************************************************************/
103 103 GModelSpatialRadialDisk::GModelSpatialRadialDisk(const GXmlElement& xml) :
104 104 GModelSpatialRadial()
... ... @@ -189,7 +189,7 @@ GModelSpatialRadialDisk&amp; GModelSpatialRadialDisk::operator=(const GModelSpatialR
189 189 ==========================================================================*/
190 190  
191 191 /***********************************************************************//**
192   - * @brief Clear instance
  192 + * @brief Clear radial disk model
193 193 ***************************************************************************/
194 194 void GModelSpatialRadialDisk::clear(void)
195 195 {
... ... @@ -209,10 +209,11 @@ void GModelSpatialRadialDisk::clear(void)
209 209  
210 210  
211 211 /***********************************************************************//**
212   - * @brief Clone instance
  212 + * @brief Clone radial disk model
213 213 ***************************************************************************/
214 214 GModelSpatialRadialDisk* GModelSpatialRadialDisk::clone(void) const
215 215 {
  216 + // Clone radial disk model
216 217 return new GModelSpatialRadialDisk(*this);
217 218 }
218 219  
... ... @@ -221,21 +222,24 @@ GModelSpatialRadialDisk* GModelSpatialRadialDisk::clone(void) const
221 222 * @brief Evaluate function (in units of sr^-1)
222 223 *
223 224 * @param[in] theta Angular distance from disk centre (radians).
  225 + * @param[in] energy Photon energy.
  226 + * @param[in] time Photon arrival time.
  227 + * @return Model value.
224 228 *
225 229 * Evaluates the spatial part for a disk source model. The disk source
226 230 * model is a radial function \f$f(\theta)\f$, where \f$\theta\f$ is the
227 231 * angular separation between disk centre and the actual location:
228 232 * \f[
229   - * f(\theta) = \left \{
230   - * \begin{array}{l l}
231   - * \displaystyle
232   - * {\tt m\_norm}
233   - * & \mbox{if $\theta \le $ radius} \\
234   - * \\
235   - * \displaystyle
236   - * 0 & \mbox{if $\theta > $ radius}
237   - * \end{array}
238   - * \right .
  233 + * f(\theta) = \left \{
  234 + * \begin{array}{l l}
  235 + * \displaystyle
  236 + * {\tt m\_norm}
  237 + * & \mbox{if $\theta \le $ radius} \\
  238 + * \\
  239 + * \displaystyle
  240 + * 0 & \mbox{if $\theta > $ radius}
  241 + * \end{array}
  242 + * \right .
239 243 * \f]
240 244 *
241 245 * where
... ... @@ -244,7 +248,9 @@ GModelSpatialRadialDisk* GModelSpatialRadialDisk::clone(void) const
244 248 *
245 249 * is a normalization constant (see the update() method).
246 250 ***************************************************************************/
247   -double GModelSpatialRadialDisk::eval(const double& theta) const
  251 +double GModelSpatialRadialDisk::eval(const double& theta,
  252 + const GEnergy& energy,
  253 + const GTime& time) const
248 254 {
249 255 // Update precomputation cache
250 256 update();
... ... @@ -273,6 +279,9 @@ double GModelSpatialRadialDisk::eval(const double&amp; theta) const
273 279 * @brief Evaluate function and gradients (in units of sr^-1)
274 280 *
275 281 * @param[in] theta Angular distance from disk centre (radians).
  282 + * @param[in] energy Photon energy.
  283 + * @param[in] time Photon arrival time.
  284 + * @return Model value.
276 285 *
277 286 * Evaluates the function value. No gradient computation is implemented as
278 287 * radial models will be convolved with the instrument response and thus
... ... @@ -280,21 +289,29 @@ double GModelSpatialRadialDisk::eval(const double&amp; theta) const
280 289 *
281 290 * See the eval() method for more information.
282 291 ***************************************************************************/
283   -double GModelSpatialRadialDisk::eval_gradients(const double& theta) const
  292 +double GModelSpatialRadialDisk::eval_gradients(const double& theta,
  293 + const GEnergy& energy,
  294 + const GTime& time) const
284 295 {
285 296 // Return value
286   - return (eval(theta));
  297 + return (eval(theta, energy, time));
287 298 }
288 299  
289 300  
290 301 /***********************************************************************//**
291   - * @brief Returns MC sky direction
  302 + * @brief Return MC sky direction
292 303 *
293   - * @param[in] ran Random number generator.
  304 + * @param[in] energy Photon energy.
  305 + * @param[in] time Photon arrival time.
  306 + * @param[in,out] ran Random number generator.
  307 + * @return Sky direction.
294 308 *
295   - * Draws an arbitrary sky position from the 2D disk distribution.
  309 + * Draws an arbitrary sky direction from the 2D disk distribution as function
  310 + * of the photon @p energy and arrival @p time.
296 311 ***************************************************************************/
297   -GSkyDir GModelSpatialRadialDisk::mc(GRan& ran) const
  312 +GSkyDir GModelSpatialRadialDisk::mc(const GEnergy& energy,
  313 + const GTime& time,
  314 + GRan& ran) const
298 315 {
299 316 // Simulate offset from photon arrival direction
300 317 double cosrad = std::cos(radius()*deg2rad);
... ...
src/model/GModelSpatialRadialGauss.cpp
  View file @cc742ac
... ... @@ -224,6 +224,9 @@ GModelSpatialRadialGauss* GModelSpatialRadialGauss::clone(void) const
224 224 * @brief Evaluate function
225 225 *
226 226 * @param[in] theta Angular distance from Gaussian centre (radians).
  227 + * @param[in] energy Photon energy.
  228 + * @param[in] time Photon arrival time.
  229 + * @return Model value.
227 230 *
228 231 * Evaluates the spatial part for a Gaussian source model. The Gaussian
229 232 * source model is defined as
... ... @@ -236,7 +239,9 @@ GModelSpatialRadialGauss* GModelSpatialRadialGauss::clone(void) const
236 239 * @todo The Gaussian function is only correct in the small angle
237 240 * approximation.
238 241 ***************************************************************************/
239   -double GModelSpatialRadialGauss::eval(const double& theta) const
  242 +double GModelSpatialRadialGauss::eval(const double& theta,
  243 + const GEnergy& energy,
  244 + const GTime& time) const
240 245 {
241 246 // Compute value
242 247 double sigma_rad = sigma() * deg2rad;
... ... @@ -266,28 +271,38 @@ double GModelSpatialRadialGauss::eval(const double&amp; theta) const
266 271 * @brief Evaluate function and gradients
267 272 *
268 273 * @param[in] theta Angular distance from Gaussian centre (radians).
  274 + * @param[in] energy Photon energy.
  275 + * @param[in] time Photon arrival time.
  276 + * @return Model value.
269 277 *
270   - * This method simply calls GModelRadialShell::eval() as no analytical
271   - * gradients will be computed. See GModelRadialShell::eval() for details
272   - * about the implemented method.
  278 + * This method simply calls the eval() method as no analytical gradients
  279 + * will be computed. See the eval() method for more details.
273 280 ***************************************************************************/
274   -double GModelSpatialRadialGauss::eval_gradients(const double& theta) const
  281 +double GModelSpatialRadialGauss::eval_gradients(const double& theta,
  282 + const GEnergy& energy,
  283 + const GTime& time) const
275 284 {
276 285 // Return value
277   - return (eval(theta));
  286 + return (eval(theta, energy, time));
278 287 }
279 288  
280 289  
281 290 /***********************************************************************//**
282 291 * @brief Returns MC sky direction
283 292 *
284   - * @param[in] ran Random number generator.
  293 + * @param[in] energy Photon energy.
  294 + * @param[in] time Photon arrival time.
  295 + * @param[in,out] ran Random number generator.
  296 + * @return Sky direction.
285 297 *
286   - * Draws an arbitrary sky position from the 2D Gaussian distribution.
  298 + * Draws an arbitrary sky direction from the 2D Gaussian distribution as
  299 + * function of the photon @p energy and arrival @p time.
287 300 *
288 301 * @todo This method is only valid in the small angle approximation.
289 302 ***************************************************************************/
290   -GSkyDir GModelSpatialRadialGauss::mc(GRan& ran) const
  303 +GSkyDir GModelSpatialRadialGauss::mc(const GEnergy& energy,
  304 + const GTime& time,
  305 + GRan& ran) const
291 306 {
292 307 // Simulate offset from photon arrival direction
293 308 double theta = sigma() * ran.chisq2();
... ...
src/model/GModelSpatialRadialShell.cpp
  View file @cc742ac
... ... @@ -233,6 +233,9 @@ GModelSpatialRadialShell* GModelSpatialRadialShell::clone(void) const
233 233 * @brief Evaluate function
234 234 *
235 235 * @param[in] theta Angular distance from shell centre (radians).
  236 + * @param[in] energy Photon energy.
  237 + * @param[in] time Photon arrival time.
  238 + * @return Model value.
236 239 *
237 240 * Evaluates the spatial part for a shell source model. The shell source
238 241 * model is a radial function \f$f(\theta)\f$, where \f$\theta\f$ is the
... ... @@ -281,7 +284,9 @@ GModelSpatialRadialShell* GModelSpatialRadialShell::clone(void) const
281 284 * Here, \f$\theta_{\rm in}\f$ and \f$\theta_{\rm out}\f$ are the shell inner
282 285 * and outer radius.
283 286 ***************************************************************************/
284   -double GModelSpatialRadialShell::eval(const double& theta) const
  287 +double GModelSpatialRadialShell::eval(const double& theta,
  288 + const GEnergy& energy,
  289 + const GTime& time) const
285 290 {
286 291 // Update precomputation cache
287 292 update();
... ... @@ -333,26 +338,35 @@ double GModelSpatialRadialShell::eval(const double&amp; theta) const
333 338 * @brief Evaluate function and gradients
334 339 *
335 340 * @param[in] theta Angular distance from shell centre (radians).
  341 + * @param[in] energy Photon energy.
  342 + * @param[in] time Photon arrival time.
  343 + * @return Model value.
336 344 *
337   - * This method simply calls GModelSpatialRadialShell::eval() as no analytical
338   - * gradients will be computed. See GModelSpatialRadialShell::eval() for details
339   - * about the implemented method.
  345 + * This method simply calls the eval() method as no analytical gradients will
  346 + * be computed. See the eval() method for details.
340 347 ***************************************************************************/
341   -double GModelSpatialRadialShell::eval_gradients(const double& theta) const
  348 +double GModelSpatialRadialShell::eval_gradients(const double& theta,
  349 + const GEnergy& energy,
  350 + const GTime& time) const
342 351 {
343 352 // Return value
344   - return (eval(theta));
  353 + return (eval(theta, energy, time));
345 354 }
346 355  
347 356  
348 357 /***********************************************************************//**
349 358 * @brief Returns MC sky direction
350 359 *
351   - * @param[in] ran Random number generator.
  360 + * @param[in] energy Photon energy.
  361 + * @param[in] time Photon arrival time.
  362 + * @param[in,out] ran Random number generator.
  363 + * @return Sky direction.
352 364 *
353 365 * Draws an arbitrary sky position from the 2D shell distribution.
354 366 ***************************************************************************/
355   -GSkyDir GModelSpatialRadialShell::mc(GRan& ran) const
  367 +GSkyDir GModelSpatialRadialShell::mc(const GEnergy& energy,
  368 + const GTime& time,
  369 + GRan& ran) const
356 370 {
357 371 // Update precomputation cache
358 372 update();
... ... @@ -370,7 +384,7 @@ GSkyDir GModelSpatialRadialShell::mc(GRan&amp; ran) const
370 384 double theta = 0.0;
371 385 do {
372 386 theta = ran.uniform() * theta_max;
373   - value = eval(theta) * std::sin(theta);
  387 + value = eval(theta, energy, time) * std::sin(theta);
374 388 u = ran.uniform() * u_max;
375 389 #if defined(G_DEBUG_MC)
376 390 n_samples++;
... ...