Merge pull request #130 from fact-project/fix_unit

Fix unit conversion in power_law_integral

Author: maxnoe

fact/VERSION

@@ -1 +1 @@
-0.25.0
+0.25.1

fact/analysis/statistics.py

@@ -407,8 +407,8 @@ def power_law_integral(
     res = flux_normalization * e_ref / int_index * e_term
 
     if flux_normalization.unit.is_equivalent(FLUX_UNIT):
-        return res.to(FLUX_UNIT)
-    return res.to(POINT_SOURCE_FLUX_UNIT)
+        return res.to(FLUX_UNIT * u.GeV)
+    return res.to(POINT_SOURCE_FLUX_UNIT * u.GeV)
 
 
 @u.quantity_input

tests/test_analysis.py

@@ -26,3 +26,20 @@ def test_power():
 
     with raises(ValueError):
         random_power(2.7, 5 * u.GeV, 10 * u.GeV, e_ref=1 * u.GeV, size=1)
+
+
+def test_power_law_integral():
+    from fact.analysis.statistics import power_law_integral, FLUX_UNIT
+    import astropy.units as u
+    import numpy as np
+
+    # wolfram alpha result https://www.wolframalpha.com/input/?i=int_100%5E1000+x%5E-2
+    result = power_law_integral(
+        flux_normalization=1 * FLUX_UNIT,
+        spectral_index=-2,
+        e_min=100 * u.GeV,
+        e_max=1000 * u.GeV,
+        e_ref=1 * u.GeV,
+    )
+    assert np.isclose(result.value, 0.009)
+    assert result.unit == (FLUX_UNIT * u.GeV)