Merge pull request #108 from cta-observatory/add_particle_types

Add counts / weighted counts per particle type

docs/notebooks/comparison_with_EventDisplay.ipynb

@@ -268,13 +268,11 @@
     "sensitivity['reco_energy_low'].info.format = '.3g'\n",
     "sensitivity['reco_energy_high'].info.format = '.3g'\n",
     "sensitivity['reco_energy_center'].info.format = '.3g'\n",
-    "sensitivity['n_signal'].info.format = '.1f'\n",
-    "sensitivity['n_signal_weighted'].info.format = '.1f'\n",
-    "sensitivity['n_background_weighted'].info.format = '.1f'\n",
-    "sensitivity['n_background'].info.format = '.1f'\n",
     "sensitivity['relative_sensitivity'].info.format = '.2g'\n",
     "sensitivity['flux_sensitivity'].info.format = '.3g'\n",
     "\n",
+    "for k in filter(lambda k: k.startswith('n_'), sensitivity.colnames):\n",
+    "    sensitivity[k].info.format = '.1f'\n",
     "\n",
     "sensitivity"
    ]
@@ -314,8 +312,9 @@
     "\n",
     "\n",
     "e = sensitivity['reco_energy_center']\n",
-    "s = (e**2 * sensitivity['flux_sensitivity'])\n",
     "w = (sensitivity['reco_energy_high'] - sensitivity['reco_energy_low'])\n",
+    "s = (e**2 * sensitivity['flux_sensitivity'])\n",
+    "\n",
     "\n",
     "ax_sens.errorbar(\n",
     "    e.to_value(u.TeV),\n",
@@ -698,7 +697,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.7"
+   "version": "3.7.9"
   }
  },
  "nbformat": 4,

environment.yml

@@ -11,6 +11,7 @@ dependencies:
   - scipy
   - astropy>=4.0.2,<5
   - setuptools
+  - tqdm
   # tests
   - pytest
   - pytest-cov

examples/calculate_eventdisplay_irfs.py

@@ -87,6 +87,7 @@ def main():
     for k, p in particles.items():
         log.info(f"Simulated {k.title()} Events:")
         p["events"], p["simulation_info"] = read_eventdisplay_fits(p["file"])
+        p["events"]["particle_type"] = k
 
         p["simulated_spectrum"] = PowerLaw.from_simulation(p["simulation_info"], T_OBS)
         p["events"]["weight"] = calculate_event_weights(
@@ -283,7 +284,7 @@ def main():
         fov_offset_bins=np.arange(0, 11) * u.deg,
     ))
     hdus.append(create_psf_table_hdu(
-            psf, true_energy_bins, source_offset_bins, fov_offset_bins,
+        psf, true_energy_bins, source_offset_bins, fov_offset_bins,
     ))
     hdus.append(create_rad_max_hdu(
         theta_cuts_opt["cut"][:, np.newaxis], theta_bins, fov_offset_bins

pyirf/binning.py

@@ -132,6 +132,27 @@ def create_histogram_table(events, bins, key="reco_energy"):
         hist["n_weighted"], _ = np.histogram(
             events[key], bins, weights=events["weight"]
         )
+    else:
+        hist["n_weighted"] = hist["n"]
+
+    # create counts per particle type, only works if there is at least 1 event
+    if 'particle_type' in events.colnames and len(events) > 0:
+        by_particle = events.group_by('particle_type')
+
+        for group_key, group in zip(by_particle.groups.keys, by_particle.groups):
+            particle = group_key['particle_type']
+
+            hist["n_" + particle], _ = np.histogram(group[key], bins)
+
+            # also calculate weighted number of events
+            col = "n_" + particle
+            if "weight" in events.colnames:
+                hist[col + "_weighted"], _ = np.histogram(
+                    group[key], bins, weights=group["weight"]
+                )
+            else:
+                hist[col + "_weighted"] = hist[col]
+
     return hist
 
 

pyirf/sensitivity.py

@@ -181,10 +181,14 @@ def calculate_sensitivity(
         s[k] = signal_hist[k]
 
     s["n_signal"] = signal_hist["n"] * rel_sens
-    s["n_background"] = background_hist["n"]
     s["n_signal_weighted"] = signal_hist["n_weighted"] * rel_sens
+    s["n_background"] = background_hist["n"]
     s["n_background_weighted"] = background_hist["n_weighted"]
 
+    # copy also "n_proton" / "n_electron_weighted" etc. if available
+    for k in filter(lambda c: c.startswith('n_') and c != 'n_weighted', background_hist.colnames):
+        s[k] = background_hist[k]
+
     # perform safety checks
     # we use the number of signal events at the flux level that yields
     # the target significance
@@ -262,7 +266,7 @@ def estimate_background(
         / cone_solid_angle(background_radius)
     ).to_value(u.one)
 
-    for key in ['n', 'n_weighted']:
+    for key in filter(lambda col: col.startswith('n'), bg.colnames):
         # *= not possible due to upcast from int to float
         bg[key] = bg[key] * size_ratio / alpha
 

pyirf/tests/test_binning.py

@@ -1,4 +1,5 @@
 import astropy.units as u
+from astropy.table import QTable
 import numpy as np
 
 
@@ -56,6 +57,50 @@ def test_bins_per_decade():
     assert np.allclose(np.diff(np.log10(bins.to_value(u.GeV))), 0.1)
 
 
+def test_create_histogram_table():
+    '''Test create histogram table'''
+    from pyirf.binning import create_histogram_table
+
+    events = QTable({
+        'reco_energy': [1, 1, 10, 100, 100, 100] * u.TeV,
+    })
+    bins = [0.5, 5, 50, 500] * u.TeV
+
+    # test without weights
+    hist = create_histogram_table(events, bins, key='reco_energy')
+    assert np.all(hist['n'] == [2, 1, 3])
+    assert np.all(hist['n_weighted'] == [2, 1, 3])
+
+    # test with weights
+    events['weight'] = [0.5, 2, 2.5, 0.5, 0.2, 0.3]
+    hist = create_histogram_table(events, bins, key='reco_energy')
+    assert np.all(hist['n'] == [2, 1, 3])
+    assert np.all(hist['n_weighted'] == [2.5, 2.5, 1.0])
+
+    # test with particle_types
+    types = np.array(['gamma', 'electron', 'gamma', 'electron', 'gamma', 'proton'])
+    events['particle_type'] = types
+    hist = create_histogram_table(events, bins, key='reco_energy')
+
+    assert np.all(hist['n'] == [2, 1, 3])
+    assert np.all(hist['n_weighted'] == [2.5, 2.5, 1.0])
+
+    assert np.all(hist['n_gamma'] == [1, 1, 1])
+    assert np.all(hist['n_electron'] == [1, 0, 1])
+    assert np.all(hist['n_proton'] == [0, 0, 1])
+
+    assert np.allclose(hist['n_gamma_weighted'], [0.5, 2.5, 0.2])
+    assert np.allclose(hist['n_electron_weighted'], [2, 0, 0.5])
+    assert np.allclose(hist['n_proton_weighted'], [0, 0, 0.3])
+
+    # test with empty table
+    empty = events[:0]
+    hist = create_histogram_table(empty, bins, key='reco_energy')
+    zeros = np.zeros(3)
+    assert np.all(hist['n'] == zeros)
+    assert np.all(hist['n_weighted'] == zeros)
+
+
 def test_calculate_bin_indices():
     from pyirf.binning import calculate_bin_indices
 

pyirf/version.py

@@ -1 +1 @@
-__version__ = "0.3.1"
+__version__ = "0.4.0.dev1"

setup.cfg

@@ -28,3 +28,12 @@ requires-dist = setuptools
 
 [options]
 python_requires = >=3.6
+
+[flake8]
+exclude=
+  build,
+  docs,
+  .eggs
+max-line-length=90
+select = C,E,F,W,B,B950
+ignore = E501,W503,E203

setup.py

@@ -32,7 +32,6 @@
         "astropy~=4.0,>=4.0.2",
         "matplotlib",
         "numpy>=1.18",
-        "pandas",
         "scipy",
         "tqdm",
         "tables",