Skip to content

Commit a82e725

Browse files
maxnoemaxnoe
authored
Merge pull request #13 from fact-project/irf_writing
Write IRFs
2 parents bc36aed + 03a1fdc commit a82e725

21 files changed

+1500
-133
lines changed

‎README.md

+30-2
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,31 @@
1-
# irf
1+
# FACT irf
2+
Tools to calculate IRFs (instrument response functions) for the FACT telescope.
23

3-
Tools to calculate IRFs.
4+
Install pyeventio first `pip install eventio==0.3.0`.
5+
6+
# Joint Crab Data and the FACT Open Data Release
7+
8+
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)
9+
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
10+
cherenkov telescope protoyping silicon-photomultipliers.
11+
FACT data is stored in FITS files from the raw data level up to higher levels containing the results of our analysis.
12+
All of the data used here is publicly available, including the simulated data used for gauging the instrument response. For more information about the data visit [https://fact-project.org/data](https://fact-project.org/data) or read the corresponding section in the paper.
13+
14+
All of FACT's software is open source and can be accessed at [https://github.com/fact-project/](https://github.com/fact-project/).
15+
16+
For more information about FACT's data do not hesitate to contact
17+
18+
* Max Nöthe ([email protected])
19+
* Kai Brügge ([email protected])
20+
21+
Or contact the fact-online mailing list at
22+
23+
* [email protected]
24+
25+
For the joint crab publication we manually selected a smaller dataset
26+
due to some runs having bad weather.
27+
28+
```
29+
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
30+
31+
```

‎irf/__init__.py

+10-1
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,13 @@
11
from .collection_area import collection_area
2+
from .energy_dispersion import energy_dispersion, energy_migration
23
from .exposure_map import estimate_exposure_time, build_exposure_map
34

4-
__all__ = ['collection_area', 'estimate_exposure_time', 'build_exposure_map']
5+
__all__ = [
6+
'collection_area',
7+
'collection_area_to_irf_table',
8+
'energy_dispersion',
9+
'energy_dispersion_to_irf_table',
10+
'energy_migration',
11+
'estimate_exposure_time',
12+
'build_exposure_map'
13+
]

‎irf/collection_area.py

+16-20
Original file line numberDiff line numberDiff line change
@@ -1,15 +1,17 @@
11
import numpy as np
2+
23
from astropy.stats import binom_conf_interval
34
import astropy.units as u
45

6+
from scipy.ndimage.filters import gaussian_filter
7+
58

69
def histograms(
710
all_events,
811
selected_events,
912
bins,
1013
range=None,
11-
log=True,
12-
):
14+
):
1315
'''
1416
Create histograms in the given bins for two vectors.
1517
@@ -21,22 +23,12 @@ def histograms(
2123
Quantity which should be histogrammed for all selected events
2224
bins: int or array-like
2325
either number of bins or bin edges for the histogram
24-
range: (float, float)
25-
The lower and upper range of the bins
26-
log: bool
27-
flag indicating whether log10 should be applied to the values.
28-
2926
returns: hist_all, hist_selected, bin_edges
3027
'''
3128

32-
if log is True:
33-
all_events = np.log10(all_events)
34-
selected_events = np.log10(selected_events)
35-
3629
hist_all, bin_edges = np.histogram(
3730
all_events,
3831
bins=bins,
39-
range=range,
4032
)
4133

4234
hist_selected, _ = np.histogram(
@@ -53,10 +45,9 @@ def collection_area(
5345
selected_events,
5446
impact,
5547
bins,
56-
range=None,
57-
log=True,
5848
sample_fraction=1.0,
59-
):
49+
smoothing=0,
50+
):
6051
'''
6152
Calculate the collection area for the given events.
6253
@@ -70,19 +61,20 @@ def collection_area(
7061
either number of bins or bin edges for the histogram
7162
impact: astropy Quantity of type length
7263
The maximal simulated impact parameter
73-
log: bool
74-
flag indicating whether log10 should be applied to the quantity.
7564
sample_fraction: float
7665
The fraction of `all_events` that was analysed
7766
to create `selected_events`
67+
sample_fraction: float
68+
The fraction of `all_events` that was analysed
69+
to create `selected_events`
70+
smoothing: float
71+
The amount of smoothing to apply to the resulting matrix
7872
'''
7973

8074
hist_all, hist_selected, bin_edges = histograms(
8175
all_events,
8276
selected_events,
8377
bins,
84-
range=range,
85-
log=log
8678
)
8779

8880
hist_selected = (hist_selected / sample_fraction).astype(int)
@@ -99,6 +91,10 @@ def collection_area(
9991
lower_conf = lower_conf * np.pi * impact**2
10092
upper_conf = upper_conf * np.pi * impact**2
10193

102-
area = hist_selected / hist_all * np.pi * impact**2
94+
area = (hist_selected / hist_all) * np.pi * impact**2
95+
96+
if smoothing > 0:
97+
a = area.copy()
98+
area = gaussian_filter(a.value, sigma=smoothing, ) * area.unit
10399

104100
return area, bin_center, bin_width, lower_conf, upper_conf

‎irf/energy_dispersion.py

+134
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,134 @@
1+
import astropy.units as u
2+
import numpy as np
3+
from scipy.ndimage.filters import gaussian_filter
4+
5+
6+
7+
def _make_energy_bins(energy_true, energy_prediction, bins):
8+
e_min = min(
9+
min(energy_true),
10+
min(energy_prediction)
11+
)
12+
13+
e_max = max(
14+
max(energy_true),
15+
max(energy_prediction)
16+
)
17+
18+
low = np.log10(e_min.value)
19+
high = np.log10(e_max.value)
20+
bin_edges = np.logspace(low, high, endpoint=True, num=bins + 1)
21+
22+
return bin_edges
23+
24+
25+
26+
@u.quantity_input(energy_true=u.TeV, energy_prediction=u.TeV)
27+
def energy_dispersion(energy_true, energy_prediction, bins=10, normalize=False, smoothing=0):
28+
'''
29+
Creates energy dispersion matrix i.e. a histogram of e_reco vs e_true.
30+
31+
Parameters
32+
----------
33+
34+
energy_true : astropy.unit.Quantity (TeV)
35+
the true event energy
36+
energy_prediction: astropy.unit.Quantity (TeV)
37+
the predicted event energy
38+
bins_energy : int or arraylike
39+
the energy bins.
40+
normalize : bool
41+
Whether to normalize the matrix. The sum of each column will be equal to 1
42+
smoothing : float
43+
Amount of smoothing to apply to the generated matrices.
44+
Equivalent to the sigma parameter in
45+
https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.gaussian_filter.html
46+
Returns
47+
-------
48+
49+
array
50+
the migration matrix
51+
astropy.unit.Quantity
52+
the bin edges for the true energy axis
53+
astropy.unit.Quantity
54+
the bin edges for the predicted energy axis
55+
56+
'''
57+
if np.isscalar(bins):
58+
bins = _make_energy_bins(energy_true, energy_prediction, bins)
59+
60+
hist, bins_e_true, bins_e_prediction = np.histogram2d(
61+
energy_true.value,
62+
energy_prediction,
63+
bins=bins,
64+
)
65+
66+
if smoothing > 0:
67+
hist = gaussian_filter(hist, sigma=smoothing)
68+
69+
if normalize:
70+
with np.errstate(invalid='ignore'):
71+
h = hist.T
72+
h = h / h.sum(axis=0)
73+
hist = np.nan_to_num(h).T
74+
75+
return hist, bins_e_true * energy_true.unit, bins_e_prediction * energy_true.unit
76+
77+
78+
@u.quantity_input(energy_true=u.TeV, energy_prediction=u.TeV)
79+
def energy_migration(energy_true, energy_prediction, bins_energy=10, bins_mu=10, normalize=True, smoothing=0):
80+
'''
81+
Creates energy migration matrix i.e. a histogram of e_reco/e_true vs e_trueself.
82+
83+
Parameters
84+
----------
85+
86+
energy_true : astropy.unit.Quantity (TeV)
87+
the true event energy
88+
energy_prediction: astropy.unit.Quantity (TeV)
89+
the predicted event energy
90+
bins_energy : int or arraylike
91+
the energy bins.
92+
bins_mu : int or arraylike
93+
the bins to use for the y axis.
94+
normalize : bool
95+
Whether to normalize the matrix. The sum of each column will be equal to 1
96+
smoothing : float
97+
Amount of smoothing to apply to the generated matrices.
98+
Equivalent to the sigma parameter in
99+
https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.gaussian_filter.html
100+
Returns
101+
-------
102+
103+
array
104+
the migration matrix
105+
astropy.unit.Quantity
106+
the bin edges for the energy axis
107+
array
108+
the bin edges for the migration axis
109+
110+
'''
111+
if np.isscalar(bins_energy):
112+
bins_energy = _make_energy_bins(energy_true, energy_prediction, bins_energy)
113+
114+
migra = (energy_prediction / energy_true).si.value
115+
116+
if np.isscalar(bins_mu):
117+
bins_mu = np.linspace(0, 6, endpoint=True, num=bins_mu + 1)
118+
119+
hist, bins_e_true, bins_mu = np.histogram2d(
120+
energy_true.value,
121+
migra,
122+
bins=[bins_energy, bins_mu],
123+
)
124+
125+
if smoothing > 0:
126+
hist = gaussian_filter(hist, sigma=smoothing, )
127+
128+
if normalize:
129+
with np.errstate(invalid='ignore'):
130+
h = hist.T
131+
h = h / h.sum(axis=0)
132+
hist = np.nan_to_num(h).T
133+
134+
return hist, bins_energy * energy_true.unit, bins_mu

‎irf/gadf/__init__.py

+17
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,17 @@
1+
'''
2+
This module contains some helpful methods tp create insturmetn response functions and
3+
eventlists according to the open gamma-ray astro data format (gadf)
4+
5+
See https://gamma-astro-data-formats.readthedocs.io/en/latest/ for some specs.
6+
7+
'''
8+
9+
from . import time
10+
from . import response
11+
from . import hdus
12+
13+
__all__ = [
14+
'time',
15+
'response',
16+
'hdus'
17+
]

0 commit comments

Comments
 (0)