Simon Bonnefoy / gammalib

Commit 2c5c39b529a3d743f8deeecfc1e31d60fbe9f411

Authored by Michael Mayer
1 parent a447c53f

Fixed problem of wrong units in flux() and eflux() methods, now use MeV constently

Inline Side-by-side
Showing 2 changed files with 47 additions and 19 deletions
include/GModelSpectralMultiplicative.hpp
  View file @2c5c39b
... ... @@ -107,12 +107,12 @@ protected:
107 107  
108 108 // Method
109 109 double eval(const double& x) {
110   - GEnergy e(x, "TeV");
111   - double flux = m_spec[0]->eval(e);
  110 + GEnergy energy(x, "MeV");
  111 + double value = m_spec[0]->eval(energy);
112 112 for(int i=1;i<m_spec.size();++i){
113   - flux*=m_spec[i]->eval(e);
  113 + value*=m_spec[i]->eval(energy);
114 114 }
115   - return flux;
  115 + return value;
116 116 }
117 117  
118 118 // Members
... ...
src/model/GModelSpectralMultiplicative.cpp
  View file @2c5c39b
... ... @@ -41,8 +41,8 @@
41 41 /* __ Constants __________________________________________________________ */
42 42  
43 43 /* __ Globals ____________________________________________________________ */
44   -const GModelSpectralMultiplicative g_spectral_comp_seed;
45   -const GModelSpectralRegistry g_spectral_comp_registry(&g_spectral_comp_seed);
  44 +const GModelSpectralMultiplicative g_spectral_multi_seed;
  45 +const GModelSpectralRegistry g_spectral_multi_registry(&g_spectral_multi_seed);
46 46  
47 47 /* __ Method name definitions ____________________________________________ */
48 48 #define G_MC "GModelSpectralMultiplicative::mc(GEnergy&, GEnergy&, GTime&, GRan&)"
... ... @@ -81,7 +81,7 @@ GModelSpectralMultiplicative::GModelSpectralMultiplicative(void) : GModelSpectra
81 81 *
82 82 * @param[in] xml XML element containing position information.
83 83 *
84   - * Constructs a power law spectral model by extracting information from an
  84 + * Constructs a multiplicative spectral model by extracting information from an
85 85 * XML element. See the read() method for more information about the expected
86 86 * structure of the XML element.
87 87 ***************************************************************************/
... ... @@ -101,7 +101,7 @@ GModelSpectralMultiplicative::GModelSpectralMultiplicative(const GXmlElement&amp; xm
101 101 /***********************************************************************//**
102 102 * @brief Copy constructor
103 103 *
104   - * @param[in] model Spectral power law model.
  104 + * @param[in] model Multiplicative spectral model.
105 105 ***************************************************************************/
106 106 GModelSpectralMultiplicative::GModelSpectralMultiplicative(const GModelSpectralMultiplicative& model)
107 107 : GModelSpectral(model)
... ... @@ -173,7 +173,7 @@ GModelSpectralMultiplicative&amp; GModelSpectralMultiplicative::operator=(const GMod
173 173 ==========================================================================*/
174 174  
175 175 /***********************************************************************//**
176   - * @brief Clear spectral power law model
  176 + * @brief Clear multiplicative spectral model
177 177 ***************************************************************************/
178 178 void GModelSpectralMultiplicative::clear(void)
179 179 {
... ... @@ -191,9 +191,9 @@ void GModelSpectralMultiplicative::clear(void)
191 191  
192 192  
193 193 /***********************************************************************//**
194   - * @brief Clone spectral power law model
  194 + * @brief Clone multiplicative spectral model model
195 195 *
196   - * @return Pointer to deep copy of spectral power law model.
  196 + * @return Pointer to deep copy of multiplicative spectral model.
197 197 ***************************************************************************/
198 198 GModelSpectralMultiplicative* GModelSpectralMultiplicative::clone(void) const
199 199 {
... ... @@ -213,7 +213,7 @@ GModelSpectralMultiplicative* GModelSpectralMultiplicative::clone(void) const
213 213 * Evaluates
214 214 *
215 215 * \f[
216   - * \sum_{i=0}^{N} {M_{\rm i}}(\rm srcEng, srcTime)
  216 + * \prod_{i=0}^{N} {M_{\rm i}}(\rm srcEng, srcTime)
217 217 * \f]
218 218 *
219 219 * where \f${M_{\rm i}}\f$ is the model evaluation of the i-th model
... ... @@ -243,6 +243,36 @@ double GModelSpectralMultiplicative::eval(const GEnergy&amp; srcEng,
243 243  
244 244 } // endfor: loop over model components
245 245  
  246 + // Modify gradients if requested
  247 + if (gradients) {
  248 +
  249 + // Loop over models
  250 + for(int i = 0; i < m_spectral.size(); ++i) {
  251 +
  252 + // Initialise scaling factor
  253 + double factor = 1.0;
  254 +
  255 + //Loop over other models and compute factor
  256 + for(int j = 0; j < m_spectral.size(); ++j) {
  257 + if (i != j) {
  258 + factor *= m_spectral[j]->eval(srcEng, srcTime, false);
  259 + }
  260 + }
  261 +
  262 + // Loop over model parameters
  263 + for(int ipar = 0; ipar < m_spectral[i]->size(); ++ipar) {
  264 +
  265 + // Get reference to model parameter
  266 + GModelPar& par = m_spectral[i]->operator[](ipar);
  267 +
  268 + // Scale parameter gradient
  269 + par.gradient(par.gradient()*factor);
  270 +
  271 + } // endfor: loop over model parameters
  272 + } // endfor: loop over models
  273 + } //endif: compute gradients
  274 +
  275 +
246 276 // Compile option: Check for NaN/Inf
247 277 #if defined(G_NAN_CHECK)
248 278 if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
... ... @@ -267,8 +297,7 @@ double GModelSpectralMultiplicative::eval(const GEnergy&amp; srcEng,
267 297 * @param[in] emax Maximum photon energy.
268 298 * @return Photon flux (ph/cm2/s).
269 299 *
270   - * Computes the sum of all photon fluxes of each individual component
271   - *
  300 + * Computes the photon flux of multiplicative spectral model
272 301 ***************************************************************************/
273 302 double GModelSpectralMultiplicative::flux(const GEnergy& emin,
274 303 const GEnergy& emax) const
... ... @@ -289,7 +318,7 @@ double GModelSpectralMultiplicative::flux(const GEnergy&amp; emin,
289 318 integral.eps(1.0e-8);
290 319  
291 320 // Calculate integral between emin and emax
292   - flux = integral.romberg(emin.TeV(), emax.TeV());
  321 + flux = integral.romberg(emin.MeV(), emax.MeV());
293 322  
294 323 } // endif: integration range was valid
295 324  
... ... @@ -305,8 +334,7 @@ double GModelSpectralMultiplicative::flux(const GEnergy&amp; emin,
305 334 * @param[in] emax Maximum photon energy.
306 335 * @return Energy flux (erg/cm2/s).
307 336 *
308   - * Computes the sum of all energy fluxes of each individual component
309   - *
  337 + * Computes the energy flux of multiplicative spectral model
310 338 ***************************************************************************/
311 339 double GModelSpectralMultiplicative::eflux(const GEnergy& emin,
312 340 const GEnergy& emax) const
... ... @@ -327,7 +355,7 @@ double GModelSpectralMultiplicative::eflux(const GEnergy&amp; emin,
327 355 integral.eps(1.0e-8);
328 356  
329 357 // Calculate integral between emin and emax
330   - eflux = integral.romberg(emin.TeV(), emax.TeV());
  358 + eflux = integral.romberg(emin.MeV(), emax.MeV());
331 359  
332 360 } // endif: integration range was valid
333 361  
... ... @@ -348,7 +376,7 @@ double GModelSpectralMultiplicative::eflux(const GEnergy&amp; emin,
348 376 * @exception GException::erange_invalid
349 377 * Energy range is invalid (emin < emax required).
350 378 *
351   - * Returns Monte Carlo energy by randomly drawing from a power law.
  379 + * Returns Monte Carlo energy by randomly drawing from a multiplicative spectral model.
352 380 ***************************************************************************/
353 381 GEnergy GModelSpectralMultiplicative::mc(const GEnergy& emin,
354 382 const GEnergy& emax,
... ...