Domenico Tiziani / gammalib

Commit b64bc21827e87b258cac7587289f9d37da133265

Authored by Jürgen Knödlseder
1 parent 6b0f4241

Code formating before merging into trunk

Inline Side-by-side
Showing 5 changed files with 102 additions and 100 deletions
include/GModelSpectralGauss.hpp
  View file @b64bc21
1 1 /***************************************************************************
2 2 * GModelSpectralGauss.hpp - Spectral Gaussian model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2013 by Juergen Knoedlseder *
  4 + * copyright (C) 2014 by Christoph Deil & Ellis Owen *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -60,9 +60,9 @@ public:
60 60 // Constructors and destructors
61 61 GModelSpectralGauss(void);
62 62 explicit GModelSpectralGauss(const GXmlElement& xml);
63   - explicit GModelSpectralGauss(const double& prefactor,
64   - const GEnergy& mean,
65   - const GEnergy& sigma);
  63 + GModelSpectralGauss(const double& prefactor,
  64 + const GEnergy& mean,
  65 + const GEnergy& sigma);
66 66 GModelSpectralGauss(const GModelSpectralGauss& model);
67 67 virtual ~GModelSpectralGauss(void);
68 68  
... ... @@ -218,5 +218,4 @@ void GModelSpectralGauss::sigma(const GEnergy& sigma)
218 218 return;
219 219 }
220 220  
221   -
222 221 #endif /* GMODELSPECTRALGAUSS_HPP */
... ...
pyext/GModelSpectralGauss.i
  View file @b64bc21
1 1 /***************************************************************************
2 2 * GModelSpectralGauss.i - Spectral gaussian model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2013 by Juergen Knoedlseder *
  4 + * copyright (C) 2014 by Christoph Deil & Ellis Owen *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -41,9 +41,9 @@ public:
41 41 // Constructors and destructors
42 42 GModelSpectralGauss(void);
43 43 explicit GModelSpectralGauss(const GXmlElement& xml);
44   - explicit GModelSpectralGauss(const double& prefactor,
45   - const GEnergy& mean,
46   - const GEnergy& sigma);
  44 + GModelSpectralGauss(const double& prefactor,
  45 + const GEnergy& mean,
  46 + const GEnergy& sigma);
47 47 GModelSpectralGauss(const GModelSpectralGauss& model);
48 48 virtual ~GModelSpectralGauss(void);
49 49  
... ... @@ -58,12 +58,11 @@ public:
58 58 virtual double flux(const GEnergy& emin,
59 59 const GEnergy& emax) const;
60 60 virtual double eflux(const GEnergy& emin,
61   -
62 61 const GEnergy& emax) const;
63   - virtual GEnergy mc(const GEnergy& emin,
64   - const GEnergy& emax,
65   - const GTime& time,
66   - GRan& ran) const;
  62 + virtual GEnergy mc(const GEnergy& emin,
  63 + const GEnergy& emax,
  64 + const GTime& time,
  65 + GRan& ran) const;
67 66 virtual void read(const GXmlElement& xml);
68 67 virtual void write(GXmlElement& xml) const;
69 68  
... ...
src/model/GModelSpectralGauss.cpp
  View file @b64bc21
1 1 /***************************************************************************
2 2 * GModelSpectralGauss.cpp - Spectral Gaussian model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2013 by Juergen Knoedlseder *
  4 + * copyright (C) 2014 by Christoph Deil & Ellis Owen *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -108,7 +108,7 @@ GModelSpectralGauss::GModelSpectralGauss(const GXmlElement& xml) :
108 108 GModelSpectralGauss::GModelSpectralGauss(const double& prefactor,
109 109 const GEnergy& mean,
110 110 const GEnergy& sigma) :
111   - GModelSpectral()
  111 + GModelSpectral()
112 112 {
113 113 // Initialise members
114 114 init_members();
... ... @@ -240,21 +240,22 @@ 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 + * \frac{dN}{dE} = \frac{\norm}{\sqrt{2\pi}\sigma}
  244 + * \exp(\frac{(E-\bar{E})^2}{2\sigma^2})
244 245 * \f]
245 246 ***************************************************************************/
246 247 double GModelSpectralGauss::eval(const GEnergy& srcEng,
247 248 const GTime& srcTime) const
248 249 {
249 250 double energy = srcEng.MeV();
250   - double norm = m_norm.value();
251   - double mean = m_mean.value();
252   - double sigma = m_sigma.value();
  251 + double norm = m_norm.value();
  252 + double mean = m_mean.value();
  253 + double sigma = m_sigma.value();
253 254  
254 255 // Compute function value
255 256 double term1 = (norm / sigma) * gammalib::inv_sqrt2pi;
256   - double term2 = (energy - mean) * (energy - mean) / (2 * sigma * sigma);
257   - double value = term1 * std::exp(- term2);
  257 + double term2 = (energy - mean) * (energy - mean) / (2.0 * sigma * sigma);
  258 + double value = term1 * std::exp(-term2);
258 259  
259 260 // Return
260 261 return value;
... ... @@ -297,28 +298,27 @@ double GModelSpectralGauss::eval_gradients(const GEnergy& srcEng,
297 298 * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
298 299 * The integration is done analytically.
299 300 ***************************************************************************/
300   -
301 301 double GModelSpectralGauss::flux(const GEnergy& emin,
302   - const GEnergy& emax) const{
303   -
304   - double energy_min = emin.MeV();
305   - double energy_max = emax.MeV();
306   -
307   -
308   - double norm = m_norm.value();
309   - double mean = m_mean.value();
310   - double sigma = m_sigma.value();
311   -
312   - double zmin = (energy_min - mean)/sqrt(2.)/sigma;
313   - double zmax = (energy_max - mean)/sqrt(2.)/sigma;
314   -
  302 + const GEnergy& emax) const
  303 +{
  304 + // Initialise flux
315 305 double flux = 0.0;
316 306  
317 307 // Compute only if integration range is valid
318 308 if (emin < emax) {
319 309  
  310 + // Precomputations
  311 + double energy_min = emin.MeV();
  312 + double energy_max = emax.MeV();
  313 + double norm = m_norm.value();
  314 + double mean = m_mean.value();
  315 + double sigma = m_sigma.value();
  316 + double denom = 1.0 / (gammalib::sqrt_two*sigma);
  317 + double zmin = (energy_min - mean) * denom;
  318 + double zmax = (energy_max - mean) * denom;
  319 +
320 320 // Compute flux for a constant model
321   - flux = norm*(gammalib::erfcc(zmin) - gammalib::erfcc(zmax))/2.;
  321 + flux = norm*(gammalib::erfcc(zmin) - gammalib::erfcc(zmax))/2.0;
322 322  
323 323 } // endif: integration range was valid
324 324  
... ... @@ -413,21 +413,24 @@ GEnergy GModelSpectralGauss::mc(const GEnergy&amp; emin,
413 413  
414 414 // Get uniform value between 0 and 1
415 415 double u = ran.uniform();
416   - float x1, x2, w;
417 416  
  417 + // Do ...
418 418 do {
  419 + double x1;
  420 + double w;
419 421 do{
420   - x1 = 2.0 * u - 1.0;
421   - x2 = 2.0 * u - 1.0;
422   - w = x1 * x1 + x2 * x2;
  422 + x1 = 2.0 * u - 1.0;
  423 + double x2 = 2.0 * u - 1.0;
  424 + w = x1 * x1 + x2 * x2;
423 425 } while (w >= 1.0);
424 426  
425   - w = sqrt( (-2.0 * log( w ) ) / w );
  427 + // Do ...
  428 + w = std::sqrt( (-2.0 * std::log( w ) ) / w );
426 429  
  430 + // Compute ...
427 431 double val = x1 * w;
428 432  
429 433 // Map into [emin,emax] range
430   -
431 434 energy = (xmax - xmin) * val + xmin;
432 435  
433 436 } while ((xmin <= energy) && (energy < xmax));
... ... @@ -459,50 +462,51 @@ GEnergy GModelSpectralGauss::mc(const GEnergy&amp; emin,
459 462 ***************************************************************************/
460 463 void GModelSpectralGauss::read(const GXmlElement& xml)
461 464 {
462   - const int n_pars = 3;
463   -
464   - // Verify that XML element has exactly 3 parameters
465   - if (xml.elements() != n_pars || xml.elements("parameter") != n_pars) {
466   - throw GException::model_invalid_parnum(G_READ, xml,
467   - "Gaussian model requires exactly 3 parameters.");
468   - }
469   -
470   - // Extract model parameters
471   - int npar[] = {0, 0, 0};
472   - for (int i = 0; i < n_pars; ++i) {
473   -
474   - // Get parameter element
475   - const GXmlElement* par = xml.element("parameter", i);
476   -
477   - // Handle prefactor
478   - if (par->attribute("name") == "Prefactor") {
479   - m_norm.read(*par);
480   - npar[0]++;
481   - }
482   -
483   - // Handle mean
484   - else if (par->attribute("name") == "Mean") {
485   - m_mean.read(*par);
486   - npar[1]++;
487   - }
488   -
489   - // Handle sigma
490   - else if (par->attribute("name") == "Sigma") {
491   - m_sigma.read(*par);
492   - npar[2]++;
493   - }
494   -
495   - } // endfor: looped over all parameters
496   -
497   - // Verify that all parameters were found
498   - if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1) {
499   - throw GException::model_invalid_parnames(G_READ, xml,
500   - "Require \"Prefactor\", \"Mean\", and \"Sigma\""
501   - " parameters.");
502   - }
503   -
504   - // Return
505   - return;
  465 + // Set number of parameters
  466 + const int n_pars = 3;
  467 +
  468 + // Verify that XML element has exactly 3 parameters
  469 + if (xml.elements() != n_pars || xml.elements("parameter") != n_pars) {
  470 + throw GException::model_invalid_parnum(G_READ, xml,
  471 + "Gaussian model requires exactly 3 parameters.");
  472 + }
  473 +
  474 + // Extract model parameters
  475 + int npar[] = {0, 0, 0};
  476 + for (int i = 0; i < n_pars; ++i) {
  477 +
  478 + // Get parameter element
  479 + const GXmlElement* par = xml.element("parameter", i);
  480 +
  481 + // Handle prefactor
  482 + if (par->attribute("name") == "Prefactor") {
  483 + m_norm.read(*par);
  484 + npar[0]++;
  485 + }
  486 +
  487 + // Handle mean
  488 + else if (par->attribute("name") == "Mean") {
  489 + m_mean.read(*par);
  490 + npar[1]++;
  491 + }
  492 +
  493 + // Handle sigma
  494 + else if (par->attribute("name") == "Sigma") {
  495 + m_sigma.read(*par);
  496 + npar[2]++;
  497 + }
  498 +
  499 + } // endfor: looped over all parameters
  500 +
  501 + // Verify that all parameters were found
  502 + if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1) {
  503 + throw GException::model_invalid_parnames(G_READ, xml,
  504 + "Require \"Prefactor\", \"Mean\", and \"Sigma\""
  505 + " parameters.");
  506 + }
  507 +
  508 + // Return
  509 + return;
506 510 }
507 511  
508 512  
... ... @@ -687,8 +691,8 @@ void GModelSpectralGauss::init_members(void)
687 691 void GModelSpectralGauss::copy_members(const GModelSpectralGauss& model)
688 692 {
689 693 // Copy members
690   - m_norm = model.m_norm;
691   - m_mean = model.m_mean;
  694 + m_norm = model.m_norm;
  695 + m_mean = model.m_mean;
692 696 m_sigma = model.m_sigma;
693 697  
694 698 // Set parameter pointer(s)
... ... @@ -720,7 +724,7 @@ double GModelSpectralGauss::eflux_kernel::eval(const double&amp; energy)
720 724 {
721 725 // Evaluate function value
722 726 double term1 = (m_norm / m_sigma) * gammalib::inv_sqrt2pi;
723   - double term2 = (energy - m_mean) * (energy - m_mean) / (2 * m_sigma * m_sigma);
  727 + double term2 = (energy - m_mean) * (energy - m_mean) / (2.0 * m_sigma * m_sigma);
724 728 double value = term1 * std::exp(- term2);
725 729  
726 730 // Return value
... ...
test/data/model_point_gauss.xml
  View file @b64bc21
1 1 <?xml version="1.0" standalone="no"?>
2 2 <source_library title="source library">
3   - <source name="Crab" type="PointSource">
  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="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"/>
8 8 </spectrum>
9   - <spatialModel type="SkyDirFunction">
10   - <parameter free="0" max="360" min="-360" name="RA" scale="1" value="83.6331" />
11   - <parameter free="0" max="90" min="-90" name="DEC" scale="1" value="22.0145" />
12   - </spatialModel>
13   - </source>
  9 + <spatialModel type="SkyDirFunction">
  10 + <parameter free="0" max="360" min="-360" name="RA" scale="1" value="83.6331" />
  11 + <parameter free="0" max="90" min="-90" name="DEC" scale="1" value="22.0145" />
  12 + </spatialModel>
  13 + </source>
14 14 </source_library>
... ...
test/test_GModel.cpp
  View file @b64bc21
1 1 /***************************************************************************
2 2 * test_GModel.cpp - test GModel class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2013 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2014 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ...