Merge branch 'master' into remove_mongo

Author: maxnoe

Diff for: fact/VERSION

@@ -1 +1 @@
-0.26.0
+0.26.0

Diff for: 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

Diff for: fact/io.py

@@ -546,10 +546,26 @@ def read_simulated_spectrum(corsika_headers_path):
     runs = read_h5py(corsika_headers_path, key='corsika_runs')
 
     summary = {}
-    try:
-        summary['n_showers'] = runs['n_showers'].sum()
-    except KeyError:
-        summary['n_showers'] = runs['n_events'].sum()
+    if 'n_showers' in runs.columns:
+        n_showers = runs['n_showers']
+    else:
+        n_showers = runs['n_events']
+
+    summary['n_showers'] = n_showers.sum()
+
+    with h5py.File(corsika_headers_path, 'r') as f:
+        summary['n_reuse'] = f['corsika_runs'].attrs.get('n_reuse', 1)
+
+        if 'n_reuse' in runs.columns:
+            # if reuse is not the same for all runs, multply n_showers
+            # and set reuse to 1
+            unique = runs['reuse'].unique()
+            if len(unique) > 1:
+                summary['n_showers'] = (n_showers * runs['n_reuse']).sum()
+                summary['n_reuse'] = 1
+
+            else:
+                summary['n_reuse'] = runs['n_reuse'].iloc[0]
 
     keys = {'energy_min': u.GeV, 'energy_max': u.GeV, 'energy_spectrum_slope': None}
     if 'x_scatter' in runs.columns:

Diff for: 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)