Simon Bonnefoy / gammalib

Commit c135e7ceca61f203c46c583ce1d1bc3e31ad7ffa

Authored by Jürgen Knödlseder
1 parent 25342e05

Remove eval_gradients() methods and remplace by eval() method with optional grad… 

…ients argument (#1847)

The eval_gradients() methods duplicate much of the code of the eval() methods. To increase the maintainability of the code, the eval_gradients() methods have been removed, and the computation of partial derivatives is now done by the eval() methods if the optional gradients argument is set to true.
Inline Side-by-side
Showing 127 changed files with 1079 additions and 2351 deletions

Too many changes.
Reload with full diff Plain diff Email patch

To preserve performance only 100 of 127 files are displayed.

ChangeLog
  View file @c135e7c
... ... @@ -4,6 +4,7 @@
4 4 =========================
5 5  
6 6 Check CTA response table dimension for consistency (#1848)
  7 + Remove eval_gradients() methods and remplace by eval() method with optional gradients argument (#1847)
7 8 Add GTestSuite::test_value() method with automatic precision computation
8 9 Add GSkyMap::is_empty() method (#1828)
9 10 Add GApplication access by index operator (#1838)
... ...
... ... @@ -36,14 +36,20 @@ The following methods have been added:
36 36 - GSkyMap::is_empty()
37 37  
38 38 The following methods have been removed:
39   -- ...
  39 +- <all classes>::eval_gradients()
  40 + The eval_gradients() methods duplicate much of the code of the eval()
  41 + methods. To increase the maintainability of the code, the eval_gradients()
  42 + methods have been removed, and the computation of partial derivatives is
  43 + now done by the eval() methods if the optional gradients argument is set
  44 + to true.
40 45  
41 46 The following methods have been renamed:
42 47 - GModelSpectralPlawPhotonFlux::flux(), was GModelSpectralPlaw2::integral()
43 48 - GLog::buffer_size(), was GLog::max_size()
44 49  
45   -The parameters for the following methods have been changed:
  50 +The arguments for the following methods have been changed:
46 51 - GApplication::log_parameters() now has a mandatory GChatter argument
  52 +- <all classes>::eval() has now an optional gradients parameter (see above)
47 53  
48 54 The return value of the following methods has been changed:
49 55 - GSkyMap::write() now returns a pointer to GFitsHDU instead of void
... ...
doc/source/users/user_manual/modules/model.rst
  View file @c135e7c
... ... @@ -35,7 +35,7 @@ methods.
35 35  
36 36 For model evaluation, the :doxy:``GObservation::model`` method loops over all
37 37 models in the :doxy:`GModels` container and calls the
38   -:doxy:`GModel::eval_gradients` method of each model.
  38 +:doxy:`GModel::eval` method of each model.
39 39 For :doxy:`GModelSky` models this method invokes the :doxy:`GResponse::convolve`
40 40 method that calls among others the :doxy:`GResponse::irf` method.
41 41 That method needs to be implement for each instruments.
... ... @@ -46,7 +46,7 @@ The integrate over observed time is for the moment a simple multiplication
46 46 with ontime.
47 47 The integral over observed energy is performed by the
48 48 ``GObservation::npred_spec`` method that calls the
49   -:doxy:`GModel::eval_gradients` method.
  49 +:doxy:`GModel::eval` method.
50 50 For :doxy:`GModelSky` models this method calls the :doxy:`GResponse::nroi`
51 51 method that needs to be implement for each instruments.
52 52  
... ...
doc/source/users/user_manual/modules/obs.rst
  View file @c135e7c
... ... @@ -233,9 +233,9 @@ takes the :doxy:`GSource` class instead of the :doxy:`GPhoton` class as argument
233 233 :doxy:`GSource` differs from :doxy:`GPhoton` in that the photon arrival direction
234 234 :math:`p` is replaced by the spatial component :doxy:`GModelSpatial` of a source
235 235 model.
236   -Equation :eq:`pirf` is used by the :doxy:`GModelSky::eval` and
237   -:doxy:`GModelSky::eval_gradients` methods for computation of the instrument
238   -response to a source model (see :ref:`fig_calltree_model`).
  236 +Equation :eq:`pirf` is used by the :doxy:`GModelSky::eval` method for
  237 +computation of the instrument response to a source model
  238 +(see :ref:`fig_calltree_model`).
239 239  
240 240 A maximum likelihood analysis of the data generally needs the computation of the
241 241 predicted number of events within the selection region for each source model.
... ...
include/GModel.hpp
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GModel.hpp - Abstract virtual model base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -46,11 +46,10 @@ class GObservation;
46 46 *
47 47 * @brief Abstract model class
48 48 *
49   - * This class implements a parametric model. The class has two methods for
50   - * model evaluation: eval() and eval_gradients().
51   - * The eval() method evaluates the model for a given event and observation.
52   - * In addition, eval_gradients() also sets the parameter gradients of the
53   - * model.
  49 + * This class implements a parametric model. The eval() method evaluates the
  50 + * model for a given event and observation. If the gradients parameter is set
  51 + * to true, the eval() method also computes analytical parameter gradients if
  52 + * they exist.
54 53 *
55 54 * A model has the following attributes:
56 55 * - @p name
... ... @@ -118,9 +117,8 @@ public:
118 117 virtual std::string type(void) const = 0;
119 118 virtual bool is_constant(void) const = 0;
120 119 virtual double eval(const GEvent& event,
121   - const GObservation& obs) const = 0;
122   - virtual double eval_gradients(const GEvent& event,
123   - const GObservation& obs) const = 0;
  120 + const GObservation& obs,
  121 + const bool& gradients = false) const = 0;
124 122 virtual double npred(const GEnergy& obsEng, const GTime& obsTime,
125 123 const GObservation& obs) const = 0;
126 124 virtual void read(const GXmlElement& xml) = 0;
... ...
include/GModelData.hpp
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GModelData.hpp - Abstract virtual data model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2013 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -69,9 +69,8 @@ public:
69 69 virtual std::string type(void) const = 0;
70 70 virtual bool is_constant(void) const = 0;
71 71 virtual double eval(const GEvent& event,
72   - const GObservation& obs) const = 0;
73   - virtual double eval_gradients(const GEvent& event,
74   - const GObservation& obs) const = 0;
  72 + const GObservation& obs,
  73 + const bool& gradients = false) const = 0;
75 74 virtual double npred(const GEnergy& obsEng, const GTime& obsTime,
76 75 const GObservation& obs) const = 0;
77 76 virtual GEvents* mc(const GObservation& obs, GRan& ran) const = 0;
... ...
include/GModelSky.hpp
  View file @c135e7c
... ... @@ -92,10 +92,8 @@ class GObservation;
92 92 *
93 93 * The class has two methods for model evaluation that evaluate the model
94 94 * for a specific event, given an observation. The eval() method returns
95   - * the model value, the eval_gradients() returns the model value and sets
96   - * the analytical gradients for all model parameters.
97   - * Note that the eval() and eval_gradients() methods call protected
98   - * methods that handle time dispersion, energy dispersion and the point
  95 + * the model value. Note that the eval() method calls protected methods that
  96 + * handle time dispersion, energy dispersion and the point
99 97 * spread function (spatial dispersion). Dispersion is handled by
100 98 * integrating over the relevant intervals of the source properties.
101 99 * Integration of the model is done by three nested integrals.
... ... @@ -149,9 +147,8 @@ public:
149 147 virtual std::string type(void) const;
150 148 virtual bool is_constant(void) const;
151 149 virtual double eval(const GEvent& event,
152   - const GObservation& obs) const;
153   - virtual double eval_gradients(const GEvent& event,
154   - const GObservation& obs) const;
  150 + const GObservation& obs,
  151 + const bool& gradients = false) const;
155 152 virtual double npred(const GEnergy& obsEng,
156 153 const GTime& obsTime,
157 154 const GObservation& obs) const;
... ...
include/GModelSpatial.hpp
  View file @c135e7c
... ... @@ -87,8 +87,8 @@ public:
87 87 virtual std::string classname(void) const = 0;
88 88 virtual std::string type(void) const = 0;
89 89 virtual GClassCode code(void) const = 0;
90   - virtual double eval(const GPhoton& photon) const = 0;
91   - virtual double eval_gradients(const GPhoton& photon) const = 0;
  90 + virtual double eval(const GPhoton& photon,
  91 + const bool& gradients = false) const = 0;
92 92 virtual GSkyDir mc(const GEnergy& energy, const GTime& time,
93 93 GRan& ran) const = 0;
94 94 virtual double mc_norm(const GSkyDir& dir,
... ...
include/GModelSpatialDiffuse.hpp
  View file @c135e7c
... ... @@ -62,8 +62,8 @@ public:
62 62 virtual GModelSpatialDiffuse* clone(void) const = 0;
63 63 virtual std::string classname(void) const = 0;
64 64 virtual std::string type(void) const = 0;
65   - virtual double eval(const GPhoton& photon) const = 0;
66   - virtual double eval_gradients(const GPhoton& photon) const = 0;
  65 + virtual double eval(const GPhoton& photon,
  66 + const bool& gradients = false) const = 0;
67 67 virtual GSkyDir mc(const GEnergy& energy, const GTime& time,
68 68 GRan& ran) const = 0;
69 69 virtual double mc_norm(const GSkyDir& dir,
... ...
include/GModelSpatialDiffuseConst.hpp
  View file @c135e7c
... ... @@ -63,8 +63,8 @@ public:
63 63 virtual GModelSpatialDiffuseConst* clone(void) const;
64 64 virtual std::string classname(void) const;
65 65 virtual std::string type(void) const;
66   - virtual double eval(const GPhoton& photon) const;
67   - virtual double eval_gradients(const GPhoton& photon) const;
  66 + virtual double eval(const GPhoton& photon,
  67 + const bool& gradients = false) const;
68 68 virtual GSkyDir mc(const GEnergy& energy,
69 69 const GTime& time,
70 70 GRan& ran) const;
... ...
include/GModelSpatialDiffuseCube.hpp
  View file @c135e7c
... ... @@ -78,8 +78,8 @@ public:
78 78 virtual GModelSpatialDiffuseCube* clone(void) const;
79 79 virtual std::string classname(void) const;
80 80 virtual std::string type(void) const;
81   - virtual double eval(const GPhoton& photon) const;
82   - virtual double eval_gradients(const GPhoton& photon) const;
  81 + virtual double eval(const GPhoton& photon,
  82 + const bool& gradients = false) const;
83 83 virtual GSkyDir mc(const GEnergy& energy,
84 84 const GTime& time,
85 85 GRan& ran) const;
... ...
include/GModelSpatialDiffuseMap.hpp
  View file @c135e7c
... ... @@ -70,8 +70,8 @@ public:
70 70 virtual GModelSpatialDiffuseMap* clone(void) const;
71 71 virtual std::string classname(void) const;
72 72 virtual std::string type(void) const;
73   - virtual double eval(const GPhoton& photon) const;
74   - virtual double eval_gradients(const GPhoton& photon) const;
  73 + virtual double eval(const GPhoton& photon,
  74 + const bool& gradients = false) const;
75 75 virtual GSkyDir mc(const GEnergy& energy,
76 76 const GTime& time,
77 77 GRan& ran) const;
... ...
include/GModelSpatialElliptical.hpp
  View file @c135e7c
... ... @@ -69,14 +69,11 @@ public:
69 69 virtual double eval(const double& theta,
70 70 const double& posangle,
71 71 const GEnergy& energy,
72   - const GTime& time) const = 0;
73   - virtual double eval_gradients(const double& theta,
74   - const double& posangle,
75   - const GEnergy& energy,
76   - const GTime& time) const = 0;
  72 + const GTime& time,
  73 + const bool& gradients = false) const = 0;
77 74 virtual GSkyDir mc(const GEnergy& energy,
78   - const GTime& time,
79   - GRan& ran) const = 0;
  75 + const GTime& time,
  76 + GRan& ran) const = 0;
80 77 virtual bool contains(const GSkyDir& dir,
81 78 const double& margin = 0.0) const = 0;
82 79 virtual double theta_max(void) const = 0;
... ... @@ -84,8 +81,8 @@ public:
84 81  
85 82 // Implemented virtual base class methods
86 83 virtual GClassCode code(void) const;
87   - virtual double eval(const GPhoton& photon) const;
88   - virtual double eval_gradients(const GPhoton& photon) const;
  84 + virtual double eval(const GPhoton& photon,
  85 + const bool& gradients = false) const;
89 86 virtual double mc_norm(const GSkyDir& dir, const double& radius) const;
90 87 virtual void read(const GXmlElement& xml);
91 88 virtual void write(GXmlElement& xml) const;
... ...
include/GModelSpatialEllipticalDisk.hpp
  View file @c135e7c
... ... @@ -68,11 +68,8 @@ public:
68 68 virtual double eval(const double& theta,
69 69 const double& posangle,
70 70 const GEnergy& energy,
71   - const GTime& time) const;
72   - virtual double eval_gradients(const double& theta,
73   - const double& posangle,
74   - const GEnergy& energy,
75   - const GTime& time) const;
  71 + const GTime& time,
  72 + const bool& gradients = false) const;
76 73 virtual GSkyDir mc(const GEnergy& energy,
77 74 const GTime& time,
78 75 GRan& ran) const;
... ...
include/GModelSpatialEllipticalGauss.hpp
  View file @c135e7c
... ... @@ -69,11 +69,8 @@ public:
69 69 virtual double eval(const double& theta,
70 70 const double& posangle,
71 71 const GEnergy& energy,
72   - const GTime& time) const;
73   - virtual double eval_gradients(const double& theta,
74   - const double& posangle,
75   - const GEnergy& energy,
76   - const GTime& time) const;
  72 + const GTime& time,
  73 + const bool& gradients = false) const;
77 74 virtual GSkyDir mc(const GEnergy& energy,
78 75 const GTime& time,
79 76 GRan& ran) const;
... ...
include/GModelSpatialPointSource.hpp
  View file @c135e7c
... ... @@ -65,8 +65,8 @@ public:
65 65 virtual std::string classname(void) const;
66 66 virtual std::string type(void) const;
67 67 virtual GClassCode code(void) const;
68   - virtual double eval(const GPhoton& photon) const;
69   - virtual double eval_gradients(const GPhoton& photon) const;
  68 + virtual double eval(const GPhoton& photon,
  69 + const bool& gradients = false) const;
70 70 virtual GSkyDir mc(const GEnergy& energy,
71 71 const GTime& time,
72 72 GRan& ran) const;
... ...
include/GModelSpatialRadial.hpp
  View file @c135e7c
... ... @@ -68,13 +68,11 @@ public:
68 68 virtual std::string type(void) const = 0;
69 69 virtual double eval(const double& theta,
70 70 const GEnergy& energy,
71   - const GTime& time) const = 0;
72   - virtual double eval_gradients(const double& theta,
73   - const GEnergy& energy,
74   - const GTime& time) const = 0;
  71 + const GTime& time,
  72 + const bool& gradients = false) const = 0;
75 73 virtual GSkyDir mc(const GEnergy& energy,
76   - const GTime& time,
77   - GRan& ran) const = 0;
  74 + const GTime& time,
  75 + GRan& ran) const = 0;
78 76 virtual bool contains(const GSkyDir& dir,
79 77 const double& margin = 0.0) const = 0;
80 78 virtual double theta_max(void) const = 0;
... ... @@ -82,8 +80,8 @@ public:
82 80  
83 81 // Implemented pure virtual base class methods
84 82 virtual GClassCode code(void) const;
85   - virtual double eval(const GPhoton& photon) const;
86   - virtual double eval_gradients(const GPhoton& photon) const;
  83 + virtual double eval(const GPhoton& photon,
  84 + const bool& gradients = false) const;
87 85 virtual double mc_norm(const GSkyDir& dir, const double& radius) const;
88 86 virtual void read(const GXmlElement& xml);
89 87 virtual void write(GXmlElement& xml) const;
... ...
include/GModelSpatialRadialDisk.hpp
  View file @c135e7c
... ... @@ -65,10 +65,8 @@ public:
65 65 virtual std::string type(void) const;
66 66 virtual double eval(const double& theta,
67 67 const GEnergy& energy,
68   - const GTime& time) const;
69   - virtual double eval_gradients(const double& theta,
70   - const GEnergy& energy,
71   - const GTime& time) const;
  68 + const GTime& time,
  69 + const bool& gradients = false) const;
72 70 virtual GSkyDir mc(const GEnergy& energy,
73 71 const GTime& time,
74 72 GRan& ran) const;
... ...
include/GModelSpatialRadialGauss.hpp
  View file @c135e7c
... ... @@ -66,10 +66,8 @@ public:
66 66 virtual std::string type(void) const;
67 67 virtual double eval(const double& theta,
68 68 const GEnergy& energy,
69   - const GTime& time) const;
70   - virtual double eval_gradients(const double& theta,
71   - const GEnergy& energy,
72   - const GTime& time) const;
  69 + const GTime& time,
  70 + const bool& gradients = false) const;
73 71 virtual GSkyDir mc(const GEnergy& energy,
74 72 const GTime& time,
75 73 GRan& ran) const;
... ...
include/GModelSpatialRadialShell.hpp
  View file @c135e7c
... ... @@ -77,10 +77,8 @@ public:
77 77 virtual std::string type(void) const;
78 78 virtual double eval(const double& theta,
79 79 const GEnergy& energy,
80   - const GTime& time) const;
81   - virtual double eval_gradients(const double& theta,
82   - const GEnergy& energy,
83   - const GTime& time) const;
  80 + const GTime& time,
  81 + const bool& gradients = false) const;
84 82 virtual GSkyDir mc(const GEnergy& energy,
85 83 const GTime& time,
86 84 GRan& ran) const;
... ...
include/GModelSpectral.hpp
  View file @c135e7c
... ... @@ -86,9 +86,8 @@ public:
86 86 virtual std::string classname(void) const = 0;
87 87 virtual std::string type(void) const = 0;
88 88 virtual double eval(const GEnergy& srcEng,
89   - const GTime& srcTime = GTime()) const = 0;
90   - virtual double eval_gradients(const GEnergy& srcEng,
91   - const GTime& srcTime = GTime()) = 0;
  89 + const GTime& srcTime = GTime(),
  90 + const bool& gradients = false) const = 0;
92 91 virtual double flux(const GEnergy& emin,
93 92 const GEnergy& emax) const = 0;
94 93 virtual double eflux(const GEnergy& emin,
... ...
include/GModelSpectralBrokenPlaw.hpp
  View file @c135e7c
... ... @@ -84,9 +84,8 @@ public:
84 84 virtual std::string classname(void) const;
85 85 virtual std::string type(void) const;
86 86 virtual double eval(const GEnergy& srcEng,
87   - const GTime& srcTime = GTime()) const;
88   - virtual double eval_gradients(const GEnergy& srcEng,
89   - const GTime& srcTime = GTime());
  87 + const GTime& srcTime = GTime(),
  88 + const bool& gradients = false) const;
90 89 virtual double flux(const GEnergy& emin,
91 90 const GEnergy& emax) const;
92 91 virtual double eflux(const GEnergy& emin,
... ...
include/GModelSpectralConst.hpp
  View file @c135e7c
... ... @@ -74,9 +74,8 @@ public:
74 74 virtual std::string classname(void) const;
75 75 virtual std::string type(void) const;
76 76 virtual double eval(const GEnergy& srcEng,
77   - const GTime& srcTime = GTime()) const;
78   - virtual double eval_gradients(const GEnergy& srcEng,
79   - const GTime& srcTime = GTime());
  77 + const GTime& srcTime = GTime(),
  78 + const bool& gradients = false) const;
80 79 virtual double flux(const GEnergy& emin,
81 80 const GEnergy& emax) const;
82 81 virtual double eflux(const GEnergy& emin,
... ...
include/GModelSpectralExpInvPlaw.hpp
  View file @c135e7c
... ... @@ -91,9 +91,8 @@ public:
91 91 virtual std::string classname(void) const;
92 92 virtual std::string type(void) const;
93 93 virtual double eval(const GEnergy& srcEng,
94   - const GTime& srcTime = GTime()) const;
95   - virtual double eval_gradients(const GEnergy& srcEng,
96   - const GTime& srcTime = GTime());
  94 + const GTime& srcTime = GTime(),
  95 + const bool& gradients = false) const;
97 96 virtual double flux(const GEnergy& emin,
98 97 const GEnergy& emax) const;
99 98 virtual double eflux(const GEnergy& emin,
... ...
include/GModelSpectralExpPlaw.hpp
  View file @c135e7c
... ... @@ -87,9 +87,8 @@ public:
87 87 virtual std::string classname(void) const;
88 88 virtual std::string type(void) const;
89 89 virtual double eval(const GEnergy& srcEng,
90   - const GTime& srcTime = GTime()) const;
91   - virtual double eval_gradients(const GEnergy& srcEng,
92   - const GTime& srcTime = GTime());
  90 + const GTime& srcTime = GTime(),
  91 + const bool& gradients = false) const;
93 92 virtual double flux(const GEnergy& emin,
94 93 const GEnergy& emax) const;
95 94 virtual double eflux(const GEnergy& emin,
... ...
include/GModelSpectralFunc.hpp
  View file @c135e7c
... ... @@ -77,9 +77,8 @@ public:
77 77 virtual std::string classname(void) const;
78 78 virtual std::string type(void) const;
79 79 virtual double eval(const GEnergy& srcEng,
80   - const GTime& srcTime = GTime()) const;
81   - virtual double eval_gradients(const GEnergy& srcEng,
82   - const GTime& srcTime = GTime());
  80 + const GTime& srcTime = GTime(),
  81 + const bool& gradients = false) const;
83 82 virtual double flux(const GEnergy& emin,
84 83 const GEnergy& emax) const;
85 84 virtual double eflux(const GEnergy& emin,
... ...
include/GModelSpectralGauss.hpp
  View file @c135e7c
... ... @@ -79,9 +79,8 @@ public:
79 79 virtual std::string classname(void) const;
80 80 virtual std::string type(void) const;
81 81 virtual double eval(const GEnergy& srcEng,
82   - const GTime& srcTime = GTime()) const;
83   - virtual double eval_gradients(const GEnergy& srcEng,
84   - const GTime& srcTime = GTime());
  82 + const GTime& srcTime = GTime(),
  83 + const bool& gradients = false) const;
85 84 virtual double flux(const GEnergy& emin,
86 85 const GEnergy& emax) const;
87 86 virtual double eflux(const GEnergy& emin,
... ...
include/GModelSpectralLogParabola.hpp
  View file @c135e7c
... ... @@ -88,9 +88,8 @@ public:
88 88 virtual std::string classname(void) const;
89 89 virtual std::string type(void) const;
90 90 virtual double eval(const GEnergy& srcEng,
91   - const GTime& srcTime = GTime()) const;
92   - virtual double eval_gradients(const GEnergy& srcEng,
93   - const GTime& srcTime = GTime());
  91 + const GTime& srcTime = GTime(),
  92 + const bool& gradients = false) const;
94 93 virtual double flux(const GEnergy& emin,
95 94 const GEnergy& emax) const;
96 95 virtual double eflux(const GEnergy& emin,
... ...
include/GModelSpectralNodes.hpp
  View file @c135e7c
... ... @@ -99,9 +99,8 @@ public:
99 99 virtual std::string classname(void) const;
100 100 virtual std::string type(void) const;
101 101 virtual double eval(const GEnergy& srcEng,
102   - const GTime& srcTime = GTime()) const;
103   - virtual double eval_gradients(const GEnergy& srcEng,
104   - const GTime& srcTime = GTime());
  102 + const GTime& srcTime = GTime(),
  103 + const bool& gradients = false) const;
105 104 virtual double flux(const GEnergy& emin,
106 105 const GEnergy& emax) const;
107 106 virtual double eflux(const GEnergy& emin,
... ...
include/GModelSpectralPlaw.hpp
  View file @c135e7c
... ... @@ -82,9 +82,8 @@ public:
82 82 virtual std::string classname(void) const;
83 83 virtual std::string type(void) const;
84 84 virtual double eval(const GEnergy& srcEng,
85   - const GTime& srcTime = GTime()) const;
86   - virtual double eval_gradients(const GEnergy& srcEng,
87   - const GTime& srcTime = GTime());
  85 + const GTime& srcTime = GTime(),
  86 + const bool& gradients = false) const;
88 87 virtual double flux(const GEnergy& emin,
89 88 const GEnergy& emax) const;
90 89 virtual double eflux(const GEnergy& emin,
... ...
include/GModelSpectralPlawEnergyFlux.hpp
  View file @c135e7c
... ... @@ -97,9 +97,8 @@ public:
97 97 virtual std::string classname(void) const;
98 98 virtual std::string type(void) const;
99 99 virtual double eval(const GEnergy& srcEng,
100   - const GTime& srcTime = GTime()) const;
101   - virtual double eval_gradients(const GEnergy& srcEng,
102   - const GTime& srcTime = GTime());
  100 + const GTime& srcTime = GTime(),
  101 + const bool& gradients = false) const;
103 102 virtual double flux(const GEnergy& emin,
104 103 const GEnergy& emax) const;
105 104 virtual double eflux(const GEnergy& emin,
... ...
include/GModelSpectralPlawPhotonFlux.hpp
  View file @c135e7c
... ... @@ -97,9 +97,8 @@ public:
97 97 virtual std::string classname(void) const;
98 98 virtual std::string type(void) const;
99 99 virtual double eval(const GEnergy& srcEng,
100   - const GTime& srcTime = GTime()) const;
101   - virtual double eval_gradients(const GEnergy& srcEng,
102   - const GTime& srcTime = GTime());
  100 + const GTime& srcTime = GTime(),
  101 + const bool& gradients = false) const;
103 102 virtual double flux(const GEnergy& emin,
104 103 const GEnergy& emax) const;
105 104 virtual double eflux(const GEnergy& emin,
... ...
include/GModelSpectralSuperExpPlaw.hpp
  View file @c135e7c
... ... @@ -92,9 +92,8 @@ public:
92 92 virtual std::string classname(void) const;
93 93 virtual std::string type(void) const;
94 94 virtual double eval(const GEnergy& srcEng,
95   - const GTime& srcTime = GTime()) const;
96   - virtual double eval_gradients(const GEnergy& srcEng,
97   - const GTime& srcTime = GTime());
  95 + const GTime& srcTime = GTime(),
  96 + const bool& gradients = false) const;
98 97 virtual double flux(const GEnergy& emin,
99 98 const GEnergy& emax) const;
100 99 virtual double eflux(const GEnergy& emin,
... ...
include/GModelTemporal.hpp
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GModelTemporal.hpp - Abstract temporal model base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -68,8 +68,8 @@ public:
68 68 virtual GModelTemporal* clone(void) const = 0;
69 69 virtual std::string classname(void) const = 0;
70 70 virtual std::string type(void) const = 0;
71   - virtual double eval(const GTime& srcTime) const = 0;
72   - virtual double eval_gradients(const GTime& srcTime) = 0;
  71 + virtual double eval(const GTime& srcTime,
  72 + const bool& gradients = false) const = 0;
73 73 virtual GTimes mc(const double& rate, const GTime& tmin,
74 74 const GTime& tmax, GRan& ran) const = 0;
75 75 virtual void read(const GXmlElement& xml) = 0;
... ...
include/GModelTemporalConst.hpp
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GModelTemporalConst.hpp - Temporal constant model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -70,8 +70,8 @@ public:
70 70 virtual GModelTemporalConst* clone(void) const;
71 71 virtual std::string classname(void) const;
72 72 virtual std::string type(void) const;
73   - virtual double eval(const GTime& srcTime) const;
74   - virtual double eval_gradients(const GTime& srcTime);
  73 + virtual double eval(const GTime& srcTime,
  74 + const bool& gradients = false) const;
75 75 virtual GTimes mc(const double& rate, const GTime& tmin,
76 76 const GTime& tmax, GRan& ran) const;
77 77 virtual void read(const GXmlElement& xml);
... ...
include/GModels.hpp
  View file @c135e7c
... ... @@ -132,13 +132,17 @@ class GFilename;
132 132 * models.read(xml); // Read models from GXml object
133 133 *
134 134 * The sum of all models in the container are evaluated for a given @p event
135   - * and @p observation using the eval() or eval_gradients() methods:
  135 + * and @p observation using the eval() method:
136 136 *
137 137 * double value = models.eval(event, observation);
138   - * double value = models.eval_gradients(event, observation);
139 138 *
140   - * The eval_gradients() method sets the parameter gradients for all free
141   - * model parameters that have an analytical parameter gradient.
  139 + * If the eval() method is called with the optional gradients parameter set
  140 + * to true, i.e.
  141 + *
  142 + * double value = models.eval(event, observation, true);
  143 + *
  144 + * the method computes also the parameter gradients for all free model
  145 + * parameters that have an analytical parameter gradient.
142 146 *
143 147 * The only member of GModels is a list of model pointers. The class handles
144 148 * the proper allocation and deallocation of the model memory.
... ... @@ -183,8 +187,9 @@ public:
183 187 void write(GXml& xml) const;
184 188 int npars(void) const;
185 189 GOptimizerPars pars(void) const;
186   - double eval(const GEvent& event, const GObservation& obs) const;
187   - double eval_gradients(const GEvent& event, const GObservation& obs) const;
  190 + double eval(const GEvent& event,
  191 + const GObservation& obs,
  192 + const bool& gradients = false) const;
188 193 std::string print(const GChatter& chatter = NORMAL) const;
189 194  
190 195 protected:
... ...
include/GOptimizerPar.hpp
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GOptimizerPar.hpp - Optimizer parameter class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2013-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2013-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -122,7 +122,7 @@ public:
122 122 const double& scale(void) const;
123 123 void factor_value(const double& value);
124 124 void factor_error(const double& error);
125   - void factor_gradient(const double& gradient);
  125 + void factor_gradient(const double& gradient) const;
126 126 void factor_min(const double& min);
127 127 void factor_max(const double& max);
128 128 void factor_range(const double& min, const double& max);
... ... @@ -162,18 +162,18 @@ protected:
162 162 void free_members(void);
163 163  
164 164 // Proteced data members
165   - std::string m_name; //!< Parameter name
166   - std::string m_unit; //!< Parameter unit
167   - double m_factor_value; //!< Parameter value factor
168   - double m_factor_error; //!< Uncertainty in parameter value factor
169   - double m_factor_gradient; //!< Function gradient factor
170   - double m_factor_min; //!< Parameter minimum factor
171   - double m_factor_max; //!< Parameter maximum factor
172   - double m_scale; //!< Parameter scaling (true = factor * scale)
173   - bool m_free; //!< Parameter is free
174   - bool m_has_min; //!< Parameter has minimum boundary
175   - bool m_has_max; //!< Parameter has maximum boundary
176   - bool m_has_grad; //!< Parameter has analytic gradient
  165 + std::string m_name; //!< Parameter name
  166 + std::string m_unit; //!< Parameter unit
  167 + double m_factor_value; //!< Parameter value factor
  168 + double m_factor_error; //!< Uncertainty in parameter value factor
  169 + double m_factor_min; //!< Parameter minimum factor
  170 + double m_factor_max; //!< Parameter maximum factor
  171 + mutable double m_factor_gradient; //!< Function gradient factor
  172 + double m_scale; //!< Parameter scaling (true = factor * scale)
  173 + bool m_free; //!< Parameter is free
  174 + bool m_has_min; //!< Parameter has minimum boundary
  175 + bool m_has_max; //!< Parameter has maximum boundary
  176 + bool m_has_grad; //!< Parameter has analytic gradient
177 177 };
178 178  
179 179  
... ... @@ -376,7 +376,7 @@ void GOptimizerPar::factor_error(const double&amp; error)
376 376 * Sets the parameter gradient factor.
377 377 ***************************************************************************/
378 378 inline
379   -void GOptimizerPar::factor_gradient(const double& gradient)
  379 +void GOptimizerPar::factor_gradient(const double& gradient) const
380 380 {
381 381 m_factor_gradient = gradient;
382 382 return;
... ...
inst/com/include/GCOMModelDRBFitting.hpp
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GCOMModelDRBFitting.hpp - COMPTEL DRB model fitting class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -67,9 +67,8 @@ public:
67 67 virtual std::string type(void) const;
68 68 virtual bool is_constant(void) const;
69 69 virtual double eval(const GEvent& event,
70   - const GObservation& obs) const;
71   - virtual double eval_gradients(const GEvent& event,
72   - const GObservation& obs) const;
  70 + const GObservation& obs,
  71 + const bool& gradients = false) const;
73 72 virtual double npred(const GEnergy& obsEng, const GTime& obsTime,
74 73 const GObservation& obs) const;
75 74 virtual GCOMEventCube* mc(const GObservation& obs, GRan& ran) const;
... ...
inst/com/pyext/GCOMModelDRBFitting.i
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GCOMModelDRBFitting.i - COMPTEL DRB model fitting class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -49,9 +49,8 @@ public:
49 49 virtual std::string type(void) const;
50 50 virtual bool is_constant(void) const;
51 51 virtual double eval(const GEvent& event,
52   - const GObservation& obs) const;
53   - virtual double eval_gradients(const GEvent& event,
54   - const GObservation& obs) const;
  52 + const GObservation& obs,
  53 + const bool& gradients = false) const;
55 54 virtual double npred(const GEnergy& obsEng, const GTime& obsTime,
56 55 const GObservation& obs) const;
57 56 virtual GCOMEventCube* mc(const GObservation& obs, GRan& ran) const;
... ...
inst/com/src/GCOMModelDRBFitting.cpp
  View file @c135e7c
... ... @@ -43,9 +43,7 @@ const GCOMModelDRBFitting g_com_drb_fitting_seed;
43 43 const GModelRegistry g_com_drb_fitting_registry(&g_com_drb_fitting_seed);
44 44  
45 45 /* __ Method name definitions ____________________________________________ */
46   -#define G_EVAL "GCOMModelDRBFitting::eval(GEvent&,GObservation&)"
47   -#define G_EVAL_GRADIENTS "GCOMModelDRBFitting::eval_gradients(GEvent&,"\
48   - "GObservation&)"
  46 +#define G_EVAL "GCOMModelDRBFitting::eval(GEvent&, GObservation&, bool&)"
49 47 #define G_READ "GCOMModelDRBFitting::read(GXmlElement&)"
50 48 #define G_WRITE "GCOMModelDRBFitting::write(GXmlElement&)"
51 49  
... ... @@ -229,6 +227,7 @@ GCOMModelDRBFitting* GCOMModelDRBFitting::clone(void) const
229 227 *
230 228 * @param[in] event Observed event.
231 229 * @param[in] obs Observation.
  230 + * @param[in] gradients Compute gradients?
232 231 * @return Background model value.
233 232 *
234 233 * @exception GException::invalid_argument
... ... @@ -238,7 +237,8 @@ GCOMModelDRBFitting* GCOMModelDRBFitting::clone(void) const
238 237 * Evaluates the COMPTEL DRB fitting model.
239 238 ***************************************************************************/
240 239 double GCOMModelDRBFitting::eval(const GEvent& event,
241   - const GObservation& obs) const
  240 + const GObservation& obs,
  241 + const bool& gradients) const
242 242 {
243 243 // Extract COMPTEL observation
244 244 const GCOMObservation* observation = dynamic_cast<const GCOMObservation*>(&obs);
... ... @@ -262,107 +262,11 @@ double GCOMModelDRBFitting::eval(const GEvent&amp; event,
262 262 // Initialise value
263 263 double value = 0.0;
264 264  
265   - // Get bin index
266   - int index = bin->index();
267   -
268   - // Get bin size.
269   - double size = bin->size();
270   -
271   - // Continue only if bin size is positive
272   - if (size > 0.0) {
273   -
274   - // Initialise scaling factor
275   - double scale = 1.0;
276   -
277   - // Get DRB model value
278   - value = observation->drb().pixels()[index] / size;
279   -
280   - // If model is a scaling factor then use the single parameter as such
281   - if (m_scale) {
282   - scale = m_values[0].value();
283   - }
284   -
285   - // ... otherwise perform a linear interpolation
286   - else {
287   -
288   - // Get Phibar value
289   - double phibar = bin->dir().phibar();
290   -
291   - // Update evaluation cache
292   - update_cache();
293   -
294   - // Set node array for linear interpolation
295   - m_nodes.set_value(phibar);
296   -
297   - // Get scale factor
298   - scale = m_values[m_nodes.inx_left()].value() * m_nodes.wgt_left() +
299   - m_values[m_nodes.inx_right()].value() * m_nodes.wgt_right();
300   -
301   - } // endelse: performed linear interpolation
302   -
303   - // Compute background value
304   - value *= scale;
305   -
306   - // Compile option: Check for NaN/Inf
307   - #if defined(G_NAN_CHECK)
308   - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
309   - std::cout << "*** ERROR: GCOMModelDRBFitting::eval";
310   - std::cout << "(index=" << index << "):";
311   - std::cout << " NaN/Inf encountered";
312   - std::cout << " (value=" << value;
313   - std::cout << ", scale=" << scale;
314   - std::cout << ")" << std::endl;
  265 + // Optionally initialise gradients
  266 + if (gradients) {
  267 + for (int i = 0; i < m_values.size(); ++i) {
  268 + m_values[i].factor_gradient(0.0);
315 269 }
316   - #endif
317   -
318   - } // endif: binsize was positive
319   -
320   - // Return
321   - return value;
322   -}
323   -
324   -
325   -/***********************************************************************//**
326   - * @brief Evaluate function and gradients
327   - *
328   - * @param[in] event Observed event.
329   - * @param[in] obs Observation.
330   - * @return Background model value.
331   - *
332   - * @exception GCOMException::bad_observation_type
333   - * Observation is not a COMPTEL observation.
334   - * @exception GCOMException::bad_event_type
335   - * Event is not a COMPTEL event bin.
336   - *
337   - * Evaluates the COMPTEL DRB fitting model value and sets the parameter
338   - * gradients.
339   - ***************************************************************************/
340   -double GCOMModelDRBFitting::eval_gradients(const GEvent& event,
341   - const GObservation& obs) const
342   -{
343   - // Extract COMPTEL observation
344   - const GCOMObservation* observation = dynamic_cast<const GCOMObservation*>(&obs);
345   - if (observation == NULL) {
346   - std::string cls = std::string(typeid(&obs).name());
347   - std::string msg = "Observation of type \""+cls+"\" is not a COMPTEL "
348   - "observations. Please specify a COMPTEL observation "
349   - "as argument.";
350   - throw GException::invalid_argument(G_EVAL_GRADIENTS, msg);
351   - }
352   -
353   - // Extract COMPTEL event bin
354   - const GCOMEventBin* bin = dynamic_cast<const GCOMEventBin*>(&event);
355   - if (bin == NULL) {
356   - std::string cls = std::string(typeid(&event).name());
357   - std::string msg = "Event of type \""+cls+"\" is not a COMPTEL event. "
358   - "Please specify a COMPTEL event as argument.";
359   - throw GException::invalid_argument(G_EVAL_GRADIENTS, msg);
360   - }
361   -
362   - // Initialise value and gradients
363   - double value = 0.0;
364   - for (int i = 0; i < m_values.size(); ++i) {
365   - const_cast<GCOMModelDRBFitting*>(this)->m_values[i].factor_gradient(0.0);
366 270 }
367 271  
368 272 // Get bin index
... ... @@ -380,20 +284,26 @@ double GCOMModelDRBFitting::eval_gradients(const GEvent&amp; event,
380 284 // Get DRB model value
381 285 value = observation->drb().pixels()[index] / size;
382 286  
383   -
384   - // If model is a scaling factor then use the single parameter as such
  287 + // If the model has a single scaling factor then use the first
  288 + // parameter as scaling factor of the model
385 289 if (m_scale) {
386 290  
387 291 // Get scale factor
388 292 scale = m_values[0].value();
389   -
390   - // Compute partial derivative
391   - double grad = (m_values[0].is_free()) ? value * m_values[0].scale() : 0.0;
392 293  
393   - // Set gradient (circumvent const correctness)
394   - const_cast<GCOMModelDRBFitting*>(this)->m_values[0].factor_gradient(grad);
  294 + // Optionally compute gradients
  295 + if (gradients) {
395 296  
396   - }
  297 + // Compute partial derivative
  298 + double grad = (m_values[0].is_free())
  299 + ? value * m_values[0].scale() : 0.0;
  300 +
  301 + // Set gradient
  302 + m_values[0].factor_gradient(grad);
  303 +
  304 + } // endif: gradient computation was requested
  305 +
  306 + } // endif: model is a scaling factor
397 307  
398 308 // ... otherwise perform a linear interpolation
399 309 else {
... ... @@ -417,17 +327,22 @@ double GCOMModelDRBFitting::eval_gradients(const GEvent&amp; event,
417 327 scale = m_values[inx_left].value() * wgt_left +
418 328 m_values[inx_right].value() * wgt_right;
419 329  
420   - // Gradient for left node
421   - if (m_values[inx_left].is_free()) {
422   - double grad = wgt_left * value * m_values[inx_left].scale();
423   - const_cast<GCOMModelDRBFitting*>(this)->m_values[inx_left].factor_gradient(grad);
424   - }
  330 + // Optionally compute gradients
  331 + if (gradients) {
425 332  
426   - // Gradient for right node
427   - if (m_values[inx_right].is_free()) {
428   - double grad = wgt_right * value * m_values[inx_right].scale();
429   - const_cast<GCOMModelDRBFitting*>(this)->m_values[inx_right].factor_gradient(grad);
430   - }
  333 + // Compute partial derivatives
  334 + double g_left = (m_values[inx_left].is_free())
  335 + ? wgt_left * value * m_values[inx_left].scale()
  336 + : 0.0;
  337 + double g_right = (m_values[inx_right].is_free())
  338 + ? wgt_right * value * m_values[inx_right].scale()
  339 + : 0.0;
  340 +
  341 + // Set gradients
  342 + m_values[inx_left].factor_gradient(g_left);
  343 + m_values[inx_right].factor_gradient(g_right);
  344 +
  345 + } // endif: gradient computation was requested
431 346  
432 347 } // endelse: performed linear interpolation
433 348  
... ... @@ -437,7 +352,7 @@ double GCOMModelDRBFitting::eval_gradients(const GEvent&amp; event,
437 352 // Compile option: Check for NaN/Inf
438 353 #if defined(G_NAN_CHECK)
439 354 if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
440   - std::cout << "*** ERROR: GCOMModelDRBFitting::eval_gradients";
  355 + std::cout << "*** ERROR: GCOMModelDRBFitting::eval";
441 356 std::cout << "(index=" << index << "):";
442 357 std::cout << " NaN/Inf encountered";
443 358 std::cout << " (value=" << value;
... ...
inst/cta/include/GCTAModelAeffBackground.hpp
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GCTAModelAeffBackground.hpp - CTA Aeff background model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2015 by Michael Mayer *
  4 + * copyright (C) 2015-2016 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -64,9 +64,8 @@ public:
64 64 virtual std::string type(void) const;
65 65 virtual bool is_constant(void) const;
66 66 virtual double eval(const GEvent& event,
67   - const GObservation& obs) const;
68   - virtual double eval_gradients(const GEvent& event,
69   - const GObservation& obs) const;
  67 + const GObservation& obs,
  68 + const bool& gradients = false) const;
70 69 virtual double npred(const GEnergy& obsEng,
71 70 const GTime& obsTime,
72 71 const GObservation& obs) const;
... ...
inst/cta/include/GCTAModelCubeBackground.hpp
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GCTAModelCubeBackground.hpp - CTA cube background model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2013-2015 by Michael Mayer *
  4 + * copyright (C) 2013-2016 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -69,9 +69,8 @@ public:
69 69 virtual std::string type(void) const;
70 70 virtual bool is_constant(void) const;
71 71 virtual double eval(const GEvent& event,
72   - const GObservation& obs) const;
73   - virtual double eval_gradients(const GEvent& event,
74   - const GObservation& obs) const;
  72 + const GObservation& obs,
  73 + const bool& gradients = false) const;
75 74 virtual double npred(const GEnergy& obsEng,
76 75 const GTime& obsTime,
77 76 const GObservation& obs) const;
... ...
inst/cta/include/GCTAModelIrfBackground.hpp
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GCTAModelIrfBackground.hpp - CTA IRF background model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2014-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -64,9 +64,8 @@ public:
64 64 virtual std::string type(void) const;
65 65 virtual bool is_constant(void) const;
66 66 virtual double eval(const GEvent& event,
67   - const GObservation& obs) const;
68   - virtual double eval_gradients(const GEvent& event,
69   - const GObservation& obs) const;
  67 + const GObservation& obs,
  68 + const bool& gradients = false) const;
70 69 virtual double npred(const GEnergy& obsEng,
71 70 const GTime& obsTime,
72 71 const GObservation& obs) const;
... ...
inst/cta/include/GCTAModelRadial.hpp
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GCTAModelRadial.hpp - Radial model abstract base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -64,8 +64,8 @@ public:
64 64 virtual GCTAModelRadial* clone(void) const = 0;
65 65 virtual std::string classname(void) const = 0;
66 66 virtual std::string type(void) const = 0;
67   - virtual double eval(const double& offset) const = 0;
68   - virtual double eval_gradients(const double& offset) const = 0;
  67 + virtual double eval(const double& offset,
  68 + const bool& gradients = false) const = 0;
69 69 virtual GCTAInstDir mc(const GCTAInstDir& dir, GRan& ran) const = 0;
70 70 virtual double omega(void) const = 0;
71 71 virtual void read(const GXmlElement& xml) = 0;
... ...
inst/cta/include/GCTAModelRadialAcceptance.hpp
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GCTAModelRadialAcceptance.hpp - Radial acceptance model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -70,9 +70,8 @@ public:
70 70 virtual std::string type(void) const;
71 71 virtual bool is_constant(void) const;
72 72 virtual double eval(const GEvent& event,
73   - const GObservation& obs) const;
74   - virtual double eval_gradients(const GEvent& event,
75   - const GObservation& obs) const;
  73 + const GObservation& obs,
  74 + const bool& gradients = false) const;
76 75 virtual double npred(const GEnergy& obsEng, const GTime& obsTime,
77 76 const GObservation& obs) const;
78 77 virtual GCTAEventList* mc(const GObservation& obs, GRan& ran) const;
... ...
inst/cta/include/GCTAModelRadialGauss.hpp
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GCTAModelRadialGauss.hpp - Radial Gaussian CTA model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -69,8 +69,8 @@ public:
69 69 virtual GCTAModelRadialGauss* clone(void) const;
70 70 virtual std::string classname(void) const;
71 71 virtual std::string type(void) const;
72   - virtual double eval(const double& offset) const;
73   - virtual double eval_gradients(const double& offset) const;
  72 + virtual double eval(const double& offset,
  73 + const bool& gradients = false) const;
74 74 virtual GCTAInstDir mc(const GCTAInstDir& dir, GRan& ran) const;
75 75 virtual double omega(void) const;
76 76 virtual void read(const GXmlElement& xml);
... ... @@ -78,8 +78,8 @@ public:
78 78 virtual std::string print(const GChatter& chatter = NORMAL) const;
79 79  
80 80 // Other methods
81   - double sigma(void) const { return m_sigma.value(); }
82   - void sigma(const double& sigma) { m_sigma.value(sigma); }
  81 + double sigma(void) const;
  82 + void sigma(const double& sigma);
83 83  
84 84 protected:
85 85 // Protected methods
... ... @@ -130,4 +130,29 @@ std::string GCTAModelRadialGauss::type(void) const
130 130 return ("Gaussian");
131 131 }
132 132  
  133 +
  134 +/***********************************************************************//**
  135 + * @brief Return Gaussian width parameter
  136 + *
  137 + * @return Gaussian width parameter.
  138 + ***************************************************************************/
  139 +inline
  140 +double GCTAModelRadialGauss::sigma(void) const
  141 +{
  142 + return (m_sigma.value());
  143 +}
  144 +
  145 +
  146 +/***********************************************************************//**
  147 + * @brief Set Gaussian width parameter
  148 + *
  149 + * @param[in] sigma Gaussian width parameter.
  150 + ***************************************************************************/
  151 +inline
  152 +void GCTAModelRadialGauss::sigma(const double& sigma)
  153 +{
  154 + m_sigma.value(sigma);
  155 + return;
  156 +}
  157 +
133 158 #endif /* GCTAMODELRADIALGAUSS_HPP */
... ...
inst/cta/include/GCTAModelRadialPolynom.hpp
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GCTAModelRadialPolynom.hpp - Radial Polynom CTA model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -70,8 +70,8 @@ public:
70 70 virtual GCTAModelRadialPolynom* clone(void) const;
71 71 virtual std::string classname(void) const;
72 72 virtual std::string type(void) const;
73   - virtual double eval(const double& offset) const;
74   - virtual double eval_gradients(const double& offset) const;
  73 + virtual double eval(const double& offset,
  74 + const bool& gradients = false) const;
75 75 virtual GCTAInstDir mc(const GCTAInstDir& dir, GRan& ran) const;
76 76 virtual double omega(void) const;
77 77 virtual void read(const GXmlElement& xml);
... ...
inst/cta/include/GCTAModelRadialProfile.hpp
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GCTAModelRadialProfile.hpp - Radial Profile CTA model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -71,8 +71,8 @@ public:
71 71 virtual GCTAModelRadialProfile* clone(void) const;
72 72 virtual std::string classname(void) const;
73 73 virtual std::string type(void) const;
74   - virtual double eval(const double& offset) const;
75   - virtual double eval_gradients(const double& offset) const;
  74 + virtual double eval(const double& offset,
  75 + const bool& gradients = false) const;
76 76 virtual GCTAInstDir mc(const GCTAInstDir& dir, GRan& ran) const;
77 77 virtual double omega(void) const;
78 78 virtual void read(const GXmlElement& xml);
... ...
inst/cta/pyext/GCTAModelAeffBackground.i
  View file @c135e7c
1 1 /***************************************************************************
2   - * GCTAModelAeffBackground.i - CTA IRF background model class *
  2 + * GCTAModelAeffBackground.i - CTA Aeff background model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2015-2016 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -44,21 +44,20 @@ public:
44 44 virtual ~GCTAModelAeffBackground(void);
45 45  
46 46 // Implemented pure virtual methods
47   - virtual void clear(void);
  47 + virtual void clear(void);
48 48 virtual GCTAModelAeffBackground* clone(void) const;
49   - virtual std::string classname(void) const;
50   - virtual std::string type(void) const;
51   - virtual bool is_constant(void) const;
52   - virtual double eval(const GEvent& event,
53   - const GObservation& obs) const;
54   - virtual double eval_gradients(const GEvent& event,
55   - const GObservation& obs) const;
56   - virtual double npred(const GEnergy& obsEng,
57   - const GTime& obsTime,
58   - const GObservation& obs) const;
59   - virtual GCTAEventList* mc(const GObservation& obs, GRan& ran) const;
60   - virtual void read(const GXmlElement& xml);
61   - virtual void write(GXmlElement& xml) const;
  49 + virtual std::string classname(void) const;
  50 + virtual std::string type(void) const;
  51 + virtual bool is_constant(void) const;
  52 + virtual double eval(const GEvent& event,
  53 + const GObservation& obs,
  54 + const bool& gradients = false) const;
  55 + virtual double npred(const GEnergy& obsEng,
  56 + const GTime& obsTime,
  57 + const GObservation& obs) const;
  58 + virtual GCTAEventList* mc(const GObservation& obs, GRan& ran) const;
  59 + virtual void read(const GXmlElement& xml);
  60 + virtual void write(GXmlElement& xml) const;
62 61  
63 62 // Other methods
64 63 GModelSpectral* spectral(void) const;
... ...
inst/cta/pyext/GCTAModelCubeBackground.i
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GCTAModelCubeBackground.i - CTA cube background model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2013-2015 by Michael Mayer *
  4 + * copyright (C) 2013-2016 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -52,9 +52,8 @@ public:
52 52 virtual std::string type(void) const;
53 53 virtual bool is_constant(void) const;
54 54 virtual double eval(const GEvent& event,
55   - const GObservation& obs) const;
56   - virtual double eval_gradients(const GEvent& event,
57   - const GObservation& obs) const;
  55 + const GObservation& obs,
  56 + const bool& gradients = false) const;
58 57 virtual double npred(const GEnergy& obsEng, const GTime& obsTime,
59 58 const GObservation& obs) const;
60 59 virtual GCTAEventList* mc(const GObservation& obs, GRan& ran) const;
... ...
inst/cta/pyext/GCTAModelIrfBackground.i
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GCTAModelIrfBackground.i - CTA IRF background model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2014-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -50,9 +50,8 @@ public:
50 50 virtual std::string type(void) const;
51 51 virtual bool is_constant(void) const;
52 52 virtual double eval(const GEvent& event,
53   - const GObservation& obs) const;
54   - virtual double eval_gradients(const GEvent& event,
55   - const GObservation& obs) const;
  53 + const GObservation& obs,
  54 + const bool& gradients = false) const;
56 55 virtual double npred(const GEnergy& obsEng,
57 56 const GTime& obsTime,
58 57 const GObservation& obs) const;
... ...
inst/cta/pyext/GCTAModelRadial.i
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GCTAModelRadial.i - Abstract radial model base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -48,8 +48,8 @@ public:
48 48 virtual GCTAModelRadial* clone(void) const = 0;
49 49 virtual std::string classname(void) const = 0;
50 50 virtual std::string type(void) const = 0;
51   - virtual double eval(const double& offset) const = 0;
52   - virtual double eval_gradients(const double& offset) const = 0;
  51 + virtual double eval(const double& offset,
  52 + const bool& gradients = false) const = 0;
53 53 virtual GCTAInstDir mc(const GCTAInstDir& dir, GRan& ran) const = 0;
54 54 virtual double omega(void) const = 0;
55 55 virtual void read(const GXmlElement& xml) = 0;
... ...
inst/cta/pyext/GCTAModelRadialAcceptance.i
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GCTAModelRadialAcceptance.i - Radial acceptance model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -53,9 +53,8 @@ public:
53 53 virtual std::string type(void) const;
54 54 virtual bool is_constant(void) const;
55 55 virtual double eval(const GEvent& event,
56   - const GObservation& obs) const;
57   - virtual double eval_gradients(const GEvent& event,
58   - const GObservation& obs) const;
  56 + const GObservation& obs,
  57 + const bool& gradients = false) const;
59 58 virtual double npred(const GEnergy& obsEng, const GTime& obsTime,
60 59 const GObservation& obs) const;
61 60 virtual GCTAEventList* mc(const GObservation& obs, GRan& ran) const;
... ...
inst/cta/pyext/GCTAModelRadialGauss.i
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GCTAModelRadialGauss.i - Radial Gaussian model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -49,8 +49,8 @@ public:
49 49 virtual GCTAModelRadialGauss* clone(void) const;
50 50 virtual std::string classname(void) const;
51 51 virtual std::string type(void) const;
52   - virtual double eval(const double& offset) const;
53   - virtual double eval_gradients(const double& offset) const;
  52 + virtual double eval(const double& offset,
  53 + const bool& gradients = false) const;
54 54 virtual GCTAInstDir mc(const GCTAInstDir& dir, GRan& ran) const;
55 55 virtual double omega(void) const;
56 56 virtual void read(const GXmlElement& xml);
... ...
inst/cta/pyext/GCTAModelRadialPolynom.i
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GCTAModelRadialPolynom.i - Radial Polynom model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -49,8 +49,8 @@ public:
49 49 virtual GCTAModelRadialPolynom* clone(void) const;
50 50 virtual std::string classname(void) const;
51 51 virtual std::string type(void) const;
52   - virtual double eval(const double& offset) const;
53   - virtual double eval_gradients(const double& offset) const;
  52 + virtual double eval(const double& offset,
  53 + const bool& gradients = false) const;
54 54 virtual GCTAInstDir mc(const GCTAInstDir& dir, GRan& ran) const;
55 55 virtual double omega(void) const;
56 56 virtual void read(const GXmlElement& xml);
... ...
inst/cta/pyext/GCTAModelRadialProfile.i
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GCTAModelRadialProfile.i - Radial Profile model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -50,8 +50,8 @@ public:
50 50 virtual GCTAModelRadialProfile* clone(void) const;
51 51 virtual std::string classname(void) const;
52 52 virtual std::string type(void) const;
53   - virtual double eval(const double& offset) const;
54   - virtual double eval_gradients(const double& offset) const;
  53 + virtual double eval(const double& offset,
  54 + const bool& gradients = false) const;
55 55 virtual GCTAInstDir mc(const GCTAInstDir& dir, GRan& ran) const;
56 56 virtual double omega(void) const;
57 57 virtual void read(const GXmlElement& xml);
... ...
inst/cta/src/GCTAModelAeffBackground.cpp
  View file @c135e7c
... ... @@ -46,9 +46,7 @@ const GCTAModelAeffBackground g_cta_aeff_background_seed;
46 46 const GModelRegistry g_cta_aeff_background_registry(&g_cta_aeff_background_seed);
47 47  
48 48 /* __ Method name definitions ____________________________________________ */
49   -#define G_EVAL "GCTAModelAeffBackground::eval(GEvent&, GObservation&)"
50   -#define G_EVAL_GRADIENTS "GCTAModelAeffBackground::eval_gradients(GEvent&,"\
51   - " GObservation&)"
  49 +#define G_EVAL "GCTAModelAeffBackground::eval(GEvent&, GObservation&, bool&)"
52 50 #define G_NPRED "GCTAModelAeffBackground::npred(GEnergy&, GTime&,"\
53 51 " GObservation&)"
54 52 #define G_MC "GCTAModelAeffBackground::mc(GObservation&, GRan&)"
... ... @@ -274,15 +272,20 @@ GCTAModelAeffBackground* GCTAModelAeffBackground::clone(void) const
274 272 *
275 273 * @param[in] event Observed event.
276 274 * @param[in] obs Observation.
  275 + * @param[in] gradients Compute gradients?
277 276 * @return Function value.
278 277 *
279 278 * @exception GException::invalid_argument
280 279 * Specified observation is not of the expected type.
281 280 *
  281 + * If the @p gradients flag is true the method will also set the parameter
  282 + * gradients of the model parameters.
  283 + *
282 284 * @todo Make sure that DETX and DETY are always set in GCTAInstDir.
283 285 ***************************************************************************/
284 286 double GCTAModelAeffBackground::eval(const GEvent& event,
285   - const GObservation& obs) const
  287 + const GObservation& obs,
  288 + const bool& gradients) const
286 289 {
287 290 // Get pointer on CTA observation
288 291 const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>(&obs);
... ... @@ -330,85 +333,10 @@ double GCTAModelAeffBackground::eval(const GEvent&amp; event,
330 333 cta->pointing().zenith(),
331 334 cta->pointing().azimuth(), false);
332 335 double spec = (spectral() != NULL)
333   - ? spectral()->eval(event.energy(), event.time()) : 1.0;
  336 + ? spectral()->eval(event.energy(), event.time(), gradients)
  337 + : 1.0;
334 338 double temp = (temporal() != NULL)
335   - ? temporal()->eval(event.time()) : 1.0;
336   -
337   - // Compute value
338   - double value = spat * spec * temp;
339   -
340   - // Apply deadtime correction
341   - value *= obs.deadc(event.time());
342   -
343   - // Return value
344   - return value;
345   -}
346   -
347   -
348   -/***********************************************************************//**
349   - * @brief Evaluate function and gradients
350   - *
351   - * @param[in] event Observed event.
352   - * @param[in] obs Observation.
353   - * @return Function value.
354   - *
355   - * @exception GException::invalid_argument
356   - * Specified observation is not of the expected type.
357   - *
358   - * @todo Make sure that DETX and DETY are always set in GCTAInstDir.
359   - ***************************************************************************/
360   -double GCTAModelAeffBackground::eval_gradients(const GEvent& event,
361   - const GObservation& obs) const
362   -{
363   - // Get pointer on CTA observation
364   - const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>(&obs);
365   - if (cta == NULL) {
366   - std::string msg = "Specified observation is not a CTA observation.\n" +
367   - obs.print();
368   - throw GException::invalid_argument(G_EVAL_GRADIENTS, msg);
369   - }
370   -
371   - // Get pointer on CTA IRF response
372   - const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(cta->response());
373   - if (rsp == NULL) {
374   - std::string msg = "Specified observation does not contain an"
375   - " IRF response.\n" + obs.print();
376   - throw GException::invalid_argument(G_EVAL_GRADIENTS, msg);
377   - }
378   -
379   - // Retrieve pointer to CTA Effective Area
380   - const GCTAAeff* aeff = rsp->aeff();
381   - if (aeff == NULL) {
382   - std::string msg = "Specified observation contains no effective area"
383   - " information.\n" + obs.print();
384   - throw GException::invalid_argument(G_EVAL, msg);
385   - }
386   -
387   - // Extract CTA instrument direction from event
388   - const GCTAInstDir* dir = dynamic_cast<const GCTAInstDir*>(&(event.dir()));
389   - if (dir == NULL) {
390   - std::string msg = "No CTA instrument direction found in event.";
391   - throw GException::invalid_argument(G_EVAL, msg);
392   - }
393   -
394   - // Set DETX and DETY in instrument direction
395   - GCTAInstDir inst_dir = cta->pointing().instdir(dir->dir());
396   -
397   - // Set theta and phi from instrument coordinates
398   - double theta = std::sqrt(inst_dir.detx() * inst_dir.detx() +
399   - inst_dir.dety() * inst_dir.dety());
400   - double phi = gammalib::atan2d(inst_dir.dety(), inst_dir.detx()) *
401   - gammalib::deg2rad;
402   -
403   - // Evaluate function
404   - double logE = event.energy().log10TeV();
405   - double spat = (*aeff)(logE, theta, phi,
406   - cta->pointing().zenith(),
407   - cta->pointing().azimuth(), false);
408   - double spec = (spectral() != NULL)
409   - ? spectral()->eval_gradients(event.energy(), event.time()) : 1.0;
410   - double temp = (temporal() != NULL)
411   - ? temporal()->eval_gradients(event.time()) : 1.0;
  339 + ? temporal()->eval(event.time(), gradients) : 1.0;
412 340  
413 341 // Compute value
414 342 double value = spat * spec * temp;
... ... @@ -417,23 +345,28 @@ double GCTAModelAeffBackground::eval_gradients(const GEvent&amp; event,
417 345 double deadc = obs.deadc(event.time());
418 346 value *= deadc;
419 347  
420   - // Multiply factors to spectral gradients
421   - if (spectral() != NULL) {
422   - double fact = spat * temp * deadc;
423   - if (fact != 1.0) {
424   - for (int i = 0; i < spectral()->size(); ++i)
425   - (*spectral())[i].factor_gradient( (*spectral())[i].factor_gradient() * fact );
  348 + // Optionally compute partial derivatives
  349 + if (gradients) {
  350 +
  351 + // Multiply factors to spectral gradients
  352 + if (spectral() != NULL) {
  353 + double fact = spat * temp * deadc;
  354 + if (fact != 1.0) {
  355 + for (int i = 0; i < spectral()->size(); ++i)
  356 + (*spectral())[i].factor_gradient((*spectral())[i].factor_gradient() * fact );
  357 + }
426 358 }
427   - }
428 359  
429   - // Multiply factors to temporal gradients
430   - if (temporal() != NULL) {
431   - double fact = spat * spec * deadc;
432   - if (fact != 1.0) {
433   - for (int i = 0; i < temporal()->size(); ++i)
434   - (*temporal())[i].factor_gradient( (*temporal())[i].factor_gradient() * fact );
  360 + // Multiply factors to temporal gradients
  361 + if (temporal() != NULL) {
  362 + double fact = spat * spec * deadc;
  363 + if (fact != 1.0) {
  364 + for (int i = 0; i < temporal()->size(); ++i)
  365 + (*temporal())[i].factor_gradient((*temporal())[i].factor_gradient() * fact );
  366 + }
435 367 }
436   - }
  368 +
  369 + } // endif: computed partial derivatives
437 370  
438 371 // Return value
439 372 return value;
... ...
inst/cta/src/GCTAModelCubeBackground.cpp
  View file @c135e7c
... ... @@ -44,9 +44,7 @@ const GCTAModelCubeBackground g_cta_inst_background_seed;
44 44 const GModelRegistry g_cta_inst_background_registry(&g_cta_inst_background_seed);
45 45  
46 46 /* __ Method name definitions ____________________________________________ */
47   -#define G_EVAL "GCTAModelCubeBackground::eval(GEvent&, GObservation&)"
48   -#define G_EVAL_GRADIENTS "GCTAModelCubeBackground::eval_gradients(GEvent&,"\
49   - " GObservation&)"
  47 +#define G_EVAL "GCTAModelCubeBackground::eval(GEvent&, GObservation&, bool&)"
50 48 #define G_NPRED "GCTAModelCubeBackground::npred(GEnergy&, GTime&,"\
51 49 " GObservation&)"
52 50 #define G_MC "GCTAModelCubeBackground::mc(GObservation&, GRan&)"
... ... @@ -251,13 +249,18 @@ GCTAModelCubeBackground* GCTAModelCubeBackground::clone(void) const
251 249 *
252 250 * @param[in] event Observed event.
253 251 * @param[in] obs Observation.
  252 + * @param[in] gradients Compute gradients?
254 253 * @return Function value.
255 254 *
256 255 * @exception GException::invalid_argument
257 256 * Specified observation is not of the expected type.
  257 + *
  258 + * If the @p gradients flag is true the method will also set the parameter
  259 + * gradients of the model parameters.
258 260 ***************************************************************************/
259 261 double GCTAModelCubeBackground::eval(const GEvent& event,
260   - const GObservation& obs) const
  262 + const GObservation& obs,
  263 + const bool& gradients) const
261 264 {
262 265 // Get pointer on CTA observation
263 266 const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>(&obs);
... ... @@ -287,89 +290,38 @@ double GCTAModelCubeBackground::eval(const GEvent&amp; event,
287 290 //double logE = event.energy().log10TeV();
288 291 double spat = bgd((*dir), event.energy());
289 292 double spec = (spectral() != NULL)
290   - ? spectral()->eval(event.energy(), event.time()) : 1.0;
291   - double temp = (temporal() != NULL)
292   - ? temporal()->eval(event.time()) : 1.0;
293   -
294   - // Compute value. Note that background rates are already per
295   - // livetime, hence no deadtime correction is needed here.
296   - double value = spat * spec * temp;
297   -
298   - // Return value
299   - return value;
300   -}
301   -
302   -
303   -/***********************************************************************//**
304   - * @brief Evaluate function and gradients
305   - *
306   - * @param[in] event Observed event.
307   - * @param[in] obs Observation.
308   - * @return Function value.
309   - *
310   - * @exception GException::invalid_argument
311   - * Specified observation is not of the expected type.
312   - ***************************************************************************/
313   -double GCTAModelCubeBackground::eval_gradients(const GEvent& event,
314   - const GObservation& obs) const
315   -{
316   - // Get pointer on CTA observation
317   - const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>(&obs);
318   - if (cta == NULL) {
319   - std::string msg = "Specified observation is not a CTA observation.\n" +
320   - obs.print();
321   - throw GException::invalid_argument(G_EVAL_GRADIENTS, msg);
322   - }
323   -
324   - // Get pointer on CTA IRF response
325   - const GCTAResponseCube* rsp = dynamic_cast<const GCTAResponseCube*>(cta->response());
326   - if (rsp == NULL) {
327   - std::string msg = "Specified observation does not contain a"
328   - " cube response.\n" + obs.print();
329   - throw GException::invalid_argument(G_EVAL_GRADIENTS, msg);
330   - }
331   -
332   - // Extract CTA instrument direction from event
333   - const GCTAInstDir* dir = dynamic_cast<const GCTAInstDir*>(&(event.dir()));
334   - if (dir == NULL) {
335   - std::string msg = "No CTA instrument direction found in event.";
336   - throw GException::invalid_argument(G_EVAL_GRADIENTS, msg);
337   - }
338   -
339   - // Retrieve reference to CTA cube background
340   - const GCTACubeBackground& bgd = rsp->background();
341   -
342   - // Evaluate function
343   - //double logE = event.energy().log10TeV();
344   - double spat = bgd((*dir), event.energy());
345   - double spec = (spectral() != NULL)
346   - ? spectral()->eval_gradients(event.energy(), event.time())
  293 + ? spectral()->eval(event.energy(), event.time(), gradients)
347 294 : 1.0;
348 295 double temp = (temporal() != NULL)
349   - ? temporal()->eval_gradients(event.time())
350   - : 1.0;
  296 + ? temporal()->eval(event.time(), gradients) : 1.0;
351 297  
352 298 // Compute value. Note that background rates are already per
353 299 // livetime, hence no deadtime correction is needed here.
354 300 double value = spat * spec * temp;
355 301  
356   - // Multiply factors to spectral gradients
357   - if (spectral() != NULL) {
358   - double fact = spat * temp;
359   - if (fact != 1.0) {
360   - for (int i = 0; i < spectral()->size(); ++i)
361   - (*spectral())[i].factor_gradient( (*spectral())[i].factor_gradient() * fact );
  302 + // If gradients were requested then multiply factors to spectral and
  303 + // temporal gradients
  304 + if (gradients) {
  305 +
  306 + // Multiply factors to spectral gradients
  307 + if (spectral() != NULL) {
  308 + double fact = spat * temp;
  309 + if (fact != 1.0) {
  310 + for (int i = 0; i < spectral()->size(); ++i)
  311 + (*spectral())[i].factor_gradient((*spectral())[i].factor_gradient() * fact );
  312 + }
362 313 }
363   - }
364 314  
365   - // Multiply factors to temporal gradients
366   - if (temporal() != NULL) {
367   - double fact = spat * spec;
368   - if (fact != 1.0) {
369   - for (int i = 0; i < temporal()->size(); ++i)
370   - (*temporal())[i].factor_gradient( (*temporal())[i].factor_gradient() * fact );
  315 + // Multiply factors to temporal gradients
  316 + if (temporal() != NULL) {
  317 + double fact = spat * spec;
  318 + if (fact != 1.0) {
  319 + for (int i = 0; i < temporal()->size(); ++i)
  320 + (*temporal())[i].factor_gradient((*temporal())[i].factor_gradient() * fact );
  321 + }
371 322 }
372   - }
  323 +
  324 + } // endif: gradients were requested
373 325  
374 326 // Return value
375 327 return value;
... ...
inst/cta/src/GCTAModelIrfBackground.cpp
  View file @c135e7c
... ... @@ -46,9 +46,7 @@ const GCTAModelIrfBackground g_cta_inst_background_seed;
46 46 const GModelRegistry g_cta_inst_background_registry(&g_cta_inst_background_seed);
47 47  
48 48 /* __ Method name definitions ____________________________________________ */
49   -#define G_EVAL "GCTAModelIrfBackground::eval(GEvent&, GObservation&)"
50   -#define G_EVAL_GRADIENTS "GCTAModelIrfBackground::eval_gradients(GEvent&,"\
51   - " GObservation&)"
  49 +#define G_EVAL "GCTAModelIrfBackground::eval(GEvent&, GObservation&, bool&)"
52 50 #define G_NPRED "GCTAModelIrfBackground::npred(GEnergy&, GTime&,"\
53 51 " GObservation&)"
54 52 #define G_MC "GCTAModelIrfBackground::mc(GObservation&, GRan&)"
... ... @@ -274,15 +272,20 @@ GCTAModelIrfBackground* GCTAModelIrfBackground::clone(void) const
274 272 *
275 273 * @param[in] event Observed event.
276 274 * @param[in] obs Observation.
  275 + * @param[in] gradients Compute gradients?
277 276 * @return Function value.
278 277 *
279 278 * @exception GException::invalid_argument
280 279 * Specified observation is not of the expected type.
281 280 *
  281 + * If the @p gradients flag is true the method will also set the parameter
  282 + * gradients of the model parameters.
  283 + *
282 284 * @todo Make sure that DETX and DETY are always set in GCTAInstDir.
283 285 ***************************************************************************/
284   -double GCTAModelIrfBackground::eval(const GEvent& event,
285   - const GObservation& obs) const
  286 +double GCTAModelIrfBackground::eval(const GEvent& event,
  287 + const GObservation& obs,
  288 + const bool& gradients) const
286 289 {
287 290 // Get pointer on CTA observation
288 291 const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>(&obs);
... ... @@ -322,79 +325,10 @@ double GCTAModelIrfBackground::eval(const GEvent&amp; event,
322 325 double logE = event.energy().log10TeV();
323 326 double spat = (*bgd)(logE, inst_dir.detx(), inst_dir.dety());
324 327 double spec = (spectral() != NULL)
325   - ? spectral()->eval(event.energy(), event.time()) : 1.0;
326   - double temp = (temporal() != NULL)
327   - ? temporal()->eval(event.time()) : 1.0;
328   -
329   - // Compute value
330   - double value = spat * spec * temp;
331   -
332   - // Apply deadtime correction
333   - value *= obs.deadc(event.time());
334   -
335   - // Return value
336   - return value;
337   -}
338   -
339   -
340   -/***********************************************************************//**
341   - * @brief Evaluate function and gradients
342   - *
343   - * @param[in] event Observed event.
344   - * @param[in] obs Observation.
345   - * @return Function value.
346   - *
347   - * @exception GException::invalid_argument
348   - * Specified observation is not of the expected type.
349   - *
350   - * @todo Make sure that DETX and DETY are always set in GCTAInstDir.
351   - ***************************************************************************/
352   -double GCTAModelIrfBackground::eval_gradients(const GEvent& event,
353   - const GObservation& obs) const
354   -{
355   - // Get pointer on CTA observation
356   - const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>(&obs);
357   - if (cta == NULL) {
358   - std::string msg = "Specified observation is not a CTA observation.\n" +
359   - obs.print();
360   - throw GException::invalid_argument(G_EVAL_GRADIENTS, msg);
361   - }
362   -
363   - // Get pointer on CTA IRF response
364   - const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(cta->response());
365   - if (rsp == NULL) {
366   - std::string msg = "Specified observation does not contain an"
367   - " IRF response.\n" + obs.print();
368   - throw GException::invalid_argument(G_EVAL_GRADIENTS, msg);
369   - }
370   -
371   - // Retrieve pointer to CTA background
372   - const GCTABackground* bgd = rsp->background();
373   - if (bgd == NULL) {
374   - std::string msg = "Specified observation contains no background"
375   - " information.\n" + obs.print();
376   - throw GException::invalid_argument(G_EVAL_GRADIENTS, msg);
377   - }
378   -
379   - // Extract CTA instrument direction from event
380   - const GCTAInstDir* dir = dynamic_cast<const GCTAInstDir*>(&(event.dir()));
381   - if (dir == NULL) {
382   - std::string msg = "No CTA instrument direction found in event.";
383   - throw GException::invalid_argument(G_EVAL_GRADIENTS, msg);
384   - }
385   -
386   - // Set DETX and DETY in instrument direction
387   - GCTAInstDir inst_dir = cta->pointing().instdir(dir->dir());
388   -
389   - // Evaluate function
390   - double logE = event.energy().log10TeV();
391   - double spat = (*bgd)(logE, inst_dir.detx(), inst_dir.dety());
392   - double spec = (spectral() != NULL)
393   - ? spectral()->eval_gradients(event.energy(), event.time())
  328 + ? spectral()->eval(event.energy(), event.time(), gradients)
394 329 : 1.0;
395 330 double temp = (temporal() != NULL)
396   - ? temporal()->eval_gradients(event.time())
397   - : 1.0;
  331 + ? temporal()->eval(event.time(), gradients) : 1.0;
398 332  
399 333 // Compute value
400 334 double value = spat * spec * temp;
... ... @@ -403,23 +337,29 @@ double GCTAModelIrfBackground::eval_gradients(const GEvent&amp; event,
403 337 double deadc = obs.deadc(event.time());
404 338 value *= deadc;
405 339  
406   - // Multiply factors to spectral gradients
407   - if (spectral() != NULL) {
408   - double fact = spat * temp * deadc;
409   - if (fact != 1.0) {
410   - for (int i = 0; i < spectral()->size(); ++i)
411   - (*spectral())[i].factor_gradient( (*spectral())[i].factor_gradient() * fact );
  340 + // If gradients were requested then multiply factors to spectral and
  341 + // temporal gradients
  342 + if (gradients) {
  343 +
  344 + // Multiply factors to spectral gradients
  345 + if (spectral() != NULL) {
  346 + double fact = spat * temp * deadc;
  347 + if (fact != 1.0) {
  348 + for (int i = 0; i < spectral()->size(); ++i)
  349 + (*spectral())[i].factor_gradient((*spectral())[i].factor_gradient() * fact );
  350 + }
412 351 }
413   - }
414 352  
415   - // Multiply factors to temporal gradients
416   - if (temporal() != NULL) {
417   - double fact = spat * spec * deadc;
418   - if (fact != 1.0) {
419   - for (int i = 0; i < temporal()->size(); ++i)
420   - (*temporal())[i].factor_gradient( (*temporal())[i].factor_gradient() * fact );
  353 + // Multiply factors to temporal gradients
  354 + if (temporal() != NULL) {
  355 + double fact = spat * spec * deadc;
  356 + if (fact != 1.0) {
  357 + for (int i = 0; i < temporal()->size(); ++i)
  358 + (*temporal())[i].factor_gradient((*temporal())[i].factor_gradient() * fact );
  359 + }
421 360 }
422   - }
  361 +
  362 + } // endif: gradients were requested
423 363  
424 364 // Return value
425 365 return value;
... ...
inst/cta/src/GCTAModelRadialAcceptance.cpp
  View file @c135e7c
... ... @@ -53,10 +53,8 @@ const GModelRegistry g_cta_radial_acceptance_registry(&amp;g_cta_radial_a
53 53  
54 54 /* __ Method name definitions ____________________________________________ */
55 55 #define G_ACCESS "GCTAModelRadialAcceptance::operator() (int)"
56   -#define G_EVAL "GCTAModelRadialAcceptance::eval(GEvent&,"\
57   - " GObservation&)"
58   -#define G_EVAL_GRADIENTS "GCTAModelRadialAcceptance::eval_gradients(GEvent&,"\
59   - " GObservation&)"
  56 +#define G_EVAL "GCTAModelRadialAcceptance::eval(GEvent&, GObservation&, "\
  57 + "bool&)"
60 58 #define G_NPRED "GCTAModelRadialAcceptance::npred(GEnergy&, GTime&,"\
61 59 " GObservation&)"
62 60 #define G_MC "GCTAModelRadialAcceptance::mc(GObservation&, GRan&)"
... ... @@ -282,77 +280,18 @@ GCTAModelRadialAcceptance* GCTAModelRadialAcceptance::clone(void) const
282 280 *
283 281 * @param[in] event Observed event.
284 282 * @param[in] obs Observation.
  283 + * @param[in] gradients Compute gradients?
  284 + * @return Function value.
285 285 *
286 286 * Evaluates tha CTA radial acceptance model which is a factorization of a
287 287 * spatial, spectral and temporal model component. This method also applies
288 288 * a deadtime correction factor, so that the normalization of the model is
289 289 * a real rate (counts/exposure time).
290 290 *
291   - * @todo Add bookkeeping of last value and evaluate only if argument
292   - * changed
293   - * @todo Verify that CTA instrument direction pointer is valid, or better,
294   - * add an offset method to GCTAPointing. Ideally, we should precompute
295   - * all offset angles (for an event cube this may only be done easily
296   - * if the pointing has been fixed; otherwise we need a structure
297   - * similar to the Fermi/LAT livetime cube that provides the effective
298   - * sky exposure as function of offset angle).
299   - ***************************************************************************/
300   -double GCTAModelRadialAcceptance::eval(const GEvent& event,
301   - const GObservation& obs) const
302   -{
303   - // Get pointer on CTA observation
304   - const GCTAObservation* ctaobs = dynamic_cast<const GCTAObservation*>(&obs);
305   - if (ctaobs == NULL) {
306   - std::string msg = "Specified observation is not a CTA observation.\n" +
307   - obs.print();
308   - throw GException::invalid_argument(G_EVAL, msg);
309   - }
310   -
311   - // Get pointer on CTA pointing
312   - const GCTAPointing& pnt = ctaobs->pointing();
313   -
314   - // Get instrument direction
315   - const GInstDir* inst_dir = &(event.dir());
316   - const GCTAInstDir* cta_dir = static_cast<const GCTAInstDir*>(inst_dir);
317   -
318   - // Compute offset angle (in degrees)
319   - double offset = cta_dir->dir().dist_deg(pnt.dir());
320   -
321   - // Evaluate function and gradients
322   - double rad = (radial() != NULL)
323   - ? radial()->eval(offset) : 1.0;
324   - double spec = (spectral() != NULL)
325   - ? spectral()->eval(event.energy(), event.time()) : 1.0;
326   - double temp = (temporal() != NULL)
327   - ? temporal()->eval(event.time()) : 1.0;
328   -
329   - // Compute value
330   - double value = rad * spec * temp;
331   -
332   - // Apply deadtime correction
333   - value *= obs.deadc(event.time());
334   -
335   - // Return
336   - return value;
337   -}
338   -
339   -
340   -/***********************************************************************//**
341   - * @brief Evaluate function and gradients
342   - *
343   - * @param[in] event Observed event.
344   - * @param[in] obs Observation (not used).
  291 + * If the @p gradients flag is true the method will also set the parameter
  292 + * gradients of the model parameters.
345 293 *
346   - * @exception GCTAException::no_pointing
347   - * No valid CTA pointing found in observation
348   - *
349   - * Evaluates tha CTA radial acceptance model and parameter gradients. The CTA
350   - * radial acceptance model is a factorization of a spatial, spectral and
351   - * temporal model component. This method also applies a deadtime correction
352   - * factor, so that the normalization of the model is a real rate
353   - * (counts/exposure time).
354   - *
355   - * @todo Add bookkeeping of last value and evaluate only if argument
  294 + * @todo Add bookkeeping of last value and evaluate only if argument
356 295 * changed
357 296 * @todo Verify that CTA instrument direction pointer is valid, or better,
358 297 * add an offset method to GCTAPointing. Ideally, we should precompute
... ... @@ -361,15 +300,16 @@ double GCTAModelRadialAcceptance::eval(const GEvent&amp; event,
361 300 * similar to the Fermi/LAT livetime cube that provides the effective
362 301 * sky exposure as function of offset angle).
363 302 ***************************************************************************/
364   -double GCTAModelRadialAcceptance::eval_gradients(const GEvent& event,
365   - const GObservation& obs) const
  303 +double GCTAModelRadialAcceptance::eval(const GEvent& event,
  304 + const GObservation& obs,
  305 + const bool& gradients) const
366 306 {
367 307 // Get pointer on CTA observation
368 308 const GCTAObservation* ctaobs = dynamic_cast<const GCTAObservation*>(&obs);
369 309 if (ctaobs == NULL) {
370 310 std::string msg = "Specified observation is not a CTA observation.\n" +
371 311 obs.print();
372   - throw GException::invalid_argument(G_EVAL_GRADIENTS, msg);
  312 + throw GException::invalid_argument(G_EVAL, msg);
373 313 }
374 314  
375 315 // Get pointer on CTA pointing
... ... @@ -384,11 +324,12 @@ double GCTAModelRadialAcceptance::eval_gradients(const GEvent&amp; event,
384 324  
385 325 // Evaluate function and gradients
386 326 double rad = (radial() != NULL)
387   - ? radial()->eval_gradients(offset) : 1.0;
  327 + ? radial()->eval(offset, gradients) : 1.0;
388 328 double spec = (spectral() != NULL)
389   - ? spectral()->eval_gradients(event.energy(), event.time()) : 1.0;
  329 + ? spectral()->eval(event.energy(), event.time(), gradients)
  330 + : 1.0;
390 331 double temp = (temporal() != NULL)
391   - ? temporal()->eval_gradients(event.time()) : 1.0;
  332 + ? temporal()->eval(event.time(), gradients) : 1.0;
392 333  
393 334 // Compute value
394 335 double value = rad * spec * temp;
... ... @@ -397,34 +338,39 @@ double GCTAModelRadialAcceptance::eval_gradients(const GEvent&amp; event,
397 338 double deadc = obs.deadc(event.time());
398 339 value *= deadc;
399 340  
400   - // Multiply factors to radial gradients
401   - if (radial() != NULL) {
402   - double fact = spec * temp * deadc;
403   - if (fact != 1.0) {
404   - for (int i = 0; i < radial()->size(); ++i)
405   - (*radial())[i].factor_gradient( (*radial())[i].factor_gradient() * fact );
  341 + // Optionally compute partial derivatives
  342 + if (gradients) {
  343 +
  344 + // Multiply factors to radial gradients
  345 + if (radial() != NULL) {
  346 + double fact = spec * temp * deadc;
  347 + if (fact != 1.0) {
  348 + for (int i = 0; i < radial()->size(); ++i)
  349 + (*radial())[i].factor_gradient((*radial())[i].factor_gradient() * fact );
  350 + }
406 351 }
407   - }
408 352  
409   - // Multiply factors to spectral gradients
410   - if (spectral() != NULL) {
411   - double fact = rad * temp * deadc;
412   - if (fact != 1.0) {
413   - for (int i = 0; i < spectral()->size(); ++i)
414   - (*spectral())[i].factor_gradient( (*spectral())[i].factor_gradient() * fact );
  353 + // Multiply factors to spectral gradients
  354 + if (spectral() != NULL) {
  355 + double fact = rad * temp * deadc;
  356 + if (fact != 1.0) {
  357 + for (int i = 0; i < spectral()->size(); ++i)
  358 + (*spectral())[i].factor_gradient((*spectral())[i].factor_gradient() * fact );
  359 + }
415 360 }
416   - }
417 361  
418   - // Multiply factors to temporal gradients
419   - if (temporal() != NULL) {
420   - double fact = rad * spec * deadc;
421   - if (fact != 1.0) {
422   - for (int i = 0; i < temporal()->size(); ++i)
423   - (*temporal())[i].factor_gradient( (*temporal())[i].factor_gradient() * fact );
  362 + // Multiply factors to temporal gradients
  363 + if (temporal() != NULL) {
  364 + double fact = rad * spec * deadc;
  365 + if (fact != 1.0) {
  366 + for (int i = 0; i < temporal()->size(); ++i)
  367 + (*temporal())[i].factor_gradient((*temporal())[i].factor_gradient() * fact );
  368 + }
424 369 }
425   - }
426 370  
427   - // Return value
  371 + } // endif: computed partial derivatives
  372 +
  373 + // Return
428 374 return value;
429 375 }
430 376  
... ...
inst/cta/src/GCTAModelRadialGauss.cpp
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GCTAModelRadialGauss.cpp - Radial Gaussian CTA model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -217,6 +217,8 @@ GCTAModelRadialGauss* GCTAModelRadialGauss::clone(void) const
217 217 * @brief Evaluate function
218 218 *
219 219 * @param[in] offset Offset angle [degrees].
  220 + * @param[in] gradients Compute gradients?
  221 + * @return Function value
220 222 *
221 223 * Evaluates the Gaussian model for a given offset. The Gaussian model is
222 224 * defined as
... ... @@ -226,66 +228,36 @@ GCTAModelRadialGauss* GCTAModelRadialGauss::clone(void) const
226 228 * \f$\theta\f$ is the offset angle (in degrees), and
227 229 * \f$\sigma\f$ is the Gaussian width (in degrees\f$^2\f$).
228 230 *
229   - * Note that this method implements a function which is unity for
230   - * \f$\theta=0\f$.
231   - ***************************************************************************/
232   -double GCTAModelRadialGauss::eval(const double& offset) const
233   -{
234   - // Compute value
235   - double arg = offset * offset / sigma();
236   - double value = std::exp(-0.5 * arg * arg);
237   -
238   - // Compile option: Check for NaN/Inf
239   - #if defined(G_NAN_CHECK)
240   - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
241   - std::cout << "*** ERROR: GCTAModelRadialGauss::eval";
242   - std::cout << "(offset=" << offset << "): NaN/Inf encountered";
243   - std::cout << " (value=" << value;
244   - std::cout << ", sigma=" << sigma();
245   - std::cout << ")" << std::endl;
246   - }
247   - #endif
248   -
249   - // Return value
250   - return value;
251   -}
252   -
253   -
254   -/***********************************************************************//**
255   - * @brief Evaluate function and gradients
256   - *
257   - * @param[in] offset Offset angle [degrees].
258   - *
259   - * Evaluates the Gaussian model for a given offset. The Gaussian model is
260   - * defined as
261   - * \f[f(\theta) = \exp \left(-\frac{1}{2}
262   - * \left( \frac{\theta^2}{\sigma} \right)^2 \right)\f]
263   - * where
264   - * \f$\theta\f$ is the offset angle (in degrees), and
265   - * \f$\sigma\f$ is the Gaussian width (in degrees\f$^2\f$).
266   - *
267   - * The partial derivative of the Gaussian width is given by
  231 + * If the @p gradients flag is true the method will also compute the partial
  232 + * derivatives of the parameters. The partial derivative of the Gaussian width
  233 + * is given by
268 234 * \f[\frac{df}{d\sigma_v} = f(\theta) \frac{\theta^4}{\sigma^3} \sigma_s\f]
269 235 * where
270 236 * \f$\sigma_v\f$ is the value part,
271 237 * \f$\sigma_s\f$ is the scaling part, and
272   - * \f$\sigma = \sigma_v \sigma_s\f$.
  238 + * \f$\sigma = \sigma_v \sigma_s\f$.
273 239 *
274 240 * Note that this method implements a function which is unity for
275 241 * \f$\theta=0\f$.
276 242 ***************************************************************************/
277   -double GCTAModelRadialGauss::eval_gradients(const double& offset) const
  243 +double GCTAModelRadialGauss::eval(const double& offset,
  244 + const bool& gradients) const
278 245 {
279 246 // Compute value
280 247 double arg = offset * offset / sigma();
281 248 double arg2 = arg * arg;
282   - double value = exp(-0.5 * arg2);
  249 + double value = std::exp(-0.5 * arg2);
  250 +
  251 + // Optionally compute partial derivatives
  252 + if (gradients) {
  253 +
  254 + // Compute partial derivatives of the sigma parameter.
  255 + double g_sigma = value * arg2 / sigma() * m_sigma.scale();
283 256  
284   - // Compute partial derivatives of the sigma parameter.
285   - double g_sigma = value * arg2 / sigma() * m_sigma.scale();
  257 + // Set factor gradient
  258 + m_sigma.factor_gradient(g_sigma);
286 259  
287   - // Set gradients (circumvent const correctness)
288   - const_cast<GCTAModelRadialGauss*>(this)->m_sigma.factor_gradient(g_sigma);
  260 + } // endif: computed partial derivatives
289 261  
290 262 // Compile option: Check for NaN/Inf
291 263 #if defined(G_NAN_CHECK)
... ...
inst/cta/src/GCTAModelRadialPolynom.cpp
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GCTAModelRadialPolynom.cpp - Radial Polynom CTA model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -243,6 +243,8 @@ GCTAModelRadialPolynom* GCTAModelRadialPolynom::clone(void) const
243 243 * @brief Evaluate function
244 244 *
245 245 * @param[in] offset Offset angle [degrees].
  246 + * @param[in] gradients Compute gradients?
  247 + * @return Function value
246 248 *
247 249 * Evaluates the radial polynomial model for a given offset. The model is
248 250 * defined as
... ... @@ -250,69 +252,18 @@ GCTAModelRadialPolynom* GCTAModelRadialPolynom::clone(void) const
250 252 * where
251 253 * \f$\theta\f$ is the offset angle (in degrees), and
252 254 * \f$c_i\f$ are the polynomial coefficients.
253   - ***************************************************************************/
254   -double GCTAModelRadialPolynom::eval(const double& offset) const
255   -{
256   - // Initialise result
257   - double value = 0.0;
258   -
259   - // Determine polynomial degree
260   - int ncoeffs = m_coeffs.size();
261   -
262   - // Compute value and gradients (only if there are coefficients)
263   - if (ncoeffs > 0) {
264   -
265   - // Compute value
266   - value = m_coeffs[ncoeffs-1].value();
267   - for (int i = ncoeffs-2; i >= 0; i--) {
268   - value = value * offset + m_coeffs[i].value();
269   - }
270   -
271   - } // endif: there were coefficients
272   -
273   - // Compile option: Check for NaN/Inf
274   - #if defined(G_NAN_CHECK)
275   - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
276   - std::cout << "*** ERROR: GCTAModelRadialPolynom::eval";
277   - std::cout << "(offset=" << offset << "): NaN/Inf encountered";
278   - std::cout << " (value=" << value;
279   - for (int i = 0; i < ncoeffs; ++i) {
280   - std::cout << ", c" << i << "=" << m_coeffs[i].value();
281   - }
282   - std::cout << ")" << std::endl;
283   - }
284   - #endif
285   -
286   - // Return value
287   - return value;
288   -}
289   -
290   -
291   -/***********************************************************************//**
292   - * @brief Evaluate function and gradients
293   - *
294   - * @param[in] offset Offset angle [degrees].
295 255 *
296   - * Evaluates the Polynom model for a given offset. The Polynom model is
297   - * defined as
298   - * Evaluates the Polynom model for a given offset. The Polynomian model is
299   - * defined as
300   - * \f[f(\theta) = \sum_{i=0}^m c_i \theta^i\f]
301   - * where
302   - * \f$\theta\f$ is the offset angle (in degrees), and
303   - * \f$c_i\f$ are the polynomial coefficients.
304   - *
305   - * The partial derivative of the Polynom are given by
  256 + * If the @p gradients flag is true the method will also compute the partial
  257 + * derivatives of the parameters. The partial derivative of the Polynom are
  258 + * given by
306 259 * \f[\frac{df}{d{c_i}_v} = {c_i}_s \theta^i\f]
307 260 * where
308 261 * \f${c_i}_v\f$ is the value part,
309 262 * \f${c_i}_s\f$ is the scaling part, and
310 263 * \f$c_i = {c_i}_v {c_i}_s\f$.
311   - *
312   - * Note that this method implements a function which is unity for
313   - * \f$\theta=0\f$.
314 264 ***************************************************************************/
315   -double GCTAModelRadialPolynom::eval_gradients(const double& offset) const
  265 +double GCTAModelRadialPolynom::eval(const double& offset,
  266 + const bool& gradients) const
316 267 {
317 268 // Initialise result
318 269 double value = 0.0;
... ... @@ -329,33 +280,38 @@ double GCTAModelRadialPolynom::eval_gradients(const double&amp; offset) const
329 280 value = value * offset + m_coeffs[i].value();
330 281 }
331 282  
332   - // Initialise theta^i for the first coefficient
333   - double offset_power = 1.0;
  283 + // Optionally compute partial derivatives
  284 + if (gradients) {
334 285  
335   - // Compute gradients for all coefficients
336   - for (int i = 0; i < ncoeffs; ++i) {
  286 + // Initialise theta^i for the first coefficient
  287 + double offset_power = 1.0;
  288 +
  289 + // Compute gradients for all coefficients
  290 + for (int i = 0; i < ncoeffs; ++i) {
337 291  
338   - // Compute gradient
339   - double grad = offset_power * m_coeffs[i].scale();
  292 + // Compute gradient
  293 + double grad = offset_power * m_coeffs[i].scale();
340 294  
341   - // Store gradient
342   - const_cast<GCTAModelRadialPolynom*>(this)->m_coeffs[i].factor_gradient(grad);
  295 + // Store gradient
  296 + m_coeffs[i].factor_gradient(grad);
343 297  
344   - // Increase offset power for next coefficient
345   - offset_power *= offset;
  298 + // Increase offset power for next coefficient
  299 + offset_power *= offset;
346 300  
347   - } // endfor: looped over all coefficients
  301 + } // endfor: looped over all coefficients
  302 +
  303 + } // endif: computed partial derivatives
348 304  
349 305 } // endif: there were coefficients
350 306  
351 307 // Compile option: Check for NaN/Inf
352 308 #if defined(G_NAN_CHECK)
353 309 if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
354   - std::cout << "*** ERROR: GCTAModelRadialPolynom::eval_gradients";
  310 + std::cout << "*** ERROR: GCTAModelRadialPolynom::eval";
355 311 std::cout << "(offset=" << offset << "): NaN/Inf encountered";
356 312 std::cout << " (value=" << value;
357 313 for (int i = 0; i < ncoeffs; ++i) {
358   - std::cout << ", c_" << i << "=" << m_coeffs[i].value();
  314 + std::cout << ", c" << i << "=" << m_coeffs[i].value();
359 315 }
360 316 std::cout << ")" << std::endl;
361 317 }
... ...
inst/cta/src/GCTAModelRadialProfile.cpp
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GCTAModelRadialProfile.cpp - Radial Profile CTA model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -224,6 +224,8 @@ GCTAModelRadialProfile* GCTAModelRadialProfile::clone(void) const
224 224 * @brief Evaluate function
225 225 *
226 226 * @param[in] offset Offset angle [degrees].
  227 + * @param[in] gradients Compute gradients?
  228 + * @return Function value
227 229 *
228 230 * Evaluates the Profile model for a given offset. The Profile model is
229 231 * defined as
... ... @@ -234,10 +236,11 @@ GCTAModelRadialProfile* GCTAModelRadialProfile::clone(void) const
234 236 * \f$c_1\f$ is the width of the central plateau (core), and
235 237 * \f$c_2\f$ is the size of the tail (tail).
236 238 *
237   - * Note that this method implements a function which is unity for
238   - * \f$\theta=0\f$.
  239 + * Note that no analytical partial derivatives are implemented for this
  240 + * function.
239 241 ***************************************************************************/
240   -double GCTAModelRadialProfile::eval(const double& offset) const
  242 +double GCTAModelRadialProfile::eval(const double& offset,
  243 + const bool& gradients) const
241 244 {
242 245 // Compute value
243 246 double arg = 1.0 + std::pow(offset / width(), core());
... ... @@ -262,22 +265,6 @@ double GCTAModelRadialProfile::eval(const double&amp; offset) const
262 265  
263 266  
264 267 /***********************************************************************//**
265   - * @brief Evaluate function and gradients
266   - *
267   - * @param[in] offset Offset angle [degrees].
268   - *
269   - * This method simply calls GCTAModelRadialProfile::eval() as no analytical
270   - * gradients will be computed. See GCTAModelRadialProfile::eval() for details
271   - * about the implemented method.
272   - ***************************************************************************/
273   -double GCTAModelRadialProfile::eval_gradients(const double& offset) const
274   -{
275   - // Return value
276   - return (eval(offset));
277   -}
278   -
279   -
280   -/***********************************************************************//**
281 268 * @brief Returns MC instrument direction
282 269 *
283 270 * @param[in] dir Pointing direction.
... ...
inst/cta/src/GCTAOnOffObservation.cpp
  View file @c135e7c
... ... @@ -1390,8 +1390,10 @@ double GCTAOnOffObservation::N_gamma(const GModels&amp; models,
1390 1390 // and the livetime (s)
1391 1391 value += spectral->flux(emin, emax) * norm_flux;
1392 1392  
1393   - // Determine the model gradients at the current energy
1394   - spectral->eval_gradients(emean);
  1393 + // Determine the model gradients at the current energy. The
  1394 + // eval() method needs a time in case that the spectral model
  1395 + // has a time dependence. We simply use a dummy time here.
  1396 + spectral->eval(emean, GTime(), true);
1395 1397  
1396 1398 // Loop over spectral model parameters
1397 1399 for (int k = 0; k < spectral->size(); ++k, ++ipar) {
... ... @@ -1445,8 +1447,8 @@ double GCTAOnOffObservation::N_gamma(const GModels&amp; models,
1445 1447 * background models; that's not very satisfactory ...
1446 1448 ***********************************************************************/
1447 1449 double GCTAOnOffObservation::N_bgd(const GModels& models,
1448   - const int& ibin,
1449   - GVector* grad) const
  1450 + const int& ibin,
  1451 + GVector* grad) const
1450 1452 {
1451 1453 // Get total number of model parameters
1452 1454 int npars = models.npars();
... ... @@ -1505,8 +1507,10 @@ double GCTAOnOffObservation::N_bgd(const GModels&amp; models,
1505 1507 // Determine the number of background events in model by
1506 1508 // computing the model normalization at the mean bin energy bin
1507 1509 // and multiplying the normalisation with the number of
1508   - // background events
1509   - value += spectral->eval_gradients(emean) * norm;
  1510 + // background events. The eval() method needs a time in case
  1511 + // that the spectral model has a time dependence. We simply
  1512 + // use a dummy time here.
  1513 + value += spectral->eval(emean, GTime(), true) * norm;
1510 1514  
1511 1515 // Compute the parameter gradients for all spectral model
1512 1516 // parameters
... ...
inst/cta/test/dev/test_gauss.py
  View file @c135e7c
... ... @@ -39,9 +39,9 @@ def gauss(dir, sigma):
39 39 t = GTime()
40 40 f1 = GModelSpatialRadialGauss(s, sigma)
41 41 f2 = GModelSpatialRadialGauss(s, sigma + dh)
42   - v1 = f1.eval_gradients(GPhoton(dir, e, t))
  42 + v1 = f1.eval(GPhoton(dir, e, t), True)
43 43 g1 = f1[2].gradient()
44   - v2 = f2.eval_gradients(GPhoton(dir, e, t))
  44 + v2 = f2.eval(GPhoton(dir, e, t), True)
45 45 g2 = f2[2].gradient()
46 46 g = (v2 - v1) / dh
47 47  
... ...
inst/cta/test/dev/test_radial_acceptance.py
  View file @c135e7c
... ... @@ -39,7 +39,7 @@ def show_model(model):
39 39  
40 40 # Get model values
41 41 values = [model.eval(theta) for theta in thetas]
42   - values_grad = [model.eval_gradients(theta) for theta in thetas]
  42 + values_grad = [model.eval(theta, True) for theta in thetas]
43 43  
44 44 # Create figure
45 45 plt.figure(1)
... ...
pyext/GModel.i
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GModel.i - Abstract virtual model base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -51,9 +51,8 @@ public:
51 51 virtual std::string type(void) const = 0;
52 52 virtual bool is_constant(void) const = 0;
53 53 virtual double eval(const GEvent& event,
54   - const GObservation& obs) const = 0;
55   - virtual double eval_gradients(const GEvent& event,
56   - const GObservation& obs) const = 0;
  54 + const GObservation& obs,
  55 + const bool& gradients = false) const = 0;
57 56 virtual double npred(const GEnergy& obsEng, const GTime& obsTime,
58 57 const GObservation& obs) const = 0;
59 58 virtual void read(const GXmlElement& xml) = 0;
... ...
pyext/GModelData.i
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GModelData.i - Abstract virtual data model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -50,9 +50,8 @@ public:
50 50 virtual std::string type(void) const = 0;
51 51 virtual bool is_constant(void) const = 0;
52 52 virtual double eval(const GEvent& event,
53   - const GObservation& obs) const = 0;
54   - virtual double eval_gradients(const GEvent& event,
55   - const GObservation& obs) const = 0;
  53 + const GObservation& obs,
  54 + const bool& gradients = false) const = 0;
56 55 virtual double npred(const GEnergy& obsEng, const GTime& obsTime,
57 56 const GObservation& obs) const = 0;
58 57 virtual GEvents* mc(const GObservation& obs, GRan& ran) const = 0;
... ...
pyext/GModelSky.i
  View file @c135e7c
... ... @@ -61,9 +61,8 @@ public:
61 61 virtual std::string type(void) const;
62 62 virtual bool is_constant(void) const;
63 63 virtual double eval(const GEvent& event,
64   - const GObservation& obs) const;
65   - virtual double eval_gradients(const GEvent& event,
66   - const GObservation& obs) const;
  64 + const GObservation& obs,
  65 + const bool& gradients = false) const;
67 66 virtual double npred(const GEnergy& obsEng,
68 67 const GTime& obsTime,
69 68 const GObservation& obs) const;
... ...
pyext/GModelSpatial.i
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GModelSpatial.i - Spatial model abstract base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -47,8 +47,8 @@ public:
47 47 virtual GModelSpatial* clone(void) const = 0;
48 48 virtual std::string classname(void) const = 0;
49 49 virtual std::string type(void) const = 0;
50   - virtual double eval(const GPhoton& photon) const = 0;
51   - virtual double eval_gradients(const GPhoton& photon) const = 0;
  50 + virtual double eval(const GPhoton& photon,
  51 + const bool& gradients = false) const = 0;
52 52 virtual GSkyDir mc(const GEnergy& energy, const GTime& time,
53 53 GRan& ran) const = 0;
54 54 virtual double mc_norm(const GSkyDir& dir,
... ...
pyext/GModelSpatialDiffuse.i
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GModelSpatialDiffuse.i - Abstract diffuse spatial model base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2013-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2013-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -46,8 +46,8 @@ public:
46 46 virtual GModelSpatialDiffuse* clone(void) const = 0;
47 47 virtual std::string classname(void) const = 0;
48 48 virtual std::string type(void) const = 0;
49   - virtual double eval(const GPhoton& photon) const = 0;
50   - virtual double eval_gradients(const GPhoton& photon) const = 0;
  49 + virtual double eval(const GPhoton& photon,
  50 + const bool& gradients = false) const = 0;
51 51 virtual GSkyDir mc(const GEnergy& energy, const GTime& time,
52 52 GRan& ran) const = 0;
53 53 virtual double mc_norm(const GSkyDir& dir,
... ...
pyext/GModelSpatialDiffuseConst.i
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GModelSpatialDiffuseConst.i - Spatial isotropic model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -48,8 +48,8 @@ public:
48 48 virtual GModelSpatialDiffuseConst* clone(void) const;
49 49 virtual std::string classname(void) const;
50 50 virtual std::string type(void) const;
51   - virtual double eval(const GPhoton& photon) const;
52   - virtual double eval_gradients(const GPhoton& photon) const;
  51 + virtual double eval(const GPhoton& photon,
  52 + const bool& gradients = false) const;
53 53 virtual GSkyDir mc(const GEnergy& energy,
54 54 const GTime& time,
55 55 GRan& ran) const;
... ...
pyext/GModelSpatialDiffuseCube.i
  View file @c135e7c
... ... @@ -53,11 +53,11 @@ public:
53 53 virtual GModelSpatialDiffuseCube* clone(void) const;
54 54 virtual std::string classname(void) const;
55 55 virtual std::string type(void) const;
56   - virtual double eval(const GPhoton& photon) const;
57   - virtual double eval_gradients(const GPhoton& photon) const;
  56 + virtual double eval(const GPhoton& photon,
  57 + const bool& gradients = false) const;
58 58 virtual GSkyDir mc(const GEnergy& energy,
59   - const GTime& time,
60   - GRan& ran) const;
  59 + const GTime& time,
  60 + GRan& ran) const;
61 61 virtual double mc_norm(const GSkyDir& dir,
62 62 const double& radius) const;
63 63 virtual bool contains(const GSkyDir& dir,
... ...
pyext/GModelSpatialDiffuseMap.i
  View file @c135e7c
... ... @@ -53,8 +53,8 @@ public:
53 53 virtual GModelSpatialDiffuseMap* clone(void) const;
54 54 virtual std::string classname(void) const;
55 55 virtual std::string type(void) const;
56   - virtual double eval(const GPhoton& photon) const;
57   - virtual double eval_gradients(const GPhoton& photon) const;
  56 + virtual double eval(const GPhoton& photon,
  57 + const bool& gradients = false) const;
58 58 virtual GSkyDir mc(const GEnergy& energy,
59 59 const GTime& time,
60 60 GRan& ran) const;
... ...
pyext/GModelSpatialElliptical.i
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GModelSpatialElliptical.i - Abstract elliptical spatial model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2013-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2013-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -50,21 +50,18 @@ public:
50 50 virtual double eval(const double& theta,
51 51 const double& posangle,
52 52 const GEnergy& energy,
53   - const GTime& time) const = 0;
54   - virtual double eval_gradients(const double& theta,
55   - const double& posangle,
56   - const GEnergy& energy,
57   - const GTime& time) const = 0;
  53 + const GTime& time,
  54 + const bool& gradients = false) const = 0;
58 55 virtual bool contains(const GSkyDir& dir,
59 56 const double& margin = 0.0) const = 0;
60 57 virtual GSkyDir mc(const GEnergy& energy,
61   - const GTime& time,
62   - GRan& ran) const = 0;
  58 + const GTime& time,
  59 + GRan& ran) const = 0;
63 60 virtual double theta_max(void) const = 0;
64 61  
65 62 // Implemented virtual methods
66   - virtual double eval(const GPhoton& photon) const;
67   - virtual double eval_gradients(const GPhoton& photon) const;
  63 + virtual double eval(const GPhoton& photon,
  64 + const bool& gradients = false) const;
68 65 virtual double mc_norm(const GSkyDir& dir, const double& radius) const;
69 66 virtual void read(const GXmlElement& xml);
70 67 virtual void write(GXmlElement& xml) const;
... ...
pyext/GModelSpatialEllipticalDisk.i
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GModelSpatialEllipticalDisk.i - Elliptical disk source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2013-2015 by Michael Mayer *
  4 + * copyright (C) 2013-2016 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -55,11 +55,8 @@ public:
55 55 virtual double eval(const double& theta,
56 56 const double& posangle,
57 57 const GEnergy& energy,
58   - const GTime& time) const;
59   - virtual double eval_gradients(const double& theta,
60   - const double& posangle,
61   - const GEnergy& energy,
62   - const GTime& time) const;
  58 + const GTime& time,
  59 + const bool& gradients = false) const;
63 60 virtual GSkyDir mc(const GEnergy& energy,
64 61 const GTime& time,
65 62 GRan& ran) const;
... ... @@ -78,10 +75,7 @@ public:
78 75 GModelSpatialEllipticalDisk copy() {
79 76 return (*self);
80 77 }
81   - double eval(const GPhoton& photon) const {
82   - return self->GModelSpatialElliptical::eval(photon);
83   - }
84   - double eval_gradients(const GPhoton& photon) const {
85   - return self->GModelSpatialElliptical::eval_gradients(photon);
  78 + double eval(const GPhoton& photon, const bool& gradients) const {
  79 + return self->GModelSpatialElliptical::eval(photon, gradients);
86 80 }
87 81 };
... ...
pyext/GModelSpatialEllipticalGauss.i
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GModelSpatialEllipticalGauss.i - Elliptical gauss source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2015 by Michael Mayer *
  4 + * copyright (C) 2015-2016 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -55,11 +55,8 @@ public:
55 55 virtual double eval(const double& theta,
56 56 const double& posangle,
57 57 const GEnergy& energy,
58   - const GTime& time) const;
59   - virtual double eval_gradients(const double& theta,
60   - const double& posangle,
61   - const GEnergy& energy,
62   - const GTime& time) const;
  58 + const GTime& time,
  59 + const bool& gradients = false) const;
63 60 virtual GSkyDir mc(const GEnergy& energy,
64 61 const GTime& time,
65 62 GRan& ran) const;
... ... @@ -78,10 +75,7 @@ public:
78 75 GModelSpatialEllipticalGauss copy() {
79 76 return (*self);
80 77 }
81   - double eval(const GPhoton& photon) const {
82   - return self->GModelSpatialElliptical::eval(photon);
83   - }
84   - double eval_gradients(const GPhoton& photon) const {
85   - return self->GModelSpatialElliptical::eval_gradients(photon);
  78 + double eval(const GPhoton& photon, const bool& gradients) const {
  79 + return self->GModelSpatialElliptical::eval(photon, gradients);
86 80 }
87 81 };
... ...
pyext/GModelSpatialPointSource.i
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GModelSpatialPointSource.i - Spatial point source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -50,8 +50,8 @@ public:
50 50 virtual GModelSpatialPointSource* clone(void) const;
51 51 virtual std::string classname(void) const;
52 52 virtual std::string type(void) const;
53   - virtual double eval(const GPhoton& photon) const;
54   - virtual double eval_gradients(const GPhoton& photon) const;
  53 + virtual double eval(const GPhoton& photon,
  54 + const bool& gradients = false) const;
55 55 virtual GSkyDir mc(const GEnergy& energy,
56 56 const GTime& time,
57 57 GRan& ran) const;
... ...
pyext/GModelSpatialRadial.i
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GModelSpatialRadial.i - Abstract radial spatial model base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -50,10 +50,8 @@ public:
50 50 virtual std::string type(void) const = 0;
51 51 virtual double eval(const double& theta,
52 52 const GEnergy& energy,
53   - const GTime& time) const = 0;
54   - virtual double eval_gradients(const double& theta,
55   - const GEnergy& energy,
56   - const GTime& time) const = 0;
  53 + const GTime& time,
  54 + const bool& gradients = false) const = 0;
57 55 virtual GSkyDir mc(const GEnergy& energy,
58 56 const GTime& time,
59 57 GRan& ran) const = 0;
... ... @@ -62,8 +60,8 @@ public:
62 60 virtual double theta_max(void) const = 0;
63 61  
64 62 // Implemented virtual base class methods
65   - virtual double eval(const GPhoton& photon) const;
66   - virtual double eval_gradients(const GPhoton& photon) const;
  63 + virtual double eval(const GPhoton& photon,
  64 + const bool& gradients = false) const;
67 65 virtual double mc_norm(const GSkyDir& dir, const double& radius) const;
68 66 virtual void read(const GXmlElement& xml);
69 67 virtual void write(GXmlElement& xml) const;
... ...
pyext/GModelSpatialRadialDisk.i
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GModelSpatialRadialDisk.i - Radial disk source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2015 by Christoph Deil *
  4 + * copyright (C) 2011-2016 by Christoph Deil *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -54,10 +54,8 @@ public:
54 54 virtual std::string type(void) const;
55 55 virtual double eval(const double& theta,
56 56 const GEnergy& energy,
57   - const GTime& time) const;
58   - virtual double eval_gradients(const double& theta,
59   - const GEnergy& energy,
60   - const GTime& time) const;
  57 + const GTime& time,
  58 + const bool& gradients = false) const;
61 59 virtual GSkyDir mc(const GEnergy& energy,
62 60 const GTime& time,
63 61 GRan& ran) const;
... ... @@ -76,9 +74,9 @@ public:
76 74 /***********************************************************************//**
77 75 * @brief GModelSpatialRadialDisk class extension
78 76 *
79   - * The eval(GSkyDir&) and eval_gradients(GSkyDir&) need to be defined in the
80   - * extension to force swig to build also the interface for these methods that
81   - * are implemented in the base class only. It's not clear to me why these
  77 + * The eval(GPhoton&, bool&) method need to be defined in the extension to
  78 + * force swig to build also the interface for these methods that are
  79 + * implemented in the base class only. It's not clear to me why these
82 80 * methods are not inherited automatically. Maybe this could also be handled
83 81 * by a %typemap(typecheck) construct.
84 82 ***************************************************************************/
... ... @@ -86,10 +84,7 @@ public:
86 84 GModelSpatialRadialDisk copy() {
87 85 return (*self);
88 86 }
89   - double eval(const GPhoton& photon) const {
90   - return self->GModelSpatialRadial::eval(photon);
91   - }
92   - double eval_gradients(const GPhoton& photon) const {
93   - return self->GModelSpatialRadial::eval_gradients(photon);
  87 + double eval(const GPhoton& photon, const bool& gradients) const {
  88 + return self->GModelSpatialRadial::eval(photon, gradients);
94 89 }
95 90 };
... ...
pyext/GModelSpatialRadialGauss.i
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GModelSpatialRadialGauss.i - Radial Gaussian source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -51,10 +51,8 @@ public:
51 51 virtual std::string type(void) const;
52 52 virtual double eval(const double& theta,
53 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;
  54 + const GTime& time,
  55 + const bool& gradients = false) const;
58 56 virtual GSkyDir mc(const GEnergy& energy,
59 57 const GTime& time,
60 58 GRan& ran) const;
... ... @@ -73,9 +71,9 @@ public:
73 71 /***********************************************************************//**
74 72 * @brief GModelSpatialRadialGauss class extension
75 73 *
76   - * The eval(GSkyDir&) and eval_gradients(GSkyDir&) need to be defined in the
77   - * extension to force swig to build also the interface for these methods that
78   - * are implemented in the base class only. It's not clear to me why these
  74 + * The eval(GPhoton&, bool&) method need to be defined in the extension to
  75 + * force swig to build also the interface for these methods that are
  76 + * implemented in the base class only. It's not clear to me why these
79 77 * methods are not inherited automatically. Maybe this could also be handled
80 78 * by a %typemap(typecheck) construct.
81 79 ***************************************************************************/
... ... @@ -83,10 +81,7 @@ public:
83 81 GModelSpatialRadialGauss copy() {
84 82 return (*self);
85 83 }
86   - double eval(const GPhoton& photon) const {
87   - return self->GModelSpatialRadial::eval(photon);
88   - }
89   - double eval_gradients(const GPhoton& photon) const {
90   - return self->GModelSpatialRadial::eval_gradients(photon);
  84 + double eval(const GPhoton& photon, const bool& gradients) const {
  85 + return self->GModelSpatialRadial::eval(photon, gradients);
91 86 }
92 87 };
... ...
pyext/GModelSpatialRadialShell.i
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GModelSpatialRadialShell.i - Radial spatial shell source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2015 by Christoph Deil *
  4 + * copyright (C) 2011-2016 by Christoph Deil *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -52,10 +52,8 @@ public:
52 52 virtual std::string type(void) const;
53 53 virtual double eval(const double& theta,
54 54 const GEnergy& energy,
55   - const GTime& time) const;
56   - virtual double eval_gradients(const double& theta,
57   - const GEnergy& energy,
58   - const GTime& time) const;
  55 + const GTime& time,
  56 + const bool& gradients = false) const;
59 57 virtual GSkyDir mc(const GEnergy& energy,
60 58 const GTime& time,
61 59 GRan& ran) const;
... ... @@ -76,9 +74,9 @@ public:
76 74 /***********************************************************************//**
77 75 * @brief GModelSpatialRadialShell class extension
78 76 *
79   - * The eval(GSkyDir&) and eval_gradients(GSkyDir&) need to be defined in the
80   - * extension to force swig to build also the interface for these methods that
81   - * are implemented in the base class only. It's not clear to me why these
  77 + * The eval(GPhoton&, bool&) method need to be defined in the extension to
  78 + * force swig to build also the interface for these methods that are
  79 + * implemented in the base class only. It's not clear to me why these
82 80 * methods are not inherited automatically. Maybe this could also be handled
83 81 * by a %typemap(typecheck) construct.
84 82 ***************************************************************************/
... ... @@ -86,10 +84,7 @@ public:
86 84 GModelSpatialRadialShell copy() {
87 85 return (*self);
88 86 }
89   - double eval(const GPhoton& photon) const {
90   - return self->GModelSpatialRadial::eval(photon);
91   - }
92   - double eval_gradients(const GPhoton& photon) const {
93   - return self->GModelSpatialRadial::eval_gradients(photon);
  87 + double eval(const GPhoton& photon, const bool& gradients) const {
  88 + return self->GModelSpatialRadial::eval(photon, gradients);
94 89 }
95 90 };
... ...
pyext/GModelSpectral.i
  View file @c135e7c
... ... @@ -48,9 +48,8 @@ public:
48 48 virtual std::string classname(void) const = 0;
49 49 virtual std::string type(void) const = 0;
50 50 virtual double eval(const GEnergy& srcEng,
51   - const GTime& srcTime = GTime()) const = 0;
52   - virtual double eval_gradients(const GEnergy& srcEng,
53   - const GTime& srcTime = GTime()) = 0;
  51 + const GTime& srcTime = GTime(),
  52 + const bool& gradients = false) const = 0;
54 53 virtual double flux(const GEnergy& emin,
55 54 const GEnergy& emax) const = 0;
56 55 virtual double eflux(const GEnergy& emin,
... ...
pyext/GModelSpectralBrokenPlaw.i
  View file @c135e7c
... ... @@ -57,9 +57,8 @@ public:
57 57 virtual std::string classname(void) const;
58 58 virtual std::string type(void) const;
59 59 virtual double eval(const GEnergy& srcEng,
60   - const GTime& srcTime = GTime()) const;
61   - virtual double eval_gradients(const GEnergy& srcEng,
62   - const GTime& srcTime = GTime());
  60 + const GTime& srcTime = GTime(),
  61 + const bool& gradients = false) const;
63 62 virtual double flux(const GEnergy& emin,
64 63 const GEnergy& emax) const;
65 64 virtual double eflux(const GEnergy& emin,
... ...
pyext/GModelSpectralConst.i
  View file @c135e7c
... ... @@ -50,9 +50,8 @@ public:
50 50 virtual std::string classname(void) const;
51 51 virtual std::string type(void) const;
52 52 virtual double eval(const GEnergy& srcEng,
53   - const GTime& srcTime = GTime()) const;
54   - virtual double eval_gradients(const GEnergy& srcEng,
55   - const GTime& srcTime = GTime());
  53 + const GTime& srcTime = GTime(),
  54 + const bool& gradients = false) const;
56 55 virtual double flux(const GEnergy& emin,
57 56 const GEnergy& emax) const;
58 57 virtual double eflux(const GEnergy& emin,
... ...
pyext/GModelSpectralExpInvPlaw.i
  View file @c135e7c
... ... @@ -56,18 +56,14 @@ public:
56 56 GModelSpectralExpInvPlaw(const GModelSpectralExpInvPlaw& model);
57 57 virtual ~GModelSpectralExpInvPlaw(void);
58 58  
59   - // Operators
60   - virtual GModelSpectralExpInvPlaw& operator=(const GModelSpectralExpInvPlaw& model);
61   -
62 59 // Implemented pure virtual methods
63 60 virtual void clear(void);
64 61 virtual GModelSpectralExpInvPlaw* clone(void) const;
65 62 virtual std::string classname(void) const;
66 63 virtual std::string type(void) const;
67 64 virtual double eval(const GEnergy& srcEng,
68   - const GTime& srcTime = GTime()) const;
69   - virtual double eval_gradients(const GEnergy& srcEng,
70   - const GTime& srcTime = GTime());
  65 + const GTime& srcTime = GTime(),
  66 + const bool& gradients = false) const;
71 67 virtual double flux(const GEnergy& emin,
72 68 const GEnergy& emax) const;
73 69 virtual double eflux(const GEnergy& emin,
... ...
pyext/GModelSpectralExpPlaw.i
  View file @c135e7c
... ... @@ -57,9 +57,8 @@ public:
57 57 virtual std::string classname(void) const;
58 58 virtual std::string type(void) const;
59 59 virtual double eval(const GEnergy& srcEng,
60   - const GTime& srcTime = GTime()) const;
61   - virtual double eval_gradients(const GEnergy& srcEng,
62   - const GTime& srcTime = GTime());
  60 + const GTime& srcTime = GTime(),
  61 + const bool& gradients = false) const;
63 62 virtual double flux(const GEnergy& emin,
64 63 const GEnergy& emax) const;
65 64 virtual double eflux(const GEnergy& emin,
... ...
pyext/GModelSpectralFunc.i
  View file @c135e7c
... ... @@ -50,9 +50,8 @@ public:
50 50 virtual std::string classname(void) const;
51 51 virtual std::string type(void) const;
52 52 virtual double eval(const GEnergy& srcEng,
53   - const GTime& srcTime = GTime()) const;
54   - virtual double eval_gradients(const GEnergy& srcEng,
55   - const GTime& srcTime = GTime());
  53 + const GTime& srcTime = GTime(),
  54 + const bool& gradients = false) const;
56 55 virtual double flux(const GEnergy& emin,
57 56 const GEnergy& emax) const;
58 57 virtual double eflux(const GEnergy& emin,
... ...
pyext/GModelSpectralGauss.i
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GModelSpectralGauss.i - Spectral gaussian model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014-2015 by Christoph Deil & Ellis Owen *
  4 + * copyright (C) 2014-2016 by Christoph Deil & Ellis Owen *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -51,9 +51,8 @@ public:
51 51 virtual std::string classname(void) const;
52 52 virtual std::string type(void) const;
53 53 virtual double eval(const GEnergy& srcEng,
54   - const GTime& srcTime = GTime()) const;
55   - virtual double eval_gradients(const GEnergy& srcEng,
56   - const GTime& srcTime = GTime());
  54 + const GTime& srcTime = GTime(),
  55 + const bool& gradients = false) const;
57 56 virtual double flux(const GEnergy& emin,
58 57 const GEnergy& emax) const;
59 58 virtual double eflux(const GEnergy& emin,
... ...
pyext/GModelSpectralLogParabola.i
  View file @c135e7c
... ... @@ -57,9 +57,8 @@ public:
57 57 virtual std::string classname(void) const;
58 58 virtual std::string type(void) const;
59 59 virtual double eval(const GEnergy& srcEng,
60   - const GTime& srcTime = GTime()) const;
61   - virtual double eval_gradients(const GEnergy& srcEng,
62   - const GTime& srcTime = GTime());
  60 + const GTime& srcTime = GTime(),
  61 + const bool& gradients = false) const;
63 62 virtual double flux(const GEnergy& emin,
64 63 const GEnergy& emax) const;
65 64 virtual double eflux(const GEnergy& emin,
... ...
pyext/GModelSpectralNodes.i
  View file @c135e7c
... ... @@ -48,9 +48,8 @@ public:
48 48 virtual std::string classname(void) const;
49 49 virtual std::string type(void) const;
50 50 virtual double eval(const GEnergy& srcEng,
51   - const GTime& srcTime = GTime()) const;
52   - virtual double eval_gradients(const GEnergy& srcEng,
53   - const GTime& srcTime = GTime());
  51 + const GTime& srcTime = GTime(),
  52 + const bool& gradients = false) const;
54 53 virtual double flux(const GEnergy& emin,
55 54 const GEnergy& emax) const;
56 55 virtual double eflux(const GEnergy& emin,
... ...
pyext/GModelSpectralPlaw.i
  View file @c135e7c
... ... @@ -55,9 +55,8 @@ public:
55 55 virtual std::string classname(void) const;
56 56 virtual std::string type(void) const;
57 57 virtual double eval(const GEnergy& srcEng,
58   - const GTime& srcTime = GTime()) const;
59   - virtual double eval_gradients(const GEnergy& srcEng,
60   - const GTime& srcTime = GTime());
  58 + const GTime& srcTime = GTime(),
  59 + const bool& gradients = false) const;
61 60 virtual double flux(const GEnergy& emin,
62 61 const GEnergy& emax) const;
63 62 virtual double eflux(const GEnergy& emin,
... ...
pyext/GModelSpectralPlawEnergyFlux.i
  View file @c135e7c
... ... @@ -57,9 +57,8 @@ public:
57 57 virtual std::string classname(void) const;
58 58 virtual std::string type(void) const;
59 59 virtual double eval(const GEnergy& srcEng,
60   - const GTime& srcTime = GTime()) const;
61   - virtual double eval_gradients(const GEnergy& srcEng,
62   - const GTime& srcTime = GTime());
  60 + const GTime& srcTime = GTime(),
  61 + const bool& gradients = false) const;
63 62 virtual double flux(const GEnergy& emin,
64 63 const GEnergy& emax) const;
65 64 virtual double eflux(const GEnergy& emin,
... ...
pyext/GModelSpectralPlawPhotonFlux.i
  View file @c135e7c
... ... @@ -57,9 +57,8 @@ public:
57 57 virtual std::string classname(void) const;
58 58 virtual std::string type(void) const;
59 59 virtual double eval(const GEnergy& srcEng,
60   - const GTime& srcTime = GTime()) const;
61   - virtual double eval_gradients(const GEnergy& srcEng,
62   - const GTime& srcTime = GTime());
  60 + const GTime& srcTime = GTime(),
  61 + const bool& gradients = false) const;
63 62 virtual double flux(const GEnergy& emin,
64 63 const GEnergy& emax) const;
65 64 virtual double eflux(const GEnergy& emin,
... ...
pyext/GModelSpectralSuperExpPlaw.i
  View file @c135e7c
... ... @@ -59,9 +59,8 @@ public:
59 59 virtual std::string classname(void) const;
60 60 virtual std::string type(void) const;
61 61 virtual double eval(const GEnergy& srcEng,
62   - const GTime& srcTime = GTime()) const;
63   - virtual double eval_gradients(const GEnergy& srcEng,
64   - const GTime& srcTime = GTime());
  62 + const GTime& srcTime = GTime(),
  63 + const bool& gradients = false) const;
65 64 virtual double flux(const GEnergy& emin,
66 65 const GEnergy& emax) const;
67 66 virtual double eflux(const GEnergy& emin,
... ...
pyext/GModelTemporal.i
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GModelTemporal.i - Abstract temporal model base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -47,10 +47,12 @@ public:
47 47 virtual GModelTemporal* clone(void) const = 0;
48 48 virtual std::string classname(void) const = 0;
49 49 virtual std::string type(void) const = 0;
50   - virtual double eval(const GTime& srcTime) const = 0;
51   - virtual double eval_gradients(const GTime& srcTime) = 0;
52   - virtual GTimes mc(const double& rate, const GTime& tmin,
53   - const GTime& tmax, GRan& ran) const = 0;
  50 + virtual double eval(const GTime& srcTime,
  51 + const bool& gradients = false) const = 0;
  52 + virtual GTimes mc(const double& rate,
  53 + const GTime& tmin,
  54 + const GTime& tmax,
  55 + GRan& ran) const = 0;
54 56 virtual void read(const GXmlElement& xml) = 0;
55 57 virtual void write(GXmlElement& xml) const = 0;
56 58  
... ...
pyext/GModelTemporalConst.i
  View file @c135e7c
1 1 /***************************************************************************
2 2 * GModelTemporalConst.i - Temporal constant model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -48,10 +48,12 @@ public:
48 48 virtual GModelTemporalConst* clone(void) const;
49 49 virtual std::string classname(void) const;
50 50 virtual std::string type(void) const;
51   - virtual double eval(const GTime& srcTime) const;
52   - virtual double eval_gradients(const GTime& srcTime);
53   - virtual GTimes mc(const double& rate, const GTime& tmin,
54   - const GTime& tmax, GRan& ran) const;
  51 + virtual double eval(const GTime& srcTime,
  52 + const bool& gradients = false) const;
  53 + virtual GTimes mc(const double& rate,
  54 + const GTime& tmin,
  55 + const GTime& tmax,
  56 + GRan& ran) const;
55 57 virtual void read(const GXmlElement& xml);
56 58 virtual void write(GXmlElement& xml) const;
57 59  
... ...
pyext/GModels.i
  View file @c135e7c
... ... @@ -68,8 +68,9 @@ public:
68 68 void write(GXml& xml) const;
69 69 int npars(void) const;
70 70 GOptimizerPars pars(void) const;
71   - double eval(const GEvent& event, const GObservation& obs) const;
72   - double eval_gradients(const GEvent& event, const GObservation& obs) const;
  71 + double eval(const GEvent& event,
  72 + const GObservation& obs,
  73 + const bool& gradients = false) const;
73 74 };
74 75  
75 76  
... ...
src/model/GModelSky.cpp
  View file @c135e7c
... ... @@ -516,10 +516,10 @@ double GModelSky::value(const GPhoton&amp; photon)
516 516 GVector GModelSky::gradients(const GPhoton& photon)
517 517 {
518 518 // Evaluate source model gradients
519   - if (m_spatial != NULL) m_spatial->eval_gradients(photon);
520   - if (m_spectral != NULL) m_spectral->eval_gradients(photon.energy(),
521   - photon.time());
522   - if (m_temporal != NULL) m_temporal->eval_gradients(photon.time());
  519 + if (m_spatial != NULL) m_spatial->eval(photon, true);
  520 + if (m_spectral != NULL) m_spectral->eval(photon.energy(), photon.time(),
  521 + true);
  522 + if (m_temporal != NULL) m_temporal->eval(photon.time(), true);
523 523  
524 524 // Set vector of gradients
525 525 GVector gradients;
... ... @@ -540,40 +540,21 @@ GVector GModelSky::gradients(const GPhoton&amp; photon)
540 540 *
541 541 * @param[in] event Observed event.
542 542 * @param[in] obs Observation.
543   - * @return Value of sky model
  543 + * @param[in] gradients Compute gradients?
  544 + * @return Value of sky model.
544 545 *
545 546 * Evalutes the value of the sky model for an @p event of a specific
546 547 * observation @p obs.
547   - ***************************************************************************/
548   -double GModelSky::eval(const GEvent& event, const GObservation& obs) const
549   -{
550   - // Evaluate function
551   - double value = obs.response()->convolve(*this, event, obs, false);
552   -
553   - // Return
554   - return value;
555   -}
556   -
557   -
558   -/***********************************************************************//**
559   - * @brief Evaluate sky model and parameter gradients for a given event of an
560   - * observation
561   - *
562   - * @param[in] event Observed event.
563   - * @param[in] obs Observation.
564   - * @return Value of sky model
565   - *
566   - * Evalutes the value of the sky model and of the parameter for an @p event
567   - * of a specific observation @p obs.
568 548 *
569   - * While the value of the sky model is returned by the method, the parameter
570   - * gradients are set as GModelPar members.
  549 + * If the @p gradients flag is true the method will also compute the
  550 + * parameter gradients for all model parameters.
571 551 ***************************************************************************/
572   -double GModelSky::eval_gradients(const GEvent& event,
573   - const GObservation& obs) const
  552 +double GModelSky::eval(const GEvent& event,
  553 + const GObservation& obs,
  554 + const bool& gradients) const
574 555 {
575 556 // Evaluate function
576   - double value = obs.response()->convolve(*this, event, obs, true);
  557 + double value = obs.response()->convolve(*this, event, obs, gradients);
577 558  
578 559 // Return
579 560 return value;
... ...
src/model/GModelSpatialDiffuseConst.cpp
  View file @c135e7c
... ... @@ -272,6 +272,7 @@ GModelSpatialDiffuseConst* GModelSpatialDiffuseConst::clone(void) const
272 272 * @brief Evaluate isotropic spatial model value
273 273 *
274 274 * @param[in] photon Incident photon.
  275 + * @param[in] gradients Compute gradients?
275 276 * @return Model value (\f$sr^{-1}\f$).
276 277 *
277 278 * Evaluates the isotropic spatial model for a given @p photon, characterised
... ... @@ -285,39 +286,9 @@ GModelSpatialDiffuseConst* GModelSpatialDiffuseConst::clone(void) const
285 286 *
286 287 * where \f$v\f$ is the value() parameter, which is divided by the solid
287 288 * angle \f$4\pi\f$ of the celestial sphere.
288   - ***************************************************************************/
289   -double GModelSpatialDiffuseConst::eval(const GPhoton& photon) const
290   -{
291   - // Set normalization constant
292   - const double norm = 1.0 / gammalib::fourpi;
293   -
294   - // Compute model value
295   - double value = norm * m_value.value();
296   -
297   - // Return model value
298   - return (value);
299   -}
300   -
301   -
302   -/***********************************************************************//**
303   - * @brief Evaluate isotropic spatial model value and parameter gradient
304 289 *
305   - * @param[in] photon Incident photon.
306   - * @return Model value (\f$sr^{-1}\f$).
307   - *
308   - * Evaluates the isotropic spatial model for a given @p photon, characterised
309   - * by a sky direction, energy and arrival time \f$(\vec{p},E,t)\f$, and
310   - * computes the gradient of the value() parameter. By definition, the
311   - * isotropic spatial model and gradient are independent of sky direction,
312   - * energy and time. The model value is given by
313   - *
314   - * \f[
315   - * M_\mathrm{S}(\vec{p}|E,t) = \frac{v}{4\pi}
316   - * \f]
317   - *
318   - * where \f$v\f$ is the value() parameter, which is divided by the solid
319   - * angle \f$4\pi\f$ of the celestial sphere. The value() gradient is given
320   - * by
  290 + * If the @p gradients flag is true the method will also evaluate the partial
  291 + * derivatives of the model. The value() gradient is given by
321 292 *
322 293 * \f[
323 294 * \frac{\partial M_\mathrm{S}(\vec{p}|E,t)}{\partial v_v} =
... ... @@ -327,7 +298,8 @@ double GModelSpatialDiffuseConst::eval(const GPhoton&amp; photon) const
327 298 * where \f$v=v_v v_s\f$ is the factorisation of the value() parameter into
328 299 * a value \f$v_v\f$ and scale \f$v_s\f$ term.
329 300 ***************************************************************************/
330   -double GModelSpatialDiffuseConst::eval_gradients(const GPhoton& photon) const
  301 +double GModelSpatialDiffuseConst::eval(const GPhoton& photon,
  302 + const bool& gradients) const
331 303 {
332 304 // Set normalization constant
333 305 const double norm = 1.0 / gammalib::fourpi;
... ... @@ -335,14 +307,19 @@ double GModelSpatialDiffuseConst::eval_gradients(const GPhoton&amp; photon) const
335 307 // Compute model value
336 308 double value = norm * m_value.value();
337 309  
338   - // Compute partial derivatives of the parameter values
339   - double g_norm = (m_value.is_free()) ? norm * m_value.scale() : 0.0;
  310 + // Optionally compute partial derivatives
  311 + if (gradients) {
  312 +
  313 + // Compute partial derivatives of the parameter values
  314 + double g_norm = (m_value.is_free()) ? norm * m_value.scale() : 0.0;
340 315  
341   - // Set gradient (circumvent const correctness)
342   - const_cast<GModelSpatialDiffuseConst*>(this)->m_value.factor_gradient(g_norm);
  316 + // Set gradient
  317 + m_value.factor_gradient(g_norm);
343 318  
344   - // Return value
345   - return value;
  319 + } // endif: computed partial derivatives
  320 +
  321 + // Return model value
  322 + return (value);
346 323 }
347 324  
348 325  
... ...