Justine Devin / gammalib

Commit 244c69d271ac4a7bb59fb92b48a555626990324d

Authored by Jürgen Knödlseder
1 parent 3a334a8c

Correct mc() and eval_gradients() methods and rename "Prefactor" to "Normalization". 

The mc() method went in endless loops as the exit condition of the while loop was inversed (the method continued looping when the energy was within the valid range and only exited when it was outside).

The eval_gradients() method computed the gradients with respect to the parameter values and not the parameter factors. All gradients have be multiplied by the scales. Also, the gradients are now set to 0 if the parameters are fixed.

Finally, the parameter "Prefactor" was renamed to "Normalization" because the parameter gives the integral flux under the Gaussian, not the prefactor.
Inline Side-by-side
Showing 5 changed files with 164 additions and 176 deletions
include/GModelSpectralGauss.hpp
  View file @244c69d
... ... @@ -45,12 +45,12 @@
45 45 *
46 46 * \f[
47 47 * S_{\rm E}(E | t) =
48   - * \frac{\tt m\_norm}{\sqrt{2 \pi} \sigma^2} \exp
  48 + * \frac{\tt m\_norm}{\sqrt{2 \pi} m\_sigma^2} \exp
49 49 * \left(\frac{-(E - m\_mean)^2}{2 m_\sigma^2} \right)
50 50 * \f]
51 51 *
52 52 * where
53   - * - \f${\tt m\_norm}\f$ is the normalization or prefactor,
  53 + * - \f${\tt m\_norm}\f$ is the normalization,
54 54 * - \f${\tt m\_mean}\f$ is the mean energy,
55 55 * - \f${\tt m\_sigma}\f$ is the energy width.
56 56 ***************************************************************************/
... ... @@ -60,7 +60,7 @@ public:
60 60 // Constructors and destructors
61 61 GModelSpectralGauss(void);
62 62 explicit GModelSpectralGauss(const GXmlElement& xml);
63   - GModelSpectralGauss(const double& prefactor,
  63 + GModelSpectralGauss(const double& norm,
64 64 const GEnergy& mean,
65 65 const GEnergy& sigma);
66 66 GModelSpectralGauss(const GModelSpectralGauss& model);
... ... @@ -90,14 +90,12 @@ public:
90 90 virtual std::string print(const GChatter& chatter = NORMAL) const;
91 91  
92 92 // Other methods
93   - double prefactor(void) const;
94   - void prefactor(const double& prefactor);
  93 + double norm(void) const;
  94 + void norm(const double& norm);
95 95 GEnergy mean(void) const;
96 96 void mean(const GEnergy& mean);
97 97 GEnergy sigma(void) const;
98 98 void sigma(const GEnergy& sigma);
99   - void update_eval_cache(const GEnergy& energy) const;
100   -
101 99  
102 100 protected:
103 101 // Protected methods
... ... @@ -126,10 +124,6 @@ protected:
126 124 GModelPar m_norm; //!< Normalization factor
127 125 GModelPar m_mean; //!< Gaussian mean energy
128 126 GModelPar m_sigma; //!< Gaussian energy width
129   -
130   - mutable double m_last_norm; //!< Last Normalization value
131   - mutable double m_last_mean; //!< Last Mean value
132   - mutable double m_last_sigma; //!< Last Sigma value
133 127 };
134 128  
135 129  
... ... @@ -148,26 +142,26 @@ std::string GModelSpectralGauss::type(void) const
148 142  
149 143  
150 144 /***********************************************************************//**
151   - * @brief Return pre factor
  145 + * @brief Return normalization
152 146 *
153   - * @return Pre factor (ph/cm2/s/MeV).
  147 + * @return Normalization (ph/cm2/s).
154 148 ***************************************************************************/
155 149 inline
156   -double GModelSpectralGauss::prefactor(void) const
  150 +double GModelSpectralGauss::norm(void) const
157 151 {
158 152 return (m_norm.value());
159 153 }
160 154  
161 155  
162 156 /***********************************************************************//**
163   - * @brief Set pre factor
  157 + * @brief Set normalization
164 158 *
165   - * @param[in] prefactor Pre factor (ph/cm2/s/MeV).
  159 + * @param[in] norm Normalization (ph/cm2/s).
166 160 ***************************************************************************/
167 161 inline
168   -void GModelSpectralGauss::prefactor(const double& prefactor)
  162 +void GModelSpectralGauss::norm(const double& norm)
169 163 {
170   - m_norm.value(prefactor);
  164 + m_norm.value(norm);
171 165 return;
172 166 }
173 167  
... ... @@ -223,5 +217,4 @@ void GModelSpectralGauss::sigma(const GEnergy&amp; sigma)
223 217 return;
224 218 }
225 219  
226   -
227 220 #endif /* GMODELSPECTRALGAUSS_HPP */
... ...
pyext/GModelSpectralGauss.i
  View file @244c69d
... ... @@ -41,7 +41,7 @@ public:
41 41 // Constructors and destructors
42 42 GModelSpectralGauss(void);
43 43 explicit GModelSpectralGauss(const GXmlElement& xml);
44   - GModelSpectralGauss(const double& prefactor,
  44 + GModelSpectralGauss(const double& norm,
45 45 const GEnergy& mean,
46 46 const GEnergy& sigma);
47 47 GModelSpectralGauss(const GModelSpectralGauss& model);
... ... @@ -67,8 +67,8 @@ public:
67 67 virtual void write(GXmlElement& xml) const;
68 68  
69 69 // Other methods
70   - double prefactor(void) const;
71   - void prefactor(const double& prefactor);
  70 + double norm(void) const;
  71 + void norm(const double& norm);
72 72 GEnergy mean(void) const;
73 73 void mean(const GEnergy& mean);
74 74 GEnergy sigma(void) const;
... ...
src/model/GModelSpectralGauss.cpp
  View file @244c69d
... ... @@ -101,11 +101,11 @@ GModelSpectralGauss::GModelSpectralGauss(const GXmlElement&amp; xml) :
101 101 /***********************************************************************//**
102 102 * @brief Constructor
103 103 *
104   - * @param[in] prefactor Power law pre factor (in ph/cm2/s/MeV).
  104 + * @param[in] norm Total flux under Gaussian (in ph/cm2/s).
105 105 * @param[in] mean Mean energy.
106 106 * @param[in] sigma Energy width.
107 107 ***************************************************************************/
108   -GModelSpectralGauss::GModelSpectralGauss(const double& prefactor,
  108 +GModelSpectralGauss::GModelSpectralGauss(const double& norm,
109 109 const GEnergy& mean,
110 110 const GEnergy& sigma) :
111 111 GModelSpectral()
... ... @@ -114,7 +114,7 @@ GModelSpectralGauss::GModelSpectralGauss(const double&amp; prefactor,
114 114 init_members();
115 115  
116 116 // Set parameters
117   - m_norm.value(prefactor);
  117 + m_norm.value(norm);
118 118 m_mean.value(mean.MeV());
119 119 m_sigma.value(sigma.MeV());
120 120  
... ... @@ -240,22 +240,42 @@ GModelSpectralGauss* GModelSpectralGauss::clone(void) const
240 240 * Evaluates
241 241 *
242 242 * \f[
243   - * \frac{dN}{dE}=\frac{\norm}{\sqrt{2\pi}\sigma}\exp(\frac{-(E-\bar{E})^2}{2\sigma^2})
  243 + * S_{\rm E}(E | t) = \frac{m\_norm}{\sqrt{2\pi}m\_sigma}
  244 + * \exp(\frac{-(E-m\_mean)^2}{2 m\_sigma^2})
244 245 * \f]
  246 + *
  247 + * where
  248 + * - \f${\tt m\_norm}\f$ is the normalization,
  249 + * - \f${\tt m\_mean}\f$ is the mean energy, and
  250 + * - \f${\tt m\_sigma}\f$ is the energy width.
245 251 ***************************************************************************/
246 252 double GModelSpectralGauss::eval(const GEnergy& srcEng,
247 253 const GTime& srcTime) const
248 254 {
  255 + // Get parameter values
249 256 double energy = srcEng.MeV();
250 257 double norm = m_norm.value();
251 258 double mean = m_mean.value();
252 259 double sigma = m_sigma.value();
253 260  
254 261 // Compute function value
  262 + double delta = energy - mean;
255 263 double term1 = (norm / sigma) * gammalib::inv_sqrt2pi;
256   - double term2 = (energy - mean) * (energy - mean) / (2.0 * sigma * sigma);
  264 + double term2 = delta * delta / (2.0 * sigma * sigma);
257 265 double value = term1 * std::exp(-term2);
258 266  
  267 + // Compile option: Check for NaN/Inf
  268 + #if defined(G_NAN_CHECK)
  269 + if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
  270 + std::cout << "*** ERROR: GModelSpectralGauss::eval";
  271 + std::cout << "(srcEng=" << srcEng;
  272 + std::cout << ", srcTime=" << srcTime << "):";
  273 + std::cout << " NaN/Inf encountered";
  274 + std::cout << " (value=" << value;
  275 + std::cout << ")" << std::endl;
  276 + }
  277 + #endif
  278 +
259 279 // Return
260 280 return value;
261 281 }
... ... @@ -268,64 +288,89 @@ double GModelSpectralGauss::eval(const GEnergy&amp; srcEng,
268 288 * @param[in] srcTime True photon arrival time.
269 289 * @return Model value (ph/cm2/s/MeV).
270 290 *
271   - * This method simply calls the eval() method as no analytical gradients will
272   - * be computed. See the eval() method for details.
  291 + * Evaluates
  292 + *
  293 + * \f[
  294 + * S_{\rm E}(E | t) = \frac{m\_norm}{\sqrt{2\pi}m\_sigma}
  295 + * \exp(\frac{-(E-m\_mean)^2}{2 m\_sigma^2})
  296 + * \f]
  297 + *
  298 + * where
  299 + * - \f${\tt m\_norm}\f$ is the normalization,
  300 + * - \f${\tt m\_mean}\f$ is the mean energy, and
  301 + * - \f${\tt m\_sigma}\f$ is the energy width.
  302 + *
  303 + * The method also evaluates the partial derivatives of the model with
  304 + * respect to the parameters using
  305 + *
  306 + * \f[
  307 + * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_norm}} =
  308 + * \frac{S_{\rm E}(E | t)}{{\tt m\_norm}}
  309 + * \f]
  310 + *
  311 + * \f[
  312 + * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_mean}} =
  313 + * S_{\rm E}(E | t) \frac{E-m\_mean}{m\_sigma^2}
  314 + * \f]
  315 + *
  316 + * \f[
  317 + * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_sigma}} =
  318 + * \frac{S_{\rm E}(E | t)}{{\tt m\_sigma}}
  319 + * \left( \frac{(E-m\_mean)^2}{m\_sigma^2} - 1 \right)
  320 + * \f]
273 321 *
274   - * TODO: update docstring.
275 322 ***************************************************************************/
276 323 double GModelSpectralGauss::eval_gradients(const GEnergy& srcEng,
277 324 const GTime& srcTime)
278 325 {
279   - double energy = srcEng.MeV();
280   - double norm = m_norm.value();
281   - double mean = m_mean.value();
282   - double sigma = m_sigma.value();
283   -
284   - // Update the evaluation cache
285   - update_eval_cache(srcEng);
286   -
287   - // Compute function terms
288   - double term1 = (norm / sigma) * gammalib::inv_sqrt2pi;
289   - double term2 = (1 / sigma) * gammalib::inv_sqrt2pi;
290   - double term3 = (energy - mean) * (energy - mean) / (2 * sigma * sigma);
291   - double term4 = (energy - mean) / (sigma * sigma);
292   - double term5 = (norm / (sigma * sigma)) * gammalib::inv_sqrt2pi;
293   -
294   - // Compute function value
295   - double value = term1 * std::exp(- term3);
296   -
297   -
298   - // Compute partial derivatives with respect to the parameter factor
299   - // values (partial differentials were determined analytically).
300   -
301   - double g_norm = term2 * std::exp(- term3);
302   - double g_mean = term1 * term4 * std::exp(- term3);
303   - double g_sigma = - term5 * std::exp(- term3) * (1 - (2 * term3));
304   -
305   - // Set gradients
306   - m_norm.factor_gradient(g_norm);
307   - m_mean.factor_gradient(g_mean);
308   - m_sigma.factor_gradient(g_sigma);
309   -
310   - // Compile option: Check for NaN/Inf
311   - #if defined(G_NAN_CHECK)
312   - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
313   - std::cout << "*** ERROR: GModelSpectralGauss::eval_gradients";
314   - std::cout << "(srcEng=" << srcEng;
315   - std::cout << ", srcTime=" << srcTime << "):";
316   - std::cout << " NaN/Inf encountered";
317   - std::cout << " (value=" << value;
318   - std::cout << ", norm=" << m_last_norm;
319   - std::cout << ", mean=" << m_last_mean;
320   - std::cout << ", Sigma=" << m_last_sigma;
321   - std::cout << ")" << std::endl;
322   - }
323   - #endif
324   -
325   - // Return
326   - return value;
327   - }
  326 + // Get parameter values
  327 + double energy = srcEng.MeV();
  328 + double norm = m_norm.value();
  329 + double mean = m_mean.value();
  330 + double sigma = m_sigma.value();
  331 +
  332 + // Compute function terms
  333 + double delta = energy - mean;
  334 + double sigma2 = sigma * sigma;
  335 + double term2 = (1.0 / sigma) * gammalib::inv_sqrt2pi;
  336 + double term1 = norm * term2;
  337 + double term3 = delta * delta / (2.0 * sigma2);
  338 + double term4 = delta / sigma2;
  339 + double term5 = (norm / sigma2) * gammalib::inv_sqrt2pi;
  340 + double eterm3 = std::exp(-term3);
  341 +
  342 + // Compute function value
  343 + double value = term1 * eterm3;
  344 +
  345 + // Compute partial derivatives with respect to the parameter factor
  346 + // values (partial differentials were determined analytically).
  347 + double g_norm = (m_norm.is_free())
  348 + ? term2 * eterm3 * m_norm.scale() : 0.0;
  349 + double g_mean = (m_mean.is_free())
  350 + ? value * term4 * m_mean.scale() : 0.0;
  351 + double g_sigma = (m_sigma.is_free())
  352 + ? -term5 * eterm3 * (1.0 - (2.0 * term3)) * m_sigma.scale()
  353 + : 0.0;
  354 +
  355 + // Set gradients
  356 + m_norm.factor_gradient(g_norm);
  357 + m_mean.factor_gradient(g_mean);
  358 + m_sigma.factor_gradient(g_sigma);
  359 +
  360 + // Compile option: Check for NaN/Inf
  361 + #if defined(G_NAN_CHECK)
  362 + if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
  363 + std::cout << "*** ERROR: GModelSpectralGauss::eval_gradients";
  364 + std::cout << "(srcEng=" << srcEng;
  365 + std::cout << ", srcTime=" << srcTime << "):";
  366 + std::cout << " NaN/Inf encountered";
  367 + std::cout << " (value=" << value;
  368 + std::cout << ")" << std::endl;
  369 + }
  370 + #endif
328 371  
  372 + // Return
  373 + return value;
329 374 }
330 375  
331 376  
... ... @@ -347,7 +392,6 @@ double GModelSpectralGauss::eval_gradients(const GEnergy&amp; srcEng,
347 392 * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
348 393 * The integration is done analytically.
349 394 ***************************************************************************/
350   -
351 395 double GModelSpectralGauss::flux(const GEnergy& emin,
352 396 const GEnergy& emax) const
353 397 {
... ... @@ -451,9 +495,12 @@ GEnergy GModelSpectralGauss::mc(const GEnergy&amp; emin,
451 495 const GTime& time,
452 496 GRan& ran) const
453 497 {
  498 + // Get energy boundaries in MeV
454 499 double xmax = emax.MeV();
455 500 double xmin = emin.MeV();
456   - double energy;
  501 +
  502 + // Initialize return energy
  503 + double energy = 0.0;
457 504  
458 505 // Throw an exception if energy range is invalid
459 506 if (xmin >= xmax) {
... ... @@ -461,33 +508,30 @@ GEnergy GModelSpectralGauss::mc(const GEnergy&amp; emin,
461 508 "Minimum energy < maximum energy required.");
462 509 }
463 510  
464   - // Get uniform value between 0 and 1
465   -
466   - // Do ...
  511 + // Sample until we find a value within the requested energy range
467 512 do {
468 513 double x1;
469 514 double w;
470   - do{
  515 + do {
471 516 x1 = 2.0 * ran.uniform() - 1.0;
472 517 double x2 = 2.0 * ran.uniform() - 1.0;
473 518 w = x1 * x1 + x2 * x2;
474 519 } while (w >= 1.0);
475 520  
476   - // Do ...
477   - w = std::sqrt( (-2.0 * std::log( w ) ) / w );
  521 + // Compute random value
  522 + double val = x1 * std::sqrt((-2.0 * std::log(w)) / w);
478 523  
479   - // Compute ...
480   - double val = x1 * w;
481   -
482   - // Map into [emin,emax] range
  524 + // Scale to specified width and shift by mean value
483 525 energy = m_sigma.value() * val + m_mean.value();
484 526  
485   - } while ((xmin <= energy) && (energy < xmax));
  527 + } while (energy < xmin || energy > xmax);
  528 +
486 529  
487 530 // Return energy
488 531 return GEnergy(energy, "MeV");
489 532 }
490 533  
  534 +
491 535 /***********************************************************************//**
492 536 * @brief Read model from XML element
493 537 *
... ... @@ -502,9 +546,9 @@ GEnergy GModelSpectralGauss::mc(const GEnergy&amp; emin,
502 546 * The format of the XML elements is:
503 547 *
504 548 * <spectrum type="Gaussian">
505   - * <parameter name="Prefactor" scale=".." value=".." min=".." max=".." free=".."/>
506   - * <parameter name="Mean" scale=".." value=".." min=".." max=".." free=".."/>
507   - * <parameter name="Sigma" scale=".." value=".." min=".." max=".." free=".."/>
  549 + * <parameter name="Normalization" scale=".." value=".." min=".." max=".." free=".."/>
  550 + * <parameter name="Mean" scale=".." value=".." min=".." max=".." free=".."/>
  551 + * <parameter name="Sigma" scale=".." value=".." min=".." max=".." free=".."/>
508 552 * </spectrum>
509 553 *
510 554 ***************************************************************************/
... ... @@ -526,8 +570,8 @@ void GModelSpectralGauss::read(const GXmlElement&amp; xml)
526 570 // Get parameter element
527 571 const GXmlElement* par = xml.element("parameter", i);
528 572  
529   - // Handle prefactor
530   - if (par->attribute("name") == "Prefactor") {
  573 + // Handle normalization
  574 + if (par->attribute("name") == "Normalization") {
531 575 m_norm.read(*par);
532 576 npar[0]++;
533 577 }
... ... @@ -549,7 +593,7 @@ void GModelSpectralGauss::read(const GXmlElement&amp; xml)
549 593 // Verify that all parameters were found
550 594 if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1) {
551 595 throw GException::model_invalid_parnames(G_READ, xml,
552   - "Require \"Prefactor\", \"Mean\", and \"Sigma\""
  596 + "Require \"Normalization\", \"Mean\", and \"Sigma\""
553 597 " parameters.");
554 598 }
555 599  
... ... @@ -557,6 +601,7 @@ void GModelSpectralGauss::read(const GXmlElement&amp; xml)
557 601 return;
558 602 }
559 603  
  604 +
560 605 /***********************************************************************//**
561 606 * @brief Write model into XML element
562 607 *
... ... @@ -573,9 +618,9 @@ void GModelSpectralGauss::read(const GXmlElement&amp; xml)
573 618 * element is
574 619 *
575 620 * <spectrum type="Gaussian">
576   - * <parameter name="Prefactor" scale=".." value=".." min=".." max=".." free=".."/>
577   - * <parameter name="Mean" scale=".." value=".." min=".." max=".." free=".."/>
578   - * <parameter name="Sigma" scale=".." value=".." min=".." max=".." free=".."/>
  621 + * <parameter name="Normalization" scale=".." value=".." min=".." max=".." free=".."/>
  622 + * <parameter name="Mean" scale=".." value=".." min=".." max=".." free=".."/>
  623 + * <parameter name="Sigma" scale=".." value=".." min=".." max=".." free=".."/>
579 624 * </spectrum>
580 625 ***************************************************************************/
581 626 void GModelSpectralGauss::write(GXmlElement& xml) const
... ... @@ -593,7 +638,7 @@ void GModelSpectralGauss::write(GXmlElement&amp; xml) const
593 638  
594 639 // If XML element has 0 nodes then append 3 parameter nodes
595 640 if (xml.elements() == 0) {
596   - xml.append(GXmlElement("parameter name=\"Prefactor\""));
  641 + xml.append(GXmlElement("parameter name=\"Normalization\""));
597 642 xml.append(GXmlElement("parameter name=\"Mean\""));
598 643 xml.append(GXmlElement("parameter name=\"Sigma\""));
599 644 }
... ... @@ -611,8 +656,8 @@ void GModelSpectralGauss::write(GXmlElement&amp; xml) const
611 656 // Get parameter element
612 657 GXmlElement* par = xml.element("parameter", i);
613 658  
614   - // Handle prefactor
615   - if (par->attribute("name") == "Prefactor") {
  659 + // Handle normalization
  660 + if (par->attribute("name") == "Normalization") {
616 661 npar[0]++;
617 662 m_norm.write(*par);
618 663 }
... ... @@ -635,7 +680,7 @@ void GModelSpectralGauss::write(GXmlElement&amp; xml) const
635 680 // Check of all required parameters are present
636 681 if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1 ) {
637 682 throw GException::model_invalid_parnames(G_WRITE, xml,
638   - "Require \"Prefactor\", \"Mean\" and \"Sigma\""
  683 + "Require \"Normalization\", \"Mean\" and \"Sigma\""
639 684 " parameters.");
640 685 }
641 686  
... ... @@ -688,8 +733,8 @@ void GModelSpectralGauss::init_members(void)
688 733 {
689 734 // Initialise normalisation
690 735 m_norm.clear();
691   - m_norm.name("Prefactor");
692   - m_norm.unit("ph/cm2/s/MeV");
  736 + m_norm.name("Normalization");
  737 + m_norm.unit("ph/cm2/s");
693 738 m_norm.scale(1.0);
694 739 m_norm.value(1.0); // default: 1.0
695 740 m_norm.min(0.0); // min: 0.0
... ... @@ -719,11 +764,6 @@ void GModelSpectralGauss::init_members(void)
719 764 m_sigma.gradient(0.0);
720 765 m_sigma.has_grad(false);
721 766  
722   - // Initialise eval cache
723   - m_last_norm = 0.0;
724   - m_last_mean = 0.0;
725   - m_last_sigma = 0.0;
726   -
727 767 // Set parameter pointer(s)
728 768 m_pars.clear();
729 769 m_pars.push_back(&m_norm);
... ... @@ -753,11 +793,6 @@ void GModelSpectralGauss::copy_members(const GModelSpectralGauss&amp; model)
753 793 m_pars.push_back(&m_mean);
754 794 m_pars.push_back(&m_sigma);
755 795  
756   - // Copy eval cache
757   - m_last_norm = model.m_last_norm;
758   - m_last_mean = model.m_last_mean;
759   - m_last_sigma = model.m_last_sigma;
760   -
761 796 // Return
762 797 return;
763 798 }
... ... @@ -780,51 +815,11 @@ void GModelSpectralGauss::free_members(void)
780 815 double GModelSpectralGauss::eflux_kernel::eval(const double& energy)
781 816 {
782 817 // Evaluate function value
  818 + double delta = energy - m_mean;
783 819 double term1 = (m_norm / m_sigma) * gammalib::inv_sqrt2pi;
784   - double term2 = (energy - m_mean) * (energy - m_mean) / (2.0 * m_sigma * m_sigma);
  820 + double term2 = delta * delta / (2.0 * m_sigma * m_sigma);
785 821 double value = term1 * std::exp(- term2);
786 822  
787 823 // Return value
788 824 return value;
789 825 }
790   -
791   -/***********************************************************************//**
792   - * @brief Update eval precomputation cache
793   - *
794   - * @param[in] energy Energy.
795   - *
796   - * Updates the precomputation cache for eval() and eval_gradients() methods.
797   - *
798   - ***************************************************************************/
799   -void GModelSpectralGauss::update_eval_cache(const GEnergy& energy) const
800   -{
801   - // Get parameter values (takes 3 multiplications which are difficult
802   - // to avoid)
803   - double norm = m_norm.value();
804   - double mean = m_mean.value();
805   - double sigma = m_sigma.value();
806   -
807   - // If the energy or one of the parameters norm, mean or sigma
808   - // energy has changed then recompute the cache
809   - if ((m_last_norm != norm) ||
810   - (m_last_mean != mean) ||
811   - (m_last_sigma != sigma)) {
812   -
813   - // Store actual energy and parameter values
814   - m_last_norm = norm;
815   - m_last_mean = mean;
816   - m_last_sigma = sigma;
817   -
818   - // Compute and store value
819   - double eng = energy.MeV();
820   - m_last_norm = eng / m_last_norm;
821   - m_last_mean = eng / m_last_mean;
822   - m_last_sigma = eng/ m_last_sigma;
823   -
824   - } // endif: recomputation was required
825   -
826   - // Return
827   - return;
828   -}
829   -
830   -
... ...
test/data/model_point_gauss.xml
  View file @244c69d
... ... @@ -2,9 +2,9 @@
2 2 <source_library title="source library">
3 3 <source name="Crab" type="PointSource">
4 4 <spectrum type="Gaussian">
5   - <parameter name="Prefactor" scale="1e-16" value="5.7" min="1e-07" max="1000.0" free="1"/>
6   - <parameter name="Mean" scale="1e6" value="2.48" min="0.0" max="+5.0" free="1"/>
7   - <parameter name="Sigma" scale="1e6" value="1.0" min="0.0001" max="1000.0" free="1"/>
  5 + <parameter name="Normalization" scale="1e-10" value="1.0" min="1e-07" max="1000.0" free="1"/>
  6 + <parameter name="Mean" scale="1e6" value="5.0" min="0.01" max="100.0" free="1"/>
  7 + <parameter name="Sigma" scale="1e6" value="1.0" min="0.01" max="100.0" free="1"/>
8 8 </spectrum>
9 9 <spatialModel type="SkyDirFunction">
10 10 <parameter free="0" max="360" min="-360" name="RA" scale="1" value="83.6331" />
... ...
test/test_GModel.cpp
  View file @244c69d
... ... @@ -1176,7 +1176,7 @@ void TestGModel::test_gauss(void)
1176 1176 test_try("Test value constructor");
1177 1177 try {
1178 1178 GModelSpectralGauss model(1.0, GEnergy(42.0, "MeV"), GEnergy(0.5, "MeV"));
1179   - test_value(model.prefactor(), 1.0);
  1179 + test_value(model.norm(), 1.0);
1180 1180 test_value(model.mean().MeV(), 42.0);
1181 1181 test_value(model.sigma().MeV(), 0.5);
1182 1182 test_try_success();
... ... @@ -1189,29 +1189,29 @@ void TestGModel::test_gauss(void)
1189 1189 test_try("Test XML constructor, value and gradients");
1190 1190 try {
1191 1191 // Test XML constructor
1192   - GXml xml(m_xml_model_point_gauss);
1193   - GXmlElement* element = xml.element(0)->element(0)->element("spectrum", 0);
  1192 + GXml xml(m_xml_model_point_gauss);
  1193 + GXmlElement* element = xml.element(0)->element(0)->element("spectrum", 0);
1194 1194 GModelSpectralGauss model(*element);
1195 1195 test_value(model.size(), 3);
1196 1196 test_assert(model.type() == "Gaussian", "Expected \"Gaussian\"");
1197   - test_value(model.prefactor(), 5.7e-16);
1198   - test_value(model.mean().MeV(), 2.48e6);
1199   - test_value(model.sigma().MeV(), 1e6);
  1197 + test_value(model.norm(), 1.0e-10);
  1198 + test_value(model.mean().MeV(), 5.0e6);
  1199 + test_value(model.sigma().MeV(), 1.0e6);
1200 1200  
1201   - // Test prefactor method
1202   - model.prefactor(2.3e-16);
1203   - test_value(model.prefactor(), 2.3e-16);
  1201 + // Test norm method
  1202 + model.norm(2.3e-10);
  1203 + test_value(model.norm(), 2.3e-10);
1204 1204  
1205 1205 // Test mean method
1206   - model.mean(GEnergy(96.0, "MeV"));
1207   - test_value(model.mean().MeV(), 96.0);
  1206 + model.mean(GEnergy(1.5678, "TeV"));
  1207 + test_value(model.mean().TeV(), 1.5678);
1208 1208  
1209 1209 // Test sigma method
1210   - model.sigma(GEnergy(0.5e7, "keV"));
1211   - test_value(model.sigma().keV(), 0.5e7);
  1210 + model.sigma(GEnergy(0.1234, "TeV"));
  1211 + test_value(model.sigma().TeV(), 0.1234);
1212 1212  
1213 1213 // Test operator access
1214   - const char* strarray[] = {"Prefactor", "Mean", "Sigma"};
  1214 + const char* strarray[] = {"Normalization", "Mean", "Sigma"};
1215 1215 for (int i = 0; i < 3; ++i) {
1216 1216 std::string keyname(strarray[i]);
1217 1217 model[keyname].remove_range(); // To allow setting of any value
... ...