Repository: gammapy/gammapy
Pull request #5661 (2025-01-07 15:48:43)

Backport PR #5645 on branch v1.3.x (Fix `expfactor` unit to `MeV**-index_2` in SuperExpCutoffPowerLaw4FGLSpectralModel)
Backport PR #5645: Fix `expfactor` unit to `MeV**-index_2` in SuperExpCutoffPowerLaw4FGLSpectralModel


Patch:
diff --git a/codemeta.json b/codemeta.json
index 642a291a5..09f305e55 100644
--- a/codemeta.json
+++ b/codemeta.json
@@ -359,7 +359,7 @@
         "active"
     ],
     "email": "GAMMAPY-COORDINATION-L@IN2P3.FR",
-    "dateModified": "2024-11-26",
+    "dateModified": "2025-01-07",
     "softwareRequirements": [
         "numpy>=1.21",
         "scipy>=1.5,!=1.10",
diff --git a/examples/models/spectral/plot_super_exp_cutoff_powerlaw_4fgl_dr1.py b/examples/models/spectral/plot_super_exp_cutoff_powerlaw_4fgl_dr1.py
index f50f8727f..4de87bb34 100644
--- a/examples/models/spectral/plot_super_exp_cutoff_powerlaw_4fgl_dr1.py
+++ b/examples/models/spectral/plot_super_exp_cutoff_powerlaw_4fgl_dr1.py
@@ -36,7 +36,7 @@ model = SuperExpCutoffPowerLaw4FGLSpectralModel(
     index_2=2,
     amplitude="1e-12 TeV-1 cm-2 s-1",
     reference="1 TeV",
-    expfactor=1e-2,
+    expfactor=1e-14,
 )
 model.plot(energy_range)
 plt.grid(which="both")
diff --git a/gammapy/modeling/models/spectral.py b/gammapy/modeling/models/spectral.py
index ac0942644..5552741c0 100644
--- a/gammapy/modeling/models/spectral.py
+++ b/gammapy/modeling/models/spectral.py
@@ -1526,8 +1526,8 @@ class SuperExpCutoffPowerLaw4FGLSpectralModel(SpectralModel):
         :math:`E_0`. Default is 1 TeV.
     expfactor : `~astropy.units.Quantity`
         :math:`a`, given as dimensionless value but
-        internally assumes unit of :math:`E_0^{-\Gamma_2}`.
-        Default is 1e-2.
+        internally assumes unit of :math:`{\rm MeV}^{-\Gamma_2}`.
+        Default is 1e-14.
     """
 
     tag = ["SuperExpCutoffPowerLaw4FGLSpectralModel", "secpl-4fgl"]
@@ -1538,7 +1538,7 @@ class SuperExpCutoffPowerLaw4FGLSpectralModel(SpectralModel):
         interp="log",
     )
     reference = Parameter("reference", "1 TeV", frozen=True)
-    expfactor = Parameter("expfactor", "1e-2")
+    expfactor = Parameter("expfactor", "1e-14")
     index_1 = Parameter("index_1", 1.5)
     index_2 = Parameter("index_2", 2)
 
@@ -1552,7 +1552,7 @@ class SuperExpCutoffPowerLaw4FGLSpectralModel(SpectralModel):
 
         pwl = amplitude * (energy / reference) ** (-index_1)
         cutoff = np.exp(
-            expfactor / reference.unit**index_2 * (reference**index_2 - energy**index_2)
+            expfactor / u.MeV**index_2 * (reference**index_2 - energy**index_2)
         )
         return pwl * cutoff
 


Test patch:
diff --git a/gammapy/catalog/tests/test_fermi.py b/gammapy/catalog/tests/test_fermi.py
index 15d8290ae..b4e60ca11 100644
--- a/gammapy/catalog/tests/test_fermi.py
+++ b/gammapy/catalog/tests/test_fermi.py
@@ -86,6 +86,7 @@ SOURCES_4FGL = [
         spec_type=PowerLawSpectralModel,
         dnde=u.Quantity(2.9476e-11, "cm-2 s-1 GeV-1"),
         dnde_err=u.Quantity(5.3318e-12, "cm-2 s-1 GeV-1"),
+        dnde_10GeV=u.Quantity(1.470835e-12, "cm-2 s-1 GeV-1"),
     ),
     dict(
         idx=3,
@@ -94,6 +95,7 @@ SOURCES_4FGL = [
         spec_type=LogParabolaSpectralModel,
         dnde=u.Quantity(2.8545e-8, "cm-2 s-1 GeV-1"),
         dnde_err=u.Quantity(1.3324e-9, "cm-2 s-1 GeV-1"),
+        dnde_10GeV=u.Quantity(1.159598e-12, "cm-2 s-1 GeV-1"),
     ),
     dict(
         idx=7,
@@ -102,6 +104,7 @@ SOURCES_4FGL = [
         spec_type=SuperExpCutoffPowerLaw4FGLSpectralModel,
         dnde=u.Quantity(2.084e-09, "cm-2 s-1 GeV-1"),
         dnde_err=u.Quantity(1.0885e-10, "cm-2 s-1 GeV-1"),
+        dnde_10GeV=u.Quantity(2.201476e-12, "cm-2 s-1 GeV-1"),
     ),
     dict(
         idx=2718,
@@ -110,6 +113,7 @@ SOURCES_4FGL = [
         spec_type=LogParabolaSpectralModel,
         dnde=u.Quantity(1.3237202133031811e-12, "cm-2 s-1 MeV-1"),
         dnde_err=u.Quantity(4.513233455580648e-14, "cm-2 s-1 MeV-1"),
+        dnde_10GeV=u.Quantity(2.234733e-10, "cm-2 s-1 GeV-1"),
     ),
 ]
 
@@ -246,9 +250,11 @@ class TestFermi4FGLObject:
 
         e_ref = model.reference.quantity
         dnde, dnde_err = model.evaluate_error(e_ref)
+        dnde_10GeV = model(10 * u.GeV)
         assert isinstance(model, ref["spec_type"])
         assert_quantity_allclose(dnde, ref["dnde"], rtol=1e-4)
         assert_quantity_allclose(dnde_err, ref["dnde_err"], rtol=1e-4)
+        assert_quantity_allclose(dnde_10GeV, ref["dnde_10GeV"], rtol=1e-4)
 
     def test_spatial_model(self):
         model = self.cat["4FGL J0000.3-7355"].spatial_model()
diff --git a/gammapy/modeling/models/tests/test_spectral.py b/gammapy/modeling/models/tests/test_spectral.py
index 89526d212..86b4e6f41 100644
--- a/gammapy/modeling/models/tests/test_spectral.py
+++ b/gammapy/modeling/models/tests/test_spectral.py
@@ -150,7 +150,7 @@ TEST_MODELS = [
             index_2=2,
             amplitude=1 / u.cm**2 / u.s / u.TeV,
             reference=1 * u.TeV,
-            expfactor=1e-2,
+            expfactor=1e-14,
         ),
         val_at_2TeV=u.Quantity(0.3431043087721737, "cm-2 s-1 TeV-1"),
         integral_1_10TeV=u.Quantity(1.2125247, "cm-2 s-1"),
@@ -683,6 +683,18 @@ def test_ecpl_integrate():
     assert_quantity_allclose(value, 8.380714e-14 * u.Unit("s-1 cm-2"))
 
 
+def test_call_plsec_4fgl_dr1():
+    model = SuperExpCutoffPowerLaw4FGLSpectralModel(
+        amplitude="2e-12 MeV-1 s-1 cm-2",
+        reference="1000 MeV",
+        expfactor=5e-03,
+        index_1=2,
+        index_2=1,
+    )
+    desired = u.Quantity(3.823779e-20, "MeV-1 s-1 cm-2")
+    assert_allclose(model(4 * u.GeV), desired, rtol=1e-3)
+
+
 def test_pwl_pivot_energy():
     pwl = PowerLawSpectralModel(amplitude="5.35510540e-11 cm-2 s-1 TeV-1")
     assert_quantity_allclose(pwl.pivot_energy, np.nan * u.TeV, rtol=1e-5)