Merge pull request #2326 from JouvinLea/fix_gaussianmodels

Fix bug in the Gaussian model evaluate method
gammapy · Aug 27, 2019 · 59f7c13 · 59f7c13
2 parents 71a0225 + 8d33eda
commit 59f7c13
Show file tree

Hide file tree

Showing 2 changed files with 10 additions and 3 deletions.
diff --git a/gammapy/spectrum/models.py b/gammapy/spectrum/models.py
@@ -1699,7 +1699,7 @@ class SpectralGaussian(SpectralModel):
 
     .. math::
 
-        \phi(E) = \frac{N_0}{\sigma \sqrt{2\pi}}  \exp{ \frac{\left( E-\bar{E} \right)^2 }{2 \sigma^2} }
+        \phi(E) = \frac{N_0}{\sigma \sqrt{2\pi}}  \exp{ \frac{- \left( E-\bar{E} \right)^2 }{2 \sigma^2} }
 
 
 
@@ -1741,7 +1741,7 @@ def evaluate(energy, norm, mean, sigma):
         return (
             norm
             / (sigma * np.sqrt(2 * np.pi))
-            * np.exp((energy - mean) ** 2 / (2 * sigma ** 2))
+            * np.exp(-(energy - mean) ** 2 / (2 * sigma ** 2))
         )
 
     def integral(self, emin, emax, **kwargs):
@@ -1806,7 +1806,7 @@ class SpectralLogGaussian(SpectralModel):
     .. math::
 
         \phi(E) = \frac{N_0}{E \, \sigma \sqrt{2\pi}}
-         \exp{ \frac{\left( \ln(\frac{E}{\bar{E}}) \right)^2 }{2 \sigma^2} }
+         \exp{ \frac{- \left( \ln(\frac{E}{\bar{E}}) \right)^2 }{2 \sigma^2} }
 
     This model was used in this CTA study for the electron spectrum: Table 3
      in https://ui.adsabs.harvard.edu/abs/2013APh....43..171B

diff --git a/gammapy/spectrum/tests/test_models.py b/gammapy/spectrum/tests/test_models.py
@@ -172,6 +172,7 @@ def table_model():
             norm=4 / u.cm ** 2 / u.s, mean=2 * u.TeV, sigma=0.2 * u.TeV
         ),
         val_at_2TeV=u.Quantity(7.978845608028654, "cm-2 s-1 TeV-1"),
+        val_at_3TeV=u.Quantity(2.973439029468601e-05, "cm-2 s-1 TeV-1"),
         integral_1_10TeV=u.Quantity(3.9999988533937123, "cm-2 s-1"),
         integral_infinity=u.Quantity(4, "cm-2 s-1"),
         eflux_1_10TeV=u.Quantity(7.999998896163037, "TeV cm-2 s-1"),
@@ -180,6 +181,7 @@ def table_model():
         name="SpectralLogGaussian",
         model=SpectralLogGaussian(norm=4 / u.cm ** 2 / u.s, mean=2 * u.TeV, sigma=0.2),
         val_at_2TeV=u.Quantity(3.98942280401, "cm-2 s-1 TeV-1"),
+        val_at_3TeV=u.Quantity(0.34066933236079916, "cm-2 s-1 TeV-1"),
         integral_1_10TeV=u.Quantity(3.994439, "cm-2 s-1"),
         eflux_1_10TeV=u.Quantity(8.151414, "TeV cm-2 s-1"),
     ),
@@ -271,6 +273,11 @@ def test_models(spectrum):
     energy = 2 * u.TeV
     value = model(energy)
     assert_quantity_allclose(value, spectrum["val_at_2TeV"])
+    if "val_at_3TeV" in spectrum:
+        energy = 3 * u.TeV
+        value = model(energy)
+        assert_quantity_allclose(value, spectrum["val_at_3TeV"])
+
     emin = 1 * u.TeV
     emax = 10 * u.TeV
     assert_quantity_allclose(