Skip to content

Commit

Permalink
Merge pull request #109 from fact-project/units
Browse files Browse the repository at this point in the history
Use astropy units consistently, fix spectral_index madness
  • Loading branch information
maxnoe authored Jun 15, 2018
2 parents 13c6619 + 59df237 commit 1a83239
Show file tree
Hide file tree
Showing 6 changed files with 176 additions and 65 deletions.
2 changes: 1 addition & 1 deletion fact/VERSION
Original file line number Diff line number Diff line change
@@ -1 +1 @@
0.19.1
0.20.0
196 changes: 145 additions & 51 deletions fact/analysis/statistics.py
Original file line number Diff line number Diff line change
@@ -1,65 +1,92 @@
import numpy as np
import astropy.units as u

POINT_SOURCE_FLUX_UNIT = (1 / u.GeV / u.s / u.m**2).unit
FLUX_UNIT = POINT_SOURCE_FLUX_UNIT / u.sr

def random_power(spectral_index, e_min, e_max, size):

@u.quantity_input
def random_power(
spectral_index,
e_min: u.TeV,
e_max: u.TeV,
size,
e_ref: u.TeV=1 * u.TeV,
) -> u.TeV:
r'''
Draw random numbers from a power law distribution

.. math::
f(E) =
\frac{1 - \gamma}
{E_{\max}^{1 - \gamma} - E_{\min}^{1 - \gamma}}
E^{-\gamma}
({E_{\max} / E_{\mathrm{ref}})^{\gamma - 1} - (E_{\min} / E_{\mathrm{ref}})^{\gamma - 1}}
(E / E_{\mathrm{ref}})^{\gamma}

Parameters
----------
spectral_index: float
The differential spectral index of the power law
e_min: float
e_min: Quantity[energy]
lower energy border
e_max: float
e_max: Quantity[energy]
upper energy border, can be np.inf
size: int or tuple[int]
Number of events to draw or a shape like (100, 2)
'''
assert spectral_index > 1.0, 'spectral_index must be > 1.0'
if spectral_index > -1.0:
raise ValueError('spectral_index must be < -1.0')

u = np.random.uniform(0, 1, size)

exponent = 1 - spectral_index
exponent = spectral_index + 1

if e_max == np.inf:
diff = -e_min**exponent
diff = -(e_min / e_ref)**exponent
else:
diff = e_max**exponent - e_min**exponent
diff = (e_max / e_ref)**exponent - (e_min / e_ref)**exponent

return (diff * u + e_min**exponent) ** (1 / exponent)
return e_ref * (diff * u + (e_min / e_ref)**exponent) ** (1 / exponent)


def power_law(energy, flux_normalization, spectral_index):
def power_law(
energy: u.TeV,
flux_normalization: (FLUX_UNIT, POINT_SOURCE_FLUX_UNIT),
spectral_index,
e_ref: u.TeV=1 * u.TeV,
):
r'''
Simple power law

.. math::
\phi = \phi_0 \cdot E ^{-\gamma}
\phi = \phi_0 \cdot (E / E_{}\mathrm{ref}})^{\gamma}

Parameters
----------
energy: number or array-like
energy: Quantity[energy]
energy points to evaluate
flux_normalization: float
flux_normalization: Quantity[m**-2 s**-2 TeV**-1]
Flux normalization
spectral_index: float
Spectral index
e_ref: Quantity[energy]
The reference energy
'''
return flux_normalization * energy**(-spectral_index)
return flux_normalization * (energy / e_ref)**(spectral_index)


def curved_power_law(energy, flux_normalization, a, b):
@u.quantity_input
def curved_power_law(
energy: u.TeV,
flux_normalization: (FLUX_UNIT, POINT_SOURCE_FLUX_UNIT),
a,
b,
e_ref: u.TeV=1 * u.TeV,
):
r'''
Curved power law

.. math::
\phi = \phi_0 \cdot E ^ {a - b \cdot \log(E)}
\phi = \phi_0 \cdot E ^ {a + b \cdot \log_{10}(E)}

Parameters
----------
Expand All @@ -70,9 +97,12 @@ def curved_power_law(energy, flux_normalization, a, b):
a: float
Parameter `a` of the curved power law
b: float
Parameter `b` of the curved power law
Parameter `b` of the curved power law
e_ref: Quantity[energy]
The reference energy
'''
return flux_normalization * energy ** (a + b * np.log10(energy))
exp = a + b * np.log10((energy / e_ref))
return flux_normalization * (energy / e_ref) ** exp


def li_ma_significance(n_on, n_off, alpha=0.2):
Expand Down Expand Up @@ -115,7 +145,14 @@ def li_ma_significance(n_on, n_off, alpha=0.2):
return significance


def calc_weight_simple_to_curved(energy, spectral_index, a, b):
@u.quantity_input
def calc_weight_simple_to_curved(
energy: u.TeV,
spectral_index,
a,
b,
e_ref: u.TeV = 1 * u.TeV,
):
'''
Reweight simulated events from a simulated power law with
spectral index `spectral_index` to a curved power law with parameters `a` and `b`
Expand All @@ -131,10 +168,16 @@ def calc_weight_simple_to_curved(energy, spectral_index, a, b):
b: float
Parameter `b` of the target curved power law
'''
return energy ** (spectral_index + a + b * np.log10(energy))
return (energy / e_ref) ** (spectral_index + a + b * np.log10(energy / e_ref))


def calc_weight_change_index(energy, simulated_index, target_index):
@u.quantity_input
def calc_weight_change_index(
energy: u.TeV,
simulated_index,
target_index,
e_ref: u.TeV=1 * u.TeV,
):
'''
Reweight simulated events from one power law index to another

Expand All @@ -147,60 +190,109 @@ def calc_weight_change_index(energy, simulated_index, target_index):
target_index: float
Spectral index of the target power law
'''
return energy ** (target_index - simulated_index)
return (energy / e_ref) ** (target_index - simulated_index)


@u.quantity_input
def calc_gamma_obstime(
n_events,
spectral_index,
flux_normalization: (FLUX_UNIT, POINT_SOURCE_FLUX_UNIT),
max_impact: u.m,
e_min: u.TeV,
e_max: u.TeV,
e_ref: u.TeV = 1 * u.TeV,
) -> u.s:
r'''
Calculate the equivalent observation time for a number of simulation enents

The number of events produced by sampling from a power law with
spectral index γ is given by

def calc_gamma_obstime(n_events, spectral_index, flux_normalization, max_impact, e_min, e_max):
'''
Calculate the equivalent observation time for a spectral_index montecarlo set
.. math::
N =
A t \Phi_0 \int_{E_{\min}}^{E_{\max}}
\left(\frac{E}{E_{\mathrm{ref}}}\right)^{\gamma}
\,\mathrm{d}E

Solving this for t yields

.. math::
t = \frac{N \cdot (\gamma - 1)}{
E_{\mathrm{ref}} A \Phi_0
\left(
\left(\frac{E_{\max}}{E_{\mathrm{ref}}}\right)^{\gamma - 1}
- \left(\frac{E_{\max}}{E_{\mathrm{ref}}}\right)^{\gamma - 1}
\right)
}

Parameters
----------
n_events: int
Number of simulated events
spectral_index: float
Spectral index of the simulated power law
Spectral index of the simulated power law, including the sign,
so typically -2.7 or -2
flux_normalization: float
Flux normalization of the simulated power law
max_impact: float
max_impact: Quantity[length]
Maximal simulated impact
e_min: float
e_min: Quantity[energy]
Mimimal simulated energy
e_max: float
e_max: Quantity[energy]
Maximal simulated energy
e_ref: Quantity[energy]
The e_ref energy
'''
if spectral_index >= -1:
raise ValueError('spectral_index must be < -1')

numerator = n_events * (1 - spectral_index)
numerator = n_events * (spectral_index + 1)

t1 = flux_normalization * max_impact**2 * np.pi
t2 = e_max**(1 - spectral_index) - e_min**(1 - spectral_index)
A = max_impact**2 * np.pi
t1 = A * flux_normalization * e_ref
t2 = (e_max / e_ref)**(spectral_index + 1)
t3 = (e_min / e_ref)**(spectral_index + 1)

denominator = t1 * t2
denominator = t1 * (t2 - t3)

return numerator / denominator


def power_law_integral(flux_normalisation, spectral_index, e_min, e_max):
@u.quantity_input
def power_law_integral(
flux_normalization: (FLUX_UNIT, POINT_SOURCE_FLUX_UNIT),
spectral_index,
e_min: u.TeV,
e_max: u.TeV,
e_ref: u.TeV = 1 * u.TeV,
):
'''
Return the integral of a power_law with normalisation
`flux_normalization` and index `spectral_index`
between `e_min` and `e_max`
'''
int_index = 1 - spectral_index
return flux_normalisation / int_index * (e_max**(int_index) - e_min**(int_index))
int_index = spectral_index + 1
e_term = (e_max / e_ref)**(int_index) - (e_min / e_ref)**(int_index)

res = flux_normalization * e_ref / int_index * e_term

if flux_normalization.unit.is_equivalent(FLUX_UNIT):
return res.to(1 / u.m**2 / u.s / u.sr)
return res.to(1 / u.m**2 / u.s)


@u.quantity_input
def calc_proton_obstime(
n_events,
spectral_index,
max_impact,
viewcone,
e_min,
e_max,
flux_normalization=1.8e4,
):
n_events,
spectral_index,
max_impact: u.m,
viewcone: u.deg,
e_min: u.TeV,
e_max: u.TeV,
flux_normalization: FLUX_UNIT = 1.8e4 * FLUX_UNIT,
e_ref: u.GeV=1 * u.GeV,
) -> u.s:
'''
Calculate the equivalent observation time for a proton montecarlo set

Expand All @@ -214,18 +306,20 @@ def calc_proton_obstime(
Maximal simulated impact in m
viewcone: float
Viewcone in degrees
e_min: float
Mimimal simulated energy in GeV
e_max: float
Maximal simulated energy in GeV
e_min: Quantity[energy]
Mimimal simulated energy
e_max: Quantity[energy]
Maximal simulated energy
flux_normalization: float
Flux normalisation of the cosmic rays in nucleons / (m² s sr GeV)
Default value is (29.2) of
http://pdg.lbl.gov/2016/reviews/rpp2016-rev-cosmic-rays.pdf
'''
area = np.pi * max_impact**2
solid_angle = 2 * np.pi * (1 - np.cos(np.deg2rad(viewcone)))
solid_angle = 2 * np.pi * (1 - np.cos(viewcone)) * u.sr

expected_integral_flux = power_law_integral(flux_normalization, 2.7, e_min, e_max)
expected_integral_flux = power_law_integral(
flux_normalization, spectral_index, e_min, e_max, e_ref
)

return n_events / area / solid_angle / expected_integral_flux
1 change: 1 addition & 0 deletions fact/io.py
Original file line number Diff line number Diff line change
Expand Up @@ -149,6 +149,7 @@ def read_h5py(

if 'index' in df.columns:
df.set_index('index', inplace=True)
df.index.name = None

return df

Expand Down
2 changes: 1 addition & 1 deletion fact/time.py
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,7 @@ def datetime_to_mjd(dt):

elif isinstance(dt, (pd.Series, pd.DatetimeIndex)):
if dt.tz is None:
dt.tz = timezone.utc
dt = dt.tz_localize(timezone.utc)

return (dt - MJD_EPOCH).total_seconds() / 24 / 3600

Expand Down
28 changes: 22 additions & 6 deletions tests/test_analysis.py
Original file line number Diff line number Diff line change
@@ -1,12 +1,28 @@
import astropy.units as u
from pytest import approx, raises


def test_proton_obstime():
from fact.analysis.statistics import calc_proton_obstime
n_simulated = 780046520
t = calc_proton_obstime(
n_events=n_simulated,
spectral_index=2.7,
max_impact=400,
viewcone=5,
e_min=100,
e_max=200e3,
spectral_index=-2.7,
max_impact=400 * u.m,
viewcone=5 * u.deg,
e_min=100 * u.GeV,
e_max=200 * u.TeV,
)
assert int(t) == 15397
assert t.to(u.s).value == approx(15397.82)


def test_power():
from fact.analysis.statistics import random_power

a = random_power(-2.7, e_min=5 * u.GeV, e_max=10 * u.TeV, size=1000)

assert a.shape == (1000, )
assert a.unit == u.TeV

with raises(ValueError):
random_power(2.7, 5 * u.GeV, 10 * u.GeV, e_ref=1 * u.GeV, size=1)
Loading

0 comments on commit 1a83239

Please sign in to comment.