Merge pull request #22 from fact-project/fixes

Improve units handling in gadf

Author: maxnoe

README.md

@@ -1,11 +1,9 @@
 # FACT irf
 Tools to calculate IRFs (instrument response functions) for the FACT telescope.
 
-Install pyeventio first `pip install eventio==0.3.0`.
-
 # Joint Crab Data and the FACT Open Data Release
 
-The FACT open data release contains 10.6 hours of Crab Nebula observations taken in October 2013 by the First G-APD Cherenkov Telescope (FACT)
+The FACT open data release contains 17.7 hours of Crab Nebula observations taken in October 2013 by the First G-APD Cherenkov Telescope (FACT)
 FACT is located on the Roque des los Muchachos on the Canary island of La Palma off the west coast of Africa. It is an imaging atmospheric
 cherenkov telescope protoyping silicon-photomultipliers.
 FACT data is stored in FITS files from the raw data level up to higher levels containing the results of our analysis.
@@ -23,9 +21,10 @@ Or contact the fact-online mailing list at
  * [email protected]
 
  For the joint crab publication we manually selected a smaller dataset 
- due to some runs having bad weather.
+ due to some runs having bad weather and getting a dataset with roughly equal 
+ number of excess events for each experiment.
 
 ```
 python irf/scripts/fact_dl3_to_gadf.py -c 0.8 -t 0.03 ../gamma_diffuse_showers.hdf5 ../gamma_diffuse_precuts.hdf5 ../crab_dl3.hdf5 fact_dl3 --start '2013/11/04 00:00' --end '2013/11/07' -e 20131104_196 -e 20131105_175 -e 20131105_176 -e 20131105_199 -e 20131105_201
 
-```
+```

irf/collection_area.py

@@ -10,7 +10,6 @@ def histograms(
         all_events,
         selected_events,
         bins,
-        range=None,
 ):
     '''
     Create histograms in the given bins for two vectors.
@@ -64,9 +63,6 @@ def collection_area(
     sample_fraction: float
         The fraction of `all_events` that was analysed
         to create `selected_events`
-    sample_fraction: float
-        The fraction of `all_events` that was analysed
-        to create `selected_events`
     smoothing: float
         The amount of smoothing to apply to the resulting matrix
     '''
@@ -94,7 +90,6 @@ def collection_area(
     area = (hist_selected / hist_all) * np.pi * impact**2
 
     if smoothing > 0:
-        a = area.copy()
-        area = gaussian_filter(a.value, sigma=smoothing, ) * area.unit
+        area = gaussian_filter(area.value, sigma=smoothing) * area.unit
 
     return area, bin_center, bin_width, lower_conf, upper_conf

irf/energy_dispersion.py

@@ -3,28 +3,27 @@
 from scipy.ndimage.filters import gaussian_filter
 
 
-
-def _make_energy_bins(energy_true, energy_prediction, bins):
+@u.quantity_input(energy_true=u.TeV, energy_prediction=u.TeV)
+def _make_energy_bins(energy_true, energy_prediction, bins, e_ref=1 * u.TeV):
     e_min = min(
-        min(energy_true),
-        min(energy_prediction)
+        np.min(energy_true),
+        np.min(energy_prediction)
     )
 
     e_max = max(
-        max(energy_true),
-        max(energy_prediction)
+        np.max(energy_true),
+        np.max(energy_prediction)
     )
 
-    low = np.log10(e_min.value)
-    high = np.log10(e_max.value)
+    low = np.log10(e_min / e_ref)
+    high = np.log10(e_max / e_ref)
     bin_edges = np.logspace(low, high, endpoint=True, num=bins + 1)
 
-    return bin_edges
-
+    return bin_edges * e_ref
 
 
 @u.quantity_input(energy_true=u.TeV, energy_prediction=u.TeV)
-def energy_dispersion(energy_true, energy_prediction, bins=10, normalize=False, smoothing=0):
+def energy_dispersion(energy_true, energy_prediction, bins=10, normalize=False, smoothing=0, e_ref=1 * u.TeV):
     '''
     Creates energy dispersion matrix i.e. a histogram of e_reco vs e_true.
 
@@ -43,6 +42,8 @@ def energy_dispersion(energy_true, energy_prediction, bins=10, normalize=False,
         Amount of smoothing to apply to the generated matrices.
         Equivalent to the sigma parameter in
         https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.gaussian_filter.html
+    e_ref: astropy.unit.Quantity[energy]
+        Reference energy
     Returns
     -------
 
@@ -55,28 +56,25 @@ def energy_dispersion(energy_true, energy_prediction, bins=10, normalize=False,
 
     '''
     if np.isscalar(bins):
-        bins = _make_energy_bins(energy_true, energy_prediction, bins)
+        bins = _make_energy_bins(energy_true, energy_prediction, bins, e_ref=e_ref)
 
     hist, bins_e_true, bins_e_prediction = np.histogram2d(
-        energy_true.value,
-        energy_prediction,
-        bins=bins,
+        (energy_true / e_ref).to_value(u.dimensionless_unscaled),
+        (energy_prediction / e_ref).to_value(u.dimensionless_unscaled),
+        bins=(bins / e_ref).to_value(u.dimensionless_unscaled),
     )
 
     if smoothing > 0:
         hist = gaussian_filter(hist, sigma=smoothing)
 
     if normalize:
-        with np.errstate(invalid='ignore'):
-            h = hist.T
-            h = h / h.sum(axis=0)
-            hist = np.nan_to_num(h).T
+        hist = _normalize_hist(hist)
 
-    return hist, bins_e_true * energy_true.unit, bins_e_prediction * energy_true.unit
+    return hist, bins_e_true, bins_e_prediction
 
 
 @u.quantity_input(energy_true=u.TeV, energy_prediction=u.TeV)
-def energy_migration(energy_true, energy_prediction, bins_energy=10, bins_mu=10, normalize=True, smoothing=0):
+def energy_migration(energy_true, energy_prediction, bins_energy=10, bins_mu=10, normalize=True, smoothing=0, e_ref=1 * u.TeV):
     '''
     Creates energy migration matrix i.e. a histogram of e_reco/e_true vs e_trueself.
 
@@ -97,6 +95,8 @@ def energy_migration(energy_true, energy_prediction, bins_energy=10, bins_mu=10,
         Amount of smoothing to apply to the generated matrices.
         Equivalent to the sigma parameter in
         https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.gaussian_filter.html
+    e_ref: astropy.unit.Quantity[energy]
+        Reference energy
     Returns
     -------
 
@@ -109,26 +109,32 @@ def energy_migration(energy_true, energy_prediction, bins_energy=10, bins_mu=10,
 
     '''
     if np.isscalar(bins_energy):
-        bins_energy = _make_energy_bins(energy_true, energy_prediction, bins_energy)
+        bins_energy = _make_energy_bins(
+            energy_true, energy_prediction, bins_energy, e_ref=e_ref
+        )
 
-    migra = (energy_prediction / energy_true).si.value
+    migra = (energy_prediction / energy_true).to_value(u.dimensionless_unscaled)
 
     if np.isscalar(bins_mu):
-        bins_mu = np.linspace(0, 6, endpoint=True, num=bins_mu + 1)
+        bins_mu = np.linspace(0, 6, bins_mu + 1)
 
     hist, bins_e_true, bins_mu = np.histogram2d(
-        energy_true.value,
+        (energy_true / e_ref).to_value(u.dimensionless_unscaled),
         migra,
-        bins=[bins_energy, bins_mu],
+        bins=[(bins_energy / e_ref).to_value(u.dimensionless_unscaled), bins_mu],
     )
 
     if smoothing > 0:
-        hist = gaussian_filter(hist, sigma=smoothing, )
+        hist = gaussian_filter(hist, sigma=smoothing)
 
     if normalize:
-        with np.errstate(invalid='ignore'):
-            h = hist.T
-            h = h / h.sum(axis=0)
-            hist = np.nan_to_num(h).T
+        hist = _normalize_hist(hist)
+
+    return hist, bins_energy, bins_mu
+
 
-    return hist, bins_energy * energy_true.unit, bins_mu
+def _normalize_hist(hist):
+    with np.errstate(invalid='ignore'):
+        h = hist.T
+        h = h / h.sum(axis=0)
+        return np.nan_to_num(h).T

irf/gadf/hdus.py

@@ -11,7 +11,6 @@
 from astropy.table import Table
 from astropy.time import Time
 import numpy as np
-import datetime
 from irf.gadf.time import timestamp_to_mjdref, ontime_info_from_runs, TIME_INFO
 
 
@@ -41,12 +40,11 @@ def create_primary_hdu():
     header['OBSERVER'] = 'The FACT collaboration'
     header['COMMENT'] = 'FACT OGA.'
     header['COMMENT'] = 'See https://gamma-astro-data-formats.readthedocs.io/en/latest/'
-    header['COMMENT'] = 'This file was created by https://github.com/fact-project/irf'
-    header['COMMENT'] = 'See our full analysis on GitHub'
+    header['COMMENT'] = 'See our Open Crab analysis on GitHub'
     header['COMMENT'] = 'https://github.com/fact-project/open_crab_sample_analysis'
 
-    now = Time.now().iso
-    header['COMMENT'] = f'This file was created on {now}'
+    header['CREATOR'] = 'This file was created by https://github.com/fact-project/irf'
+    header['DATE'] = Time.now().iso
 
     return fits.PrimaryHDU(header=header)
 
@@ -63,17 +61,23 @@ def create_events_hdu(dl3, run):
     event_id = dl3.index
 
     # convert from hour angles to deg
-    ra = dl3.ra_prediction.values * u.hourangle
+    ra = u.Quantity(dl3.ra_prediction.values, u.hourangle, copy=False)
     ra = ra.to('deg')
 
-    dec = dl3.dec_prediction.values * u.deg
+    dec = u.Quantity(dl3.dec_prediction.values, u.deg, copy=False)
 
     # convert to TeV energy
-    energy = dl3.gamma_energy_prediction.values * u.GeV
+    energy = u.Quantity(dl3.gamma_energy_prediction.values, u.GeV, copy=False)
     energy = energy.to('TeV')
 
     timestamp = timestamp_to_mjdref(dl3.timestamp)
-    t = Table({'EVENT_ID': event_id, 'ENERGY': energy, 'DEC': dec, 'RA': ra, 'TIME': timestamp})
+    t = Table({
+        'EVENT_ID': event_id,
+        'ENERGY': energy,
+        'DEC': dec,
+        'RA': ra,
+        'TIME': timestamp
+    })
     hdu = fits.table_to_hdu(t)
 
     # add information to HDU header
@@ -146,14 +150,13 @@ def create_observation_index_hdu(runs):
     returns a fits hdu object
     '''
     obs_id = _observation_ids(runs)
-    ra_pnt = runs.right_ascension.values * u.hourangle
+    ra_pnt = u.Quantity(runs.right_ascension.values, u.hourangle, copy=False)
     ra_pnt = ra_pnt.to('deg')
 
-    dec_pnt = runs.declination.values * u.deg
-    zen_pnt = runs.zenith.values * u.deg
-    alt_pnt = (90 - runs.zenith.values) * u.deg
-    az_pnt = runs.azimuth.values * u.deg
-
+    dec_pnt = u.Quantity(runs.declination.values, u.deg, copy=False)
+    zen_pnt = u.Quantity(runs.zenith.values, u.deg, copy=False)
+    alt_pnt = u.Quantity(90 - runs.zenith.values, u.deg, copy=False)
+    az_pnt = u.Quantity(runs.azimuth.values, u.deg, copy=False)
 
     n_tels = np.ones_like(obs_id)
     tellist = np.ones_like(obs_id).astype(np.str)
@@ -198,7 +201,7 @@ def add_meta_information_to_hdu(hdu, **kwargs):
     hdu.header['EQUINOX'] = '2000.0'
     hdu.header['RADECSYS'] = 'ICRS'
     hdu.header['EUNIT'] = 'TeV'
-    hdu.header['DATE'] = datetime.datetime.now().replace(microsecond=0).isoformat()
+    hdu.header['DATE'] = Time.now().iso
     _extend_hdu_header(hdu.header, TIME_INFO)
     if kwargs:
         _extend_hdu_header(hdu.header, kwargs)

setup.py

@@ -3,7 +3,7 @@
 
 setup(
     name='irf',
-    version='0.4.0',
+    version='0.4.1',
     description='Functions to do instrument response functions for FACT',
     url='http://github.com/fact-project/irf',
     author='Kai Brügge, Maximilian Nöthe',