1 | # Licensed under a 3-clause BSD style license - see LICENSE.rst |
||
0 ignored issues
–
show
coding-style
introduced
by
Loading history...
|
|||
2 | """Spectral models for Gammapy.""" |
||
3 | from __future__ import absolute_import, division, print_function, unicode_literals |
||
4 | import numpy as np |
||
5 | import copy |
||
6 | import operator |
||
7 | import astropy.units as u |
||
8 | from astropy.table import Table |
||
9 | from ..utils.energy import EnergyBounds |
||
10 | from ..utils.nddata import NDDataArray, BinnedDataAxis |
||
11 | from .utils import integrate_spectrum |
||
12 | from ..utils.scripts import make_path |
||
13 | from ..utils.fitting import Parameter, Parameters |
||
14 | from ..utils.interpolation import ScaledRegularGridInterpolator |
||
15 | |||
16 | __all__ = [ |
||
17 | "SpectralModel", |
||
18 | "ConstantModel", |
||
19 | "CompoundSpectralModel", |
||
20 | "PowerLaw", |
||
21 | "PowerLaw2", |
||
22 | "ExponentialCutoffPowerLaw", |
||
23 | "ExponentialCutoffPowerLaw3FGL", |
||
24 | "PLSuperExpCutoff3FGL", |
||
25 | "LogParabola", |
||
26 | "TableModel", |
||
27 | "AbsorbedSpectralModel", |
||
28 | "Absorption", |
||
29 | ] |
||
30 | |||
31 | |||
32 | class SpectralModel(object): |
||
33 | """Spectral model base class. |
||
34 | |||
35 | Derived classes should store their parameters as |
||
36 | `~gammapy.utils.modeling.Parameters` |
||
37 | See for example return pardict of |
||
38 | `~gammapy.spectrum.models.PowerLaw`. |
||
39 | """ |
||
40 | |||
41 | def __repr__(self): |
||
42 | fmt = "{}()" |
||
43 | return fmt.format(self.__class__.__name__) |
||
44 | |||
45 | View Code Duplication | def __str__(self): |
|
46 | ss = self.__class__.__name__ |
||
47 | ss += "\n\nParameters: \n\n\t" |
||
48 | |||
49 | table = self.parameters.to_table() |
||
50 | ss += "\n\t".join(table.pformat()) |
||
51 | |||
52 | if self.parameters.covariance is not None: |
||
53 | ss += "\n\nCovariance: \n\n\t" |
||
54 | covar = self.parameters.covariance_to_table() |
||
55 | ss += "\n\t".join(covar.pformat()) |
||
56 | return ss |
||
57 | |||
58 | def __call__(self, energy): |
||
59 | """Call evaluate method of derived classes""" |
||
60 | kwargs = dict() |
||
61 | for par in self.parameters.parameters: |
||
62 | kwargs[par.name] = par.quantity |
||
63 | |||
64 | return self.evaluate(energy, **kwargs) |
||
65 | |||
66 | def __mul__(self, model): |
||
67 | if not isinstance(model, SpectralModel): |
||
68 | model = ConstantModel(const=model) |
||
69 | return CompoundSpectralModel(self, model, operator.mul) |
||
70 | |||
71 | def __rmul__(self, model): |
||
72 | # This is needed to support e.g. 5 * model |
||
73 | return self.__mul__(model) |
||
74 | |||
75 | def __add__(self, model): |
||
76 | if not isinstance(model, SpectralModel): |
||
77 | model = ConstantModel(const=model) |
||
78 | return CompoundSpectralModel(self, model, operator.add) |
||
79 | |||
80 | def __radd__(self, model): |
||
81 | return self.__add__(model) |
||
82 | |||
83 | def __sub__(self, model): |
||
84 | if not isinstance(model, SpectralModel): |
||
85 | model = ConstantModel(const=model) |
||
86 | return CompoundSpectralModel(self, model, operator.sub) |
||
87 | |||
88 | def __rsub__(self, model): |
||
89 | return self.__sub__(model) |
||
90 | |||
91 | def __truediv__(self, model): |
||
92 | if not isinstance(model, SpectralModel): |
||
93 | model = ConstantModel(const=model) |
||
94 | return CompoundSpectralModel(self, model, operator.truediv) |
||
95 | |||
96 | def __rtruediv__(self, model): |
||
97 | return self.__div__(model) |
||
98 | |||
99 | def _parse_uarray(self, uarray): |
||
100 | from uncertainties import unumpy |
||
101 | |||
102 | values = unumpy.nominal_values(uarray) |
||
103 | errors = unumpy.std_devs(uarray) |
||
104 | return values, errors |
||
105 | |||
106 | def _convert_energy(self, energy): |
||
107 | if "reference" in self.parameters.names: |
||
108 | return energy.to(self.parameters["reference"].unit) |
||
109 | elif "emin" in self.parameters.names: |
||
110 | return energy.to(self.parameters["emin"].unit) |
||
111 | else: |
||
112 | return energy |
||
113 | |||
114 | def evaluate_error(self, energy): |
||
115 | """Evaluate spectral model with error propagation. |
||
116 | |||
117 | Parameters |
||
118 | ---------- |
||
119 | energy : `~astropy.units.Quantity` |
||
120 | Energy at which to evaluate |
||
121 | |||
122 | Returns |
||
123 | ------- |
||
124 | flux, flux_error : tuple of `~astropy.units.Quantity` |
||
125 | Tuple of flux and flux error. |
||
126 | """ |
||
127 | energy = self._convert_energy(energy) |
||
128 | |||
129 | unit = self(energy).unit |
||
130 | upars = self.parameters._ufloats |
||
131 | uarray = self.evaluate(energy.value, **upars) |
||
132 | return self._parse_uarray(uarray) * unit |
||
133 | |||
134 | def integral(self, emin, emax, **kwargs): |
||
135 | """Integrate spectral model numerically. |
||
136 | |||
137 | .. math:: |
||
138 | |||
139 | F(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}}\phi(E)dE |
||
140 | |||
141 | If array input for ``emin`` and ``emax`` is given you have to set |
||
142 | ``intervals=True`` if you want the integral in each energy bin. |
||
143 | |||
144 | Parameters |
||
145 | ---------- |
||
146 | emin, emax : `~astropy.units.Quantity` |
||
147 | Lower and upper bound of integration range. |
||
148 | **kwargs : dict |
||
149 | Keyword arguments passed to :func:`~gammapy.spectrum.integrate_spectrum` |
||
150 | """ |
||
151 | return integrate_spectrum(self, emin, emax, **kwargs) |
||
152 | |||
153 | def integral_error(self, emin, emax, **kwargs): |
||
154 | """Integrate spectral model numerically with error propagation. |
||
155 | |||
156 | Parameters |
||
157 | ---------- |
||
158 | emin, emax : `~astropy.units.Quantity` |
||
159 | Lower adn upper bound of integration range. |
||
160 | **kwargs : dict |
||
161 | Keyword arguments passed to func:`~gammapy.spectrum.integrate_spectrum` |
||
162 | |||
163 | Returns |
||
164 | ------- |
||
165 | integral, integral_error : tuple of `~astropy.units.Quantity` |
||
166 | Tuple of integral flux and integral flux error. |
||
167 | """ |
||
168 | emin = self._convert_energy(emin) |
||
169 | emax = self._convert_energy(emax) |
||
170 | unit = self.integral(emin, emax, **kwargs).unit |
||
171 | upars = self.parameters._ufloats |
||
172 | |||
173 | def f(x): |
||
174 | return self.evaluate(x, **upars) |
||
175 | |||
176 | uarray = integrate_spectrum(f, emin.value, emax.value, **kwargs) |
||
177 | return self._parse_uarray(uarray) * unit |
||
178 | |||
179 | def energy_flux(self, emin, emax, **kwargs): |
||
180 | """Compute energy flux in given energy range. |
||
181 | |||
182 | .. math:: |
||
183 | |||
184 | G(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}}E \phi(E)dE |
||
185 | |||
186 | Parameters |
||
187 | ---------- |
||
188 | emin, emax : `~astropy.units.Quantity` |
||
189 | Lower and upper bound of integration range. |
||
190 | **kwargs : dict |
||
191 | Keyword arguments passed to func:`~gammapy.spectrum.integrate_spectrum` |
||
192 | """ |
||
193 | |||
194 | def f(x): |
||
195 | return x * self(x) |
||
196 | |||
197 | return integrate_spectrum(f, emin, emax, **kwargs) |
||
198 | |||
199 | def energy_flux_error(self, emin, emax, **kwargs): |
||
200 | """Compute energy flux in given energy range with error propagation. |
||
201 | |||
202 | .. math:: |
||
203 | |||
204 | G(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}}E \phi(E)dE |
||
205 | |||
206 | Parameters |
||
207 | ---------- |
||
208 | emin, emax : `~astropy.units.Quantity` |
||
209 | Lower bound of integration range. |
||
210 | **kwargs : dict |
||
211 | Keyword arguments passed to :func:`~gammapy.spectrum.integrate_spectrum` |
||
212 | |||
213 | Returns |
||
214 | ------- |
||
215 | energy_flux, energy_flux_error : tuple of `~astropy.units.Quantity` |
||
216 | Tuple of energy flux and energy flux error. |
||
217 | """ |
||
218 | emin = self._convert_energy(emin) |
||
219 | emax = self._convert_energy(emax) |
||
220 | |||
221 | unit = self.energy_flux(emin, emax, **kwargs).unit |
||
222 | upars = self.parameters._ufloats |
||
223 | |||
224 | def f(x): |
||
225 | return x * self.evaluate(x, **upars) |
||
226 | |||
227 | uarray = integrate_spectrum(f, emin.value, emax.value, **kwargs) |
||
228 | return self._parse_uarray(uarray) * unit |
||
229 | |||
230 | def to_dict(self): |
||
231 | """Convert to dict.""" |
||
232 | retval = self.parameters.to_dict() |
||
233 | retval["name"] = self.__class__.__name__ |
||
234 | return retval |
||
235 | |||
236 | @classmethod |
||
237 | def from_dict(cls, val): |
||
238 | """Create from dict.""" |
||
239 | classname = val.pop("name") |
||
240 | parameters = Parameters.from_dict(val) |
||
241 | model = globals()[classname]() |
||
242 | model.parameters = parameters |
||
243 | model.parameters.covariance = parameters.covariance |
||
244 | return model |
||
245 | |||
246 | def plot( |
||
247 | self, |
||
248 | energy_range, |
||
249 | ax=None, |
||
250 | energy_unit="TeV", |
||
251 | flux_unit="cm-2 s-1 TeV-1", |
||
252 | energy_power=0, |
||
253 | n_points=100, |
||
254 | **kwargs |
||
255 | ): |
||
256 | """Plot spectral model curve. |
||
257 | |||
258 | kwargs are forwarded to `matplotlib.pyplot.plot` |
||
259 | |||
260 | By default a log-log scaling of the axes is used, if you want to change |
||
261 | the y axis scaling to linear you can use: |
||
262 | |||
263 | .. code: |
||
264 | |||
265 | from gammapy.spectrum.models import PowerLaw |
||
266 | from astropy import units as u |
||
267 | |||
268 | pwl = ExponentialCutoffPowerLaw() |
||
269 | ax = pwl.plot(energy_range=(0.1, 100) * u.TeV) |
||
270 | ax.set_yscale('linear') |
||
271 | |||
272 | |||
273 | Parameters |
||
274 | ---------- |
||
275 | ax : `~matplotlib.axes.Axes`, optional |
||
276 | Axis |
||
277 | energy_range : `~astropy.units.Quantity` |
||
278 | Plot range |
||
279 | energy_unit : str, `~astropy.units.Unit`, optional |
||
280 | Unit of the energy axis |
||
281 | flux_unit : str, `~astropy.units.Unit`, optional |
||
282 | Unit of the flux axis |
||
283 | energy_power : int, optional |
||
284 | Power of energy to multiply flux axis with |
||
285 | n_points : int, optional |
||
286 | Number of evaluation nodes |
||
287 | |||
288 | Returns |
||
289 | ------- |
||
290 | ax : `~matplotlib.axes.Axes`, optional |
||
291 | Axis |
||
292 | """ |
||
293 | import matplotlib.pyplot as plt |
||
294 | |||
295 | ax = plt.gca() if ax is None else ax |
||
296 | |||
297 | emin, emax = energy_range |
||
298 | energy = EnergyBounds.equal_log_spacing(emin, emax, n_points, energy_unit) |
||
299 | |||
300 | # evaluate model |
||
301 | flux = self(energy).to(flux_unit) |
||
302 | |||
303 | y = self._plot_scale_flux(energy, flux, energy_power) |
||
304 | |||
305 | ax.plot(energy.value, y.value, **kwargs) |
||
306 | |||
307 | self._plot_format_ax(ax, energy, y, energy_power) |
||
308 | return ax |
||
309 | |||
310 | def plot_error( |
||
311 | self, |
||
312 | energy_range, |
||
313 | ax=None, |
||
314 | energy_unit="TeV", |
||
315 | flux_unit="cm-2 s-1 TeV-1", |
||
316 | energy_power=0, |
||
317 | n_points=100, |
||
318 | **kwargs |
||
319 | ): |
||
320 | """Plot spectral model error band. |
||
321 | |||
322 | .. note:: |
||
323 | |||
324 | This method calls ``ax.set_yscale("log", nonposy='clip')`` and |
||
325 | ``ax.set_xscale("log", nonposx='clip')`` to create a log-log representation. |
||
326 | The additional argument ``nonposx='clip'`` avoids artefacts in the plot, |
||
327 | when the error band extends to negative values (see also |
||
328 | https://github.com/matplotlib/matplotlib/issues/8623). |
||
329 | |||
330 | When you call ``plt.loglog()`` or ``plt.semilogy()`` explicitely in your |
||
331 | plotting code and the error band extends to negative values, it is not |
||
332 | shown correctly. To circumvent this issue also use |
||
333 | ``plt.loglog(nonposx='clip', nonposy='clip')`` |
||
334 | or ``plt.semilogy(nonposy='clip')``. |
||
335 | |||
336 | Parameters |
||
337 | ---------- |
||
338 | ax : `~matplotlib.axes.Axes`, optional |
||
339 | Axis |
||
340 | energy_range : `~astropy.units.Quantity` |
||
341 | Plot range |
||
342 | energy_unit : str, `~astropy.units.Unit`, optional |
||
343 | Unit of the energy axis |
||
344 | flux_unit : str, `~astropy.units.Unit`, optional |
||
345 | Unit of the flux axis |
||
346 | energy_power : int, optional |
||
347 | Power of energy to multiply flux axis with |
||
348 | n_points : int, optional |
||
349 | Number of evaluation nodes |
||
350 | **kwargs : dict |
||
351 | Keyword arguments forwarded to `matplotlib.pyplot.fill_between` |
||
352 | |||
353 | |||
354 | Returns |
||
355 | ------- |
||
356 | ax : `~matplotlib.axes.Axes`, optional |
||
357 | Axis |
||
358 | """ |
||
359 | import matplotlib.pyplot as plt |
||
360 | |||
361 | ax = plt.gca() if ax is None else ax |
||
362 | |||
363 | kwargs.setdefault("facecolor", "black") |
||
364 | kwargs.setdefault("alpha", 0.2) |
||
365 | kwargs.setdefault("linewidth", 0) |
||
366 | |||
367 | emin, emax = energy_range |
||
368 | energy = EnergyBounds.equal_log_spacing(emin, emax, n_points, energy_unit) |
||
369 | |||
370 | flux, flux_err = self.evaluate_error(energy).to(flux_unit) |
||
371 | |||
372 | y_lo = self._plot_scale_flux(energy, flux - flux_err, energy_power) |
||
373 | y_hi = self._plot_scale_flux(energy, flux + flux_err, energy_power) |
||
374 | |||
375 | where = (energy >= energy_range[0]) & (energy <= energy_range[1]) |
||
376 | ax.fill_between(energy.value, y_lo.value, y_hi.value, where=where, **kwargs) |
||
377 | |||
378 | self._plot_format_ax(ax, energy, y_lo, energy_power) |
||
379 | return ax |
||
380 | |||
381 | def _plot_format_ax(self, ax, energy, y, energy_power): |
||
382 | ax.set_xlabel("Energy [{}]".format(energy.unit)) |
||
383 | if energy_power > 0: |
||
384 | ax.set_ylabel("E{} * Flux [{}]".format(energy_power, y.unit)) |
||
385 | else: |
||
386 | ax.set_ylabel("Flux [{}]".format(y.unit)) |
||
387 | |||
388 | ax.set_xscale("log", nonposx="clip") |
||
389 | ax.set_yscale("log", nonposy="clip") |
||
390 | |||
391 | def _plot_scale_flux(self, energy, flux, energy_power): |
||
392 | try: |
||
393 | eunit = [_ for _ in flux.unit.bases if _.physical_type == "energy"][0] |
||
394 | except IndexError: |
||
395 | eunit = energy.unit |
||
396 | y = flux * np.power(energy, energy_power) |
||
397 | return y.to(flux.unit * eunit ** energy_power) |
||
398 | |||
399 | def spectral_index(self, energy, epsilon=1e-5): |
||
400 | """Compute spectral index at given energy. |
||
401 | |||
402 | Parameters |
||
403 | ---------- |
||
404 | energy : `~astropy.units.Quantity` |
||
405 | Energy at which to estimate the index |
||
406 | epsilon : float |
||
407 | Fractional energy increment to use for determining the spectral index. |
||
408 | |||
409 | Returns |
||
410 | ------- |
||
411 | index : float |
||
412 | Estimated spectral index. |
||
413 | """ |
||
414 | f1 = self(energy) |
||
415 | f2 = self(energy * (1 + epsilon)) |
||
416 | return np.log(f1 / f2) / np.log(1 + epsilon) |
||
417 | |||
418 | def inverse(self, value, emin=0.1 * u.TeV, emax=100 * u.TeV): |
||
419 | """Return energy for a given function value of the spectral model. |
||
420 | |||
421 | Calls the `scipy.optimize.brentq` numerical root finding method. |
||
422 | |||
423 | Parameters |
||
424 | ---------- |
||
425 | value : `~astropy.units.Quantity` |
||
426 | Function value of the spectral model. |
||
427 | emin : `~astropy.units.Quantity` |
||
428 | Lower bracket value in case solution is not unique. |
||
429 | emax : `~astropy.units.Quantity` |
||
430 | Upper bracket value in case solution is not unique. |
||
431 | |||
432 | Returns |
||
433 | ------- |
||
434 | energy : `~astropy.units.Quantity` |
||
435 | Energies at which the model has the given ``value``. |
||
436 | """ |
||
437 | from scipy.optimize import brentq |
||
438 | eunit = "TeV" |
||
439 | |||
440 | energies = [] |
||
441 | for val in np.atleast_1d(value): |
||
442 | |||
443 | def f(x): |
||
444 | # scale by 1e12 to achieve better precision |
||
445 | energy = u.Quantity(x, eunit, copy=False) |
||
446 | y = self(energy).to_value(value.unit) |
||
447 | return 1e12 * (y - val.value) |
||
448 | |||
449 | energy = brentq(f, emin.to_value(eunit), emax.to_value(eunit)) |
||
450 | energies.append(energy) |
||
451 | |||
452 | return u.Quantity(energies, eunit, copy=False) |
||
453 | |||
454 | def copy(self): |
||
455 | """A deep copy.""" |
||
456 | return copy.deepcopy(self) |
||
457 | |||
458 | |||
459 | class ConstantModel(SpectralModel): |
||
460 | r"""Constant model |
||
461 | |||
462 | .. math:: |
||
463 | |||
464 | \phi(E) = k |
||
465 | |||
466 | Parameters |
||
467 | ---------- |
||
468 | const : `~astropy.units.Quantity` |
||
469 | :math:`k` |
||
470 | """ |
||
471 | |||
472 | def __init__(self, const): |
||
473 | self.parameters = Parameters([Parameter("const", const)]) |
||
474 | |||
475 | @staticmethod |
||
476 | def evaluate(energy, const): |
||
477 | """Evaluate the model (static function).""" |
||
478 | return np.ones(np.atleast_1d(energy).shape) * const |
||
479 | |||
480 | |||
481 | class CompoundSpectralModel(SpectralModel): |
||
482 | """Represents the algebraic combination of two |
||
483 | `~gammapy.spectrum.models.SpectralModel` |
||
484 | |||
485 | """ |
||
486 | |||
487 | def __init__(self, model1, model2, operator): |
||
488 | self.model1 = model1 |
||
489 | self.model2 = model2 |
||
490 | self.operator = operator |
||
491 | |||
492 | # TODO: Think about how to deal with covariance matrix |
||
493 | @property |
||
494 | def parameters(self): |
||
495 | val = self.model1.parameters.parameters + self.model2.parameters.parameters |
||
496 | return Parameters(val) |
||
497 | |||
498 | @parameters.setter |
||
499 | def parameters(self, parameters): |
||
500 | idx = len(self.model1.parameters.parameters) |
||
501 | self.model1.parameters.parameters = parameters.parameters[:idx] |
||
502 | self.model2.parameters.parameters = parameters.parameters[idx:] |
||
503 | |||
504 | def __str__(self): |
||
505 | ss = self.__class__.__name__ |
||
506 | ss += "\n Component 1 : {}".format(self.model1) |
||
507 | ss += "\n Component 2 : {}".format(self.model2) |
||
508 | ss += "\n Operator : {}".format(self.operator) |
||
509 | return ss |
||
510 | |||
511 | def __call__(self, energy): |
||
512 | val1 = self.model1(energy) |
||
513 | val2 = self.model2(energy) |
||
514 | |||
515 | return self.operator(val1, val2) |
||
516 | |||
517 | def to_dict(self): |
||
518 | retval = dict() |
||
519 | retval["model1"] = self.model1.to_dict() |
||
520 | retval["model2"] = self.model2.to_dict() |
||
521 | retval["operator"] = self.operator |
||
522 | |||
523 | |||
524 | class PowerLaw(SpectralModel): |
||
525 | r"""Spectral power-law model. |
||
526 | |||
527 | .. math:: |
||
528 | |||
529 | \phi(E) = \phi_0 \cdot \left( \frac{E}{E_0} \right)^{-\Gamma} |
||
530 | |||
531 | Parameters |
||
532 | ---------- |
||
533 | index : `~astropy.units.Quantity` |
||
534 | :math:`\Gamma` |
||
535 | amplitude : `~astropy.units.Quantity` |
||
536 | :math:`\phi_0` |
||
537 | reference : `~astropy.units.Quantity` |
||
538 | :math:`E_0` |
||
539 | |||
540 | |||
541 | Examples |
||
542 | -------- |
||
543 | This is how to plot the default `PowerLaw` model: |
||
544 | |||
545 | .. code:: python |
||
546 | |||
547 | from astropy import units as u |
||
548 | from gammapy.spectrum.models import PowerLaw |
||
549 | |||
550 | pwl = PowerLaw() |
||
551 | pwl.plot(energy_range=[0.1, 100] * u.TeV) |
||
552 | plt.show() |
||
553 | |||
554 | """ |
||
555 | |||
556 | def __init__( |
||
557 | self, index=2.0, amplitude="1e-12 cm-2 s-1 TeV-1", reference="1 TeV" |
||
558 | ): |
||
559 | self.parameters = Parameters( |
||
560 | [ |
||
561 | Parameter("index", index), |
||
562 | Parameter("amplitude", amplitude), |
||
563 | Parameter("reference", reference, min=0, frozen=True), |
||
564 | ] |
||
565 | ) |
||
566 | |||
567 | @staticmethod |
||
568 | def evaluate(energy, index, amplitude, reference): |
||
569 | """Evaluate the model (static function).""" |
||
570 | return amplitude * np.power((energy / reference), -index) |
||
571 | |||
572 | def integral(self, emin, emax, **kwargs): |
||
573 | r"""Integrate power law analytically. |
||
574 | |||
575 | .. math:: |
||
576 | |||
577 | F(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}}\phi(E)dE = \left. |
||
578 | \phi_0 \frac{E_0}{-\Gamma + 1} \left( \frac{E}{E_0} \right)^{-\Gamma + 1} |
||
579 | \right \vert _{E_{min}}^{E_{max}} |
||
580 | |||
581 | Parameters |
||
582 | ---------- |
||
583 | emin, emax : `~astropy.units.Quantity` |
||
584 | Lower and upper bound of integration range |
||
585 | """ |
||
586 | # kwargs are passed to this function but not used |
||
587 | # this is to get a consistent API with SpectralModel.integral() |
||
588 | pars = self.parameters |
||
589 | |||
590 | if np.isclose(pars["index"].value, 1): |
||
591 | e_unit = emin.unit |
||
592 | prefactor = pars["amplitude"].quantity * pars["reference"].quantity.to( |
||
593 | e_unit |
||
594 | ) |
||
595 | upper = np.log(emax.to_value(e_unit)) |
||
596 | lower = np.log(emin.to_value(e_unit)) |
||
597 | else: |
||
598 | val = -1 * pars["index"].value + 1 |
||
599 | prefactor = pars["amplitude"].quantity * pars["reference"].quantity / val |
||
600 | upper = np.power((emax / pars["reference"].quantity), val) |
||
601 | lower = np.power((emin / pars["reference"].quantity), val) |
||
602 | |||
603 | return prefactor * (upper - lower) |
||
604 | |||
605 | def integral_error(self, emin, emax, **kwargs): |
||
606 | r"""Integrate power law analytically with error propagation. |
||
607 | |||
608 | Parameters |
||
609 | ---------- |
||
610 | emin, emax : `~astropy.units.Quantity` |
||
611 | Lower and upper bound of integration range. |
||
612 | |||
613 | Returns |
||
614 | ------- |
||
615 | integral, integral_error : tuple of `~astropy.units.Quantity` |
||
616 | Tuple of integral flux and integral flux error. |
||
617 | """ |
||
618 | # kwargs are passed to this function but not used |
||
619 | # this is to get a consistent API with SpectralModel.integral() |
||
620 | emin = self._convert_energy(emin) |
||
621 | emax = self._convert_energy(emax) |
||
622 | |||
623 | unit = self.integral(emin, emax, **kwargs).unit |
||
624 | upars = self.parameters._ufloats |
||
625 | |||
626 | if np.isclose(upars["index"].nominal_value, 1): |
||
627 | prefactor = upars["amplitude"] * upars["reference"] |
||
628 | upper = np.log(emax.value) |
||
629 | lower = np.log(emin.value) |
||
630 | else: |
||
631 | val = -1 * upars["index"] + 1 |
||
632 | prefactor = upars["amplitude"] * upars["reference"] / val |
||
633 | upper = np.power((emax.value / upars["reference"]), val) |
||
634 | lower = np.power((emin.value / upars["reference"]), val) |
||
635 | |||
636 | uarray = prefactor * (upper - lower) |
||
637 | return self._parse_uarray(uarray) * unit |
||
638 | |||
639 | def energy_flux(self, emin, emax): |
||
640 | r"""Compute energy flux in given energy range analytically. |
||
641 | |||
642 | .. math:: |
||
643 | |||
644 | G(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}}E \phi(E)dE = \left. |
||
645 | \phi_0 \frac{E_0^2}{-\Gamma + 2} \left( \frac{E}{E_0} \right)^{-\Gamma + 2} |
||
646 | \right \vert _{E_{min}}^{E_{max}} |
||
647 | |||
648 | Parameters |
||
649 | ---------- |
||
650 | emin, emax : `~astropy.units.Quantity` |
||
651 | Lower and upper bound of integration range. |
||
652 | """ |
||
653 | pars = self.parameters |
||
654 | val = -1 * pars["index"].value + 2 |
||
655 | |||
656 | if np.isclose(val, 0): |
||
657 | # see https://www.wolframalpha.com/input/?i=a+*+x+*+(x%2Fb)+%5E+(-2) |
||
658 | # for reference |
||
659 | temp = pars["amplitude"].quantity * pars["reference"].quantity ** 2 |
||
660 | return temp * np.log(emax / emin) |
||
661 | else: |
||
662 | prefactor = ( |
||
663 | pars["amplitude"].quantity * pars["reference"].quantity ** 2 / val |
||
664 | ) |
||
665 | upper = (emax / pars["reference"].quantity) ** val |
||
666 | lower = (emin / pars["reference"].quantity) ** val |
||
667 | return prefactor * (upper - lower) |
||
668 | |||
669 | def energy_flux_error(self, emin, emax, **kwargs): |
||
670 | r"""Compute energy flux in given energy range analytically with error propagation. |
||
671 | |||
672 | Parameters |
||
673 | ---------- |
||
674 | emin, emax : `~astropy.units.Quantity` |
||
675 | Lower and upper bound of integration range. |
||
676 | |||
677 | Returns |
||
678 | ------- |
||
679 | energy_flux, energy_flux_error : tuple of `~astropy.units.Quantity` |
||
680 | Tuple of energy flux and energy flux error. |
||
681 | """ |
||
682 | emin = self._convert_energy(emin) |
||
683 | emax = self._convert_energy(emax) |
||
684 | |||
685 | unit = self.energy_flux(emin, emax, **kwargs).unit |
||
686 | upars = self.parameters._ufloats |
||
687 | |||
688 | val = -1 * upars["index"] + 2 |
||
689 | |||
690 | if np.isclose(val.nominal_value, 0): |
||
691 | # see https://www.wolframalpha.com/input/?i=a+*+x+*+(x%2Fb)+%5E+(-2) |
||
692 | # for reference |
||
693 | temp = upars["amplitude"] * upars["reference"] ** 2 |
||
694 | uarray = temp * np.log(emax.value / emin.value) |
||
695 | else: |
||
696 | prefactor = upars["amplitude"] * upars["reference"] ** 2 / val |
||
697 | upper = (emax.value / upars["reference"]) ** val |
||
698 | lower = (emin.value / upars["reference"]) ** val |
||
699 | uarray = prefactor * (upper - lower) |
||
700 | |||
701 | return self._parse_uarray(uarray) * unit |
||
702 | |||
703 | def inverse(self, value): |
||
704 | """Return energy for a given function value of the spectral model. |
||
705 | |||
706 | Parameters |
||
707 | ---------- |
||
708 | value : `~astropy.units.Quantity` |
||
709 | Function value of the spectral model. |
||
710 | """ |
||
711 | p = self.parameters |
||
712 | base = value / p["amplitude"].quantity |
||
713 | return p["reference"].quantity * np.power(base, -1.0 / p["index"].value) |
||
714 | |||
715 | |||
716 | class PowerLaw2(SpectralModel): |
||
717 | r"""Spectral power-law model with integral as amplitude parameter. |
||
718 | |||
719 | See also: https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/source_models.html |
||
720 | |||
721 | .. math:: |
||
722 | |||
723 | \phi(E) = F_0 \cdot \frac{\Gamma + 1}{E_{0, max}^{-\Gamma + 1} |
||
724 | - E_{0, min}^{-\Gamma + 1}} \cdot E^{-\Gamma} |
||
725 | |||
726 | Parameters |
||
727 | ---------- |
||
728 | index : `~astropy.units.Quantity` |
||
729 | Spectral index :math:`\Gamma` |
||
730 | amplitude : `~astropy.units.Quantity` |
||
731 | Integral flux :math:`F_0`. |
||
732 | emin : `~astropy.units.Quantity` |
||
733 | Lower energy limit :math:`E_{0, min}`. |
||
734 | emax : `~astropy.units.Quantity` |
||
735 | Upper energy limit :math:`E_{0, max}`. |
||
736 | |||
737 | Examples |
||
738 | -------- |
||
739 | This is how to plot the default `PowerLaw2` model: |
||
740 | |||
741 | .. code:: python |
||
742 | |||
743 | from astropy import units as u |
||
744 | from gammapy.spectrum.models import PowerLaw2 |
||
745 | |||
746 | pwl2 = PowerLaw2() |
||
747 | pwl2.plot(energy_range=[0.1, 100] * u.TeV) |
||
748 | plt.show() |
||
749 | """ |
||
750 | |||
751 | def __init__( |
||
752 | self, |
||
753 | amplitude="1e-12 cm-2 s-1", |
||
754 | index=2, |
||
755 | emin="0.1 TeV", |
||
756 | emax="100 TeV", |
||
757 | ): |
||
758 | self.parameters = Parameters( |
||
759 | [ |
||
760 | Parameter("amplitude", amplitude), |
||
761 | Parameter("index", index), |
||
762 | Parameter("emin", emin, frozen=True), |
||
763 | Parameter("emax", emax, frozen=True), |
||
764 | ] |
||
765 | ) |
||
766 | |||
767 | @staticmethod |
||
768 | def evaluate(energy, amplitude, index, emin, emax): |
||
769 | """Evaluate the model (static function).""" |
||
770 | top = -index + 1 |
||
771 | |||
772 | # to get the energies dimensionless we use a modified formula |
||
773 | bottom = emax - emin * (emin / emax) ** (-index) |
||
774 | return amplitude * (top / bottom) * np.power(energy / emax, -index) |
||
775 | |||
776 | def integral(self, emin, emax, **kwargs): |
||
777 | r"""Integrate power law analytically. |
||
778 | |||
779 | .. math:: |
||
780 | |||
781 | F(E_{min}, E_{max}) = F_0 \cdot \frac{E_{max}^{\Gamma + 1} \ |
||
782 | - E_{min}^{\Gamma + 1}}{E_{0, max}^{\Gamma + 1} \ |
||
783 | - E_{0, min}^{\Gamma + 1}} |
||
784 | |||
785 | Parameters |
||
786 | ---------- |
||
787 | emin, emax : `~astropy.units.Quantity` |
||
788 | Lower and upper bound of integration range. |
||
789 | """ |
||
790 | pars = self.parameters |
||
791 | |||
792 | temp1 = np.power(emax, -pars["index"].value + 1) |
||
793 | temp2 = np.power(emin, -pars["index"].value + 1) |
||
794 | top = temp1 - temp2 |
||
795 | |||
796 | temp1 = np.power(pars["emax"].quantity, -pars["index"].value + 1) |
||
797 | temp2 = np.power(pars["emin"].quantity, -pars["index"].value + 1) |
||
798 | bottom = temp1 - temp2 |
||
799 | |||
800 | return pars["amplitude"].quantity * top / bottom |
||
801 | |||
802 | def integral_error(self, emin, emax, **kwargs): |
||
803 | r"""Integrate power law analytically with error propagation. |
||
804 | |||
805 | Parameters |
||
806 | ---------- |
||
807 | emin, emax : `~astropy.units.Quantity` |
||
808 | Lower and upper bound of integration range. |
||
809 | |||
810 | Returns |
||
811 | ------- |
||
812 | integral, integral_error : tuple of `~astropy.units.Quantity` |
||
813 | Tuple of integral flux and integral flux error. |
||
814 | """ |
||
815 | emin = self._convert_energy(emin) |
||
816 | emax = self._convert_energy(emax) |
||
817 | |||
818 | unit = self.integral(emin, emax, **kwargs).unit |
||
819 | upars = self.parameters._ufloats |
||
820 | |||
821 | temp1 = np.power(emax.value, -upars["index"] + 1) |
||
822 | temp2 = np.power(emin.value, -upars["index"] + 1) |
||
823 | top = temp1 - temp2 |
||
824 | |||
825 | temp1 = np.power(upars["emax"], -upars["index"] + 1) |
||
826 | temp2 = np.power(upars["emin"], -upars["index"] + 1) |
||
827 | bottom = temp1 - temp2 |
||
828 | |||
829 | uarray = upars["amplitude"] * top / bottom |
||
830 | return self._parse_uarray(uarray) * unit |
||
831 | |||
832 | def inverse(self, value): |
||
833 | """Return energy for a given function value of the spectral model. |
||
834 | |||
835 | Parameters |
||
836 | ---------- |
||
837 | value : `~astropy.units.Quantity` |
||
838 | Function value of the spectral model. |
||
839 | """ |
||
840 | p = self.parameters |
||
841 | amplitude, index, emin, emax = (p["amplitude"].quantity, p["index"].value, |
||
842 | p["emin"].quantity, p["emax"].quantity) |
||
843 | |||
844 | # to get the energies dimensionless we use a modified formula |
||
845 | top = -index + 1 |
||
846 | bottom = emax - emin * (emin / emax) ** (-index) |
||
847 | term = (bottom / top) * (value / amplitude) |
||
848 | return np.power(term.to_value(""), -1.0 / index) * emax |
||
849 | |||
850 | |||
851 | class ExponentialCutoffPowerLaw(SpectralModel): |
||
852 | r"""Spectral exponential cutoff power-law model. |
||
853 | |||
854 | .. math:: |
||
855 | |||
856 | \phi(E) = \phi_0 \cdot \left(\frac{E}{E_0}\right)^{-\Gamma} \exp(-\lambda E) |
||
857 | |||
858 | Parameters |
||
859 | ---------- |
||
860 | index : `~astropy.units.Quantity` |
||
861 | :math:`\Gamma` |
||
862 | amplitude : `~astropy.units.Quantity` |
||
863 | :math:`\phi_0` |
||
864 | reference : `~astropy.units.Quantity` |
||
865 | :math:`E_0` |
||
866 | lambda_ : `~astropy.units.Quantity` |
||
867 | :math:`\lambda` |
||
868 | |||
869 | Examples |
||
870 | -------- |
||
871 | This is how to plot the default `ExponentialCutoffPowerLaw` model: |
||
872 | |||
873 | .. code:: python |
||
874 | |||
875 | from astropy import units as u |
||
876 | from gammapy.spectrum.models import ExponentialCutoffPowerLaw |
||
877 | |||
878 | ecpl = ExponentialCutoffPowerLaw() |
||
879 | ecpl.plot(energy_range=[0.1, 100] * u.TeV) |
||
880 | plt.show() |
||
881 | """ |
||
882 | |||
883 | def __init__( |
||
884 | self, |
||
885 | index=1.5, |
||
886 | amplitude="1e-12 cm-2 s-1 TeV-1", |
||
887 | reference="1 TeV", |
||
888 | lambda_="0.1 TeV-1", |
||
889 | ): |
||
890 | self.parameters = Parameters( |
||
891 | [ |
||
892 | Parameter("index", index), |
||
893 | Parameter("amplitude", amplitude), |
||
894 | Parameter("reference", reference, frozen=True), |
||
895 | Parameter("lambda_", lambda_), |
||
896 | ] |
||
897 | ) |
||
898 | |||
899 | @staticmethod |
||
900 | def evaluate(energy, index, amplitude, reference, lambda_): |
||
901 | """Evaluate the model (static function).""" |
||
902 | pwl = amplitude * (energy / reference) ** (-index) |
||
903 | try: |
||
904 | cutoff = np.exp(-energy * lambda_) |
||
905 | except AttributeError: |
||
906 | from uncertainties.unumpy import exp |
||
907 | |||
908 | cutoff = exp(-energy * lambda_) |
||
909 | return pwl * cutoff |
||
910 | |||
911 | @property |
||
912 | def e_peak(self): |
||
913 | r"""Spectral energy distribution peak energy (`~astropy.utils.Quantity`). |
||
914 | |||
915 | This is the peak in E^2 x dN/dE and is given by: |
||
916 | |||
917 | .. math:: |
||
918 | |||
919 | E_{Peak} = (2 - \Gamma) / \lambda |
||
920 | |||
921 | """ |
||
922 | p = self.parameters |
||
923 | reference = p["reference"].quantity |
||
924 | index = p["index"].quantity |
||
925 | lambda_ = p["lambda_"].quantity |
||
926 | if index >= 2: |
||
927 | return np.nan * reference.unit |
||
928 | else: |
||
929 | return (2 - index) / lambda_ |
||
930 | |||
931 | |||
932 | class ExponentialCutoffPowerLaw3FGL(SpectralModel): |
||
933 | r"""Spectral exponential cutoff power-law model used for 3FGL. |
||
934 | |||
935 | Note that the parametrization is different from `ExponentialCutoffPowerLaw`: |
||
936 | |||
937 | .. math:: |
||
938 | |||
939 | \phi(E) = \phi_0 \cdot \left(\frac{E}{E_0}\right)^{-\Gamma} |
||
940 | \exp \left( \frac{E_0 - E}{E_{C}} \right) |
||
941 | |||
942 | Parameters |
||
943 | ---------- |
||
944 | index : `~astropy.units.Quantity` |
||
945 | :math:`\Gamma` |
||
946 | amplitude : `~astropy.units.Quantity` |
||
947 | :math:`\phi_0` |
||
948 | reference : `~astropy.units.Quantity` |
||
949 | :math:`E_0` |
||
950 | ecut : `~astropy.units.Quantity` |
||
951 | :math:`E_{C}` |
||
952 | |||
953 | Examples |
||
954 | -------- |
||
955 | This is how to plot the default `ExponentialCutoffPowerLaw3FGL` model: |
||
956 | |||
957 | .. code:: python |
||
958 | |||
959 | from astropy import units as u |
||
960 | from gammapy.spectrum.models import ExponentialCutoffPowerLaw3FGL |
||
961 | |||
962 | ecpl_3fgl = ExponentialCutoffPowerLaw3FGL() |
||
963 | ecpl_3fgl.plot(energy_range=[0.1, 100] * u.TeV) |
||
964 | plt.show() |
||
965 | """ |
||
966 | |||
967 | def __init__( |
||
968 | self, |
||
969 | index=1.5, |
||
970 | amplitude="1e-12 cm-2 s-1 TeV-1", |
||
971 | reference="1 TeV", |
||
972 | ecut="10 TeV", |
||
973 | ): |
||
974 | self.parameters = Parameters( |
||
975 | [ |
||
976 | Parameter("index", index), |
||
977 | Parameter("amplitude", amplitude), |
||
978 | Parameter("reference", reference, frozen=True), |
||
979 | Parameter("ecut", ecut), |
||
980 | ] |
||
981 | ) |
||
982 | |||
983 | @staticmethod |
||
984 | def evaluate(energy, index, amplitude, reference, ecut): |
||
985 | """Evaluate the model (static function).""" |
||
986 | pwl = amplitude * (energy / reference) ** (-index) |
||
987 | try: |
||
988 | cutoff = np.exp((reference - energy) / ecut) |
||
989 | except AttributeError: |
||
990 | from uncertainties.unumpy import exp |
||
991 | |||
992 | cutoff = exp((reference - energy) / ecut) |
||
993 | return pwl * cutoff |
||
994 | |||
995 | |||
996 | class PLSuperExpCutoff3FGL(SpectralModel): |
||
997 | r"""Spectral super exponential cutoff power-law model used for 3FGL. |
||
998 | |||
999 | .. math:: |
||
1000 | |||
1001 | \phi(E) = \phi_0 \cdot \left(\frac{E}{E_0}\right)^{-\Gamma_1} |
||
1002 | \exp \left( \left(\frac{E_0}{E_{C}} \right)^{\Gamma_2} - |
||
1003 | \left(\frac{E}{E_{C}} \right)^{\Gamma_2} |
||
1004 | \right) |
||
1005 | |||
1006 | Parameters |
||
1007 | ---------- |
||
1008 | index_1 : `~astropy.units.Quantity` |
||
1009 | :math:`\Gamma_1` |
||
1010 | index_2 : `~astropy.units.Quantity` |
||
1011 | :math:`\Gamma_2` |
||
1012 | amplitude : `~astropy.units.Quantity` |
||
1013 | :math:`\phi_0` |
||
1014 | reference : `~astropy.units.Quantity` |
||
1015 | :math:`E_0` |
||
1016 | ecut : `~astropy.units.Quantity` |
||
1017 | :math:`E_{C}` |
||
1018 | |||
1019 | Examples |
||
1020 | -------- |
||
1021 | This is how to plot the default `PLSuperExpCutoff3FGL` model: |
||
1022 | |||
1023 | .. code:: python |
||
1024 | |||
1025 | from astropy import units as u |
||
1026 | from gammapy.spectrum.models import PLSuperExpCutoff3FGL |
||
1027 | |||
1028 | secpl_3fgl = PLSuperExpCutoff3FGL() |
||
1029 | secpl_3fgl.plot(energy_range=[0.1, 100] * u.TeV) |
||
1030 | plt.show() |
||
1031 | """ |
||
1032 | |||
1033 | def __init__( |
||
1034 | self, |
||
1035 | index_1=1.5, |
||
1036 | index_2=2, |
||
1037 | amplitude="1e-12 cm-2 s-1 TeV-1", |
||
1038 | reference="1 TeV", |
||
1039 | ecut="10 TeV", |
||
1040 | ): |
||
1041 | self.parameters = Parameters( |
||
1042 | [ |
||
1043 | Parameter("index_1", index_1), |
||
1044 | Parameter("index_2", index_2), |
||
1045 | Parameter("amplitude", amplitude), |
||
1046 | Parameter("reference", reference, frozen=True), |
||
1047 | Parameter("ecut", ecut), |
||
1048 | ] |
||
1049 | ) |
||
1050 | |||
1051 | @staticmethod |
||
1052 | def evaluate(energy, amplitude, reference, ecut, index_1, index_2): |
||
1053 | """Evaluate the model (static function).""" |
||
1054 | pwl = amplitude * (energy / reference) ** (-index_1) |
||
1055 | try: |
||
1056 | cutoff = np.exp((reference / ecut) ** index_2 - (energy / ecut) ** index_2) |
||
1057 | except AttributeError: |
||
1058 | from uncertainties.unumpy import exp |
||
1059 | |||
1060 | cutoff = exp((reference / ecut) ** index_2 - (energy / ecut) ** index_2) |
||
1061 | return pwl * cutoff |
||
1062 | |||
1063 | |||
1064 | class LogParabola(SpectralModel): |
||
1065 | r"""Spectral log parabola model. |
||
1066 | |||
1067 | .. math:: |
||
1068 | |||
1069 | \phi(E) = \phi_0 \left( \frac{E}{E_0} \right) ^ { |
||
1070 | - \alpha - \beta \log{ \left( \frac{E}{E_0} \right) } |
||
1071 | } |
||
1072 | |||
1073 | Note that :math:`log` refers to the natural logarithm. This is consistent |
||
1074 | with the `Fermi Science Tools |
||
1075 | <https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/source_models.html>`_ |
||
1076 | and `ctools |
||
1077 | <http://cta.irap.omp.eu/ctools-devel/users/user_manual/getting_started/models.html#log-parabola>`_. |
||
1078 | The `Sherpa <http://cxc.harvard.edu/sherpa/ahelp/logparabola.html_ |
||
1079 | package>`_ package, however, uses :math:`log_{10}`. If you have |
||
1080 | parametrization based on :math:`log_{10}` you can use the |
||
1081 | :func:`~gammapy.spectrum.models.LogParabola.from_log10` method. |
||
1082 | |||
1083 | Parameters |
||
1084 | ---------- |
||
1085 | amplitude : `~astropy.units.Quantity` |
||
1086 | :math:`\phi_0` |
||
1087 | reference : `~astropy.units.Quantity` |
||
1088 | :math:`E_0` |
||
1089 | alpha : `~astropy.units.Quantity` |
||
1090 | :math:`\alpha` |
||
1091 | beta : `~astropy.units.Quantity` |
||
1092 | :math:`\beta` |
||
1093 | |||
1094 | Examples |
||
1095 | -------- |
||
1096 | This is how to plot the default `LogParabola` model: |
||
1097 | |||
1098 | .. code:: python |
||
1099 | |||
1100 | from astropy import units as u |
||
1101 | from gammapy.spectrum.models import LogParabola |
||
1102 | |||
1103 | log_parabola = LogParabola() |
||
1104 | log_parabola.plot(energy_range=[0.1, 100] * u.TeV) |
||
1105 | plt.show() |
||
1106 | """ |
||
1107 | |||
1108 | def __init__( |
||
1109 | self, |
||
1110 | amplitude="1e-12 cm-2 s-1 TeV-1", |
||
1111 | reference="10 TeV", |
||
1112 | alpha=2, |
||
1113 | beta=1, |
||
1114 | ): |
||
1115 | self.parameters = Parameters( |
||
1116 | [ |
||
1117 | Parameter("amplitude", amplitude), |
||
1118 | Parameter("reference", reference, frozen=True), |
||
1119 | Parameter("alpha", alpha), |
||
1120 | Parameter("beta", beta), |
||
1121 | ] |
||
1122 | ) |
||
1123 | |||
1124 | @classmethod |
||
1125 | def from_log10(cls, amplitude, reference, alpha, beta): |
||
1126 | """Construct LogParabola from :math:`log_{10}` parametrization""" |
||
1127 | beta_ = beta / np.log(10) |
||
1128 | return cls(amplitude=amplitude, reference=reference, alpha=alpha, beta=beta_) |
||
1129 | |||
1130 | @staticmethod |
||
1131 | def evaluate(energy, amplitude, reference, alpha, beta): |
||
1132 | """Evaluate the model (static function).""" |
||
1133 | try: |
||
1134 | xx = (energy / reference).to("") |
||
1135 | exponent = -alpha - beta * np.log(xx) |
||
1136 | except AttributeError: |
||
1137 | from uncertainties.unumpy import log |
||
1138 | |||
1139 | xx = energy / reference |
||
1140 | exponent = -alpha - beta * log(xx) |
||
1141 | return amplitude * np.power(xx, exponent) |
||
1142 | |||
1143 | @property |
||
1144 | def e_peak(self): |
||
1145 | r"""Spectral energy distribution peak energy (`~astropy.utils.Quantity`). |
||
1146 | |||
1147 | This is the peak in E^2 x dN/dE and is given by: |
||
1148 | |||
1149 | .. math:: |
||
1150 | |||
1151 | E_{Peak} = E_{0} \exp{ (2 - \alpha) / (2 * \beta)} |
||
1152 | |||
1153 | """ |
||
1154 | p = self.parameters |
||
1155 | reference = p["reference"].quantity |
||
1156 | alpha = p["alpha"].quantity |
||
1157 | beta = p["beta"].quantity |
||
1158 | return reference * np.exp((2 - alpha) / (2 * beta)) |
||
1159 | |||
1160 | |||
1161 | class TableModel(SpectralModel): |
||
1162 | """A model generated from a table of energy and value arrays. |
||
1163 | |||
1164 | the units returned will be the units of the values array provided at |
||
1165 | initialization. The model will return values interpolated in |
||
1166 | log-space, returning 0 for energies outside of the limits of the provided |
||
1167 | energy array. |
||
1168 | |||
1169 | Class implementation follows closely what has been done in |
||
1170 | `naima.models.TableModel` |
||
1171 | |||
1172 | Parameters |
||
1173 | ---------- |
||
1174 | energy : `~astropy.units.Quantity` array |
||
1175 | Array of energies at which the model values are given |
||
1176 | values : array |
||
1177 | Array with the values of the model at energies ``energy``. |
||
1178 | norm : float |
||
1179 | Model scale that is multiplied to the supplied arrays. Defaults to 1. |
||
1180 | values_scale : {'log', 'lin', 'sqrt'} |
||
1181 | Interpolation scaling applied to values. If the values vary over many magnitudes |
||
1182 | a 'log' scaling is recommended. |
||
1183 | interp_kwargs : dict |
||
1184 | Interpolation keyword arguments pass to `scipy.interpolate.interp1d`. |
||
1185 | By default all values outside the interpolation range are set to zero. |
||
1186 | If you want to apply linear extrapolation you can pass `interp_kwargs={'fill_value': |
||
1187 | 'extrapolate', 'kind': 'linear'}` |
||
1188 | meta : dict, optional |
||
1189 | Meta information, meta['filename'] will be used for serialization |
||
1190 | """ |
||
1191 | |||
1192 | def __init__( |
||
1193 | self, energy, values, norm=1, values_scale="log", interp_kwargs=None, meta=None |
||
1194 | ): |
||
1195 | self.parameters = Parameters([Parameter("norm", norm, min=0, unit="")]) |
||
1196 | self.energy = energy |
||
1197 | self.values = values |
||
1198 | self.meta = dict() if meta is None else meta |
||
1199 | |||
1200 | interp_kwargs = interp_kwargs or {} |
||
1201 | interp_kwargs.setdefault("values_scale", "log") |
||
1202 | |||
1203 | self._evaluate = ScaledRegularGridInterpolator( |
||
1204 | points=(np.log(energy.value),), values=values.value, **interp_kwargs |
||
1205 | ) |
||
1206 | |||
1207 | @classmethod |
||
1208 | def read_xspec_model(cls, filename, param, **kwargs): |
||
1209 | """Read XSPEC table model |
||
1210 | |||
1211 | The input is a table containing absorbed values from a XSPEC model as a |
||
1212 | function of energy. |
||
1213 | |||
1214 | TODO: Format of the file should be described and discussed in |
||
1215 | https://gamma-astro-data-formats.readthedocs.io/en/latest/index.html |
||
1216 | |||
1217 | Parameters |
||
1218 | ---------- |
||
1219 | filename : str |
||
1220 | File containing the XSPEC model |
||
1221 | param : float |
||
1222 | Model parameter value |
||
1223 | |||
1224 | Examples |
||
1225 | -------- |
||
1226 | Fill table from an EBL model (Franceschini, 2008) |
||
1227 | |||
1228 | >>> from gammapy.spectrum.models import TableModel |
||
1229 | >>> filename = '$GAMMAPY_DATA/ebl/ebl_franceschini.fits.gz' |
||
1230 | >>> table_model = TableModel.read_xspec_model(filename=filename, param=0.3) |
||
1231 | """ |
||
1232 | filename = str(make_path(filename)) |
||
1233 | |||
1234 | # Check if parameter value is in range |
||
1235 | table_param = Table.read(filename, hdu="PARAMETERS") |
||
1236 | param_min = table_param["MINIMUM"] |
||
1237 | param_max = table_param["MAXIMUM"] |
||
1238 | if param < param_min or param > param_max: |
||
1239 | err = "Parameter out of range, param={}, param_min={}, param_max={}".format( |
||
1240 | param, param_min, param_max |
||
1241 | ) |
||
1242 | raise ValueError(err) |
||
1243 | |||
1244 | # Get energy values |
||
1245 | table_energy = Table.read(filename, hdu="ENERGIES") |
||
1246 | energy_lo = table_energy["ENERG_LO"] |
||
1247 | energy_hi = table_energy["ENERG_HI"] |
||
1248 | |||
1249 | # Hack while format is not fixed, energy values are in keV |
||
1250 | energy_bounds = EnergyBounds.from_lower_and_upper_bounds( |
||
1251 | lower=energy_lo, upper=energy_hi, unit=u.keV |
||
1252 | ) |
||
1253 | energy = energy_bounds.log_centers |
||
1254 | |||
1255 | # Get spectrum values (no interpolation, take closest value for param) |
||
1256 | table_spectra = Table.read(filename, hdu="SPECTRA") |
||
1257 | idx = np.abs(table_spectra["PARAMVAL"] - param).argmin() |
||
1258 | values = u.Quantity(table_spectra[idx][1], "", copy=False) # no dimension |
||
1259 | |||
1260 | kwargs.setdefault("values_scale", "lin") |
||
1261 | return cls(energy=energy, values=values, **kwargs) |
||
1262 | |||
1263 | @classmethod |
||
1264 | def read_fermi_isotropic_model(cls, filename, **kwargs): |
||
1265 | """Read Fermi isotropic diffuse model |
||
1266 | |||
1267 | see `LAT Background models <https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html>`_ |
||
1268 | |||
1269 | Parameters |
||
1270 | ---------- |
||
1271 | filename : str |
||
1272 | filename |
||
1273 | """ |
||
1274 | filename = str(make_path(filename)) |
||
1275 | vals = np.loadtxt(filename) |
||
1276 | energy = u.Quantity(vals[:, 0], "MeV", copy=False) |
||
1277 | values = u.Quantity(vals[:, 1], "MeV-1 s-1 cm-2 sr-1", copy=False) |
||
1278 | return cls(energy=energy, values=values, **kwargs) |
||
1279 | |||
1280 | def evaluate(self, energy, norm): |
||
1281 | """Evaluate the model (static function).""" |
||
1282 | x = np.log(energy.to_value(self.energy.unit)) |
||
1283 | vals = self._evaluate(x, clip=True) |
||
1284 | vals = np.reshape(vals, x.shape) |
||
1285 | return u.Quantity(norm.value * vals, self.values.unit, copy=False) |
||
1286 | |||
1287 | |||
1288 | class Absorption(object): |
||
1289 | """Gamma-ray absorption models. |
||
1290 | |||
1291 | Parameters |
||
1292 | ---------- |
||
1293 | energy_lo, energy_hi : `~astropy.units.Quantity` |
||
1294 | Lower and upper bin edges of energy axis |
||
1295 | param_lo, param_hi : `~astropy.units.Quantity` |
||
1296 | Lower and upper bin edges of parameter axis |
||
1297 | data : `~astropy.units.Quantity` |
||
1298 | Model value |
||
1299 | |||
1300 | Examples |
||
1301 | -------- |
||
1302 | Create and plot EBL absorption models for a redshift of 0.5: |
||
1303 | |||
1304 | .. plot:: |
||
1305 | :include-source: |
||
1306 | |||
1307 | import matplotlib.pyplot as plt |
||
1308 | import astropy.units as u |
||
1309 | from gammapy.spectrum.models import Absorption |
||
1310 | |||
1311 | # Load tables for z=0.5 |
||
1312 | redshift = 0.5 |
||
1313 | dominguez = Absorption.read_builtin('dominguez').table_model(redshift) |
||
1314 | franceschini = Absorption.read_builtin('franceschini').table_model(redshift) |
||
1315 | finke = Absorption.read_builtin('finke').table_model(redshift) |
||
1316 | |||
1317 | # start customised plot |
||
1318 | energy_range = [0.08, 3] * u.TeV |
||
1319 | ax = plt.gca() |
||
1320 | opts = dict(energy_range=energy_range, energy_unit='TeV', ax=ax, flux_unit='') |
||
1321 | franceschini.plot(label='Franceschini 2008', **opts) |
||
1322 | finke.plot(label='Finke 2010', **opts) |
||
1323 | dominguez.plot(label='Dominguez 2011', **opts) |
||
1324 | |||
1325 | # tune plot |
||
1326 | ax.set_ylabel(r'Absorption coefficient [$\exp{(-\tau(E))}$]') |
||
1327 | ax.set_xlim(energy_range.value) # we get ride of units |
||
1328 | ax.set_ylim([1.e-4, 2.]) |
||
1329 | ax.set_yscale('log') |
||
1330 | ax.set_title('EBL models (z=' + str(redshift) + ')') |
||
1331 | plt.grid(which='both') |
||
1332 | plt.legend(loc='best') # legend |
||
1333 | |||
1334 | # show plot |
||
1335 | plt.show() |
||
1336 | """ |
||
1337 | |||
1338 | def __init__(self, energy_lo, energy_hi, param_lo, param_hi, data): |
||
1339 | axes = [ |
||
1340 | BinnedDataAxis( |
||
1341 | param_lo, param_hi, interpolation_mode="linear", name="parameter" |
||
1342 | ), |
||
1343 | BinnedDataAxis( |
||
1344 | energy_lo, energy_hi, interpolation_mode="log", name="energy" |
||
1345 | ), |
||
1346 | ] |
||
1347 | |||
1348 | self.data = NDDataArray(axes=axes, data=data) |
||
1349 | self.data.default_interp_kwargs["fill_value"] = None |
||
1350 | |||
1351 | @classmethod |
||
1352 | def read(cls, filename): |
||
1353 | """Build object from an XSPEC model. |
||
1354 | |||
1355 | Todo: Format of XSPEC binary files should be referenced at https://gamma-astro-data-formats.readthedocs.io/en/latest/ |
||
1356 | |||
1357 | Parameters |
||
1358 | ---------- |
||
1359 | filename : str |
||
1360 | File containing the model. |
||
1361 | """ |
||
1362 | # Create EBL data array |
||
1363 | filename = str(make_path(filename)) |
||
1364 | table_param = Table.read(filename, hdu="PARAMETERS") |
||
1365 | |||
1366 | par_min = table_param["MINIMUM"] |
||
1367 | par_max = table_param["MAXIMUM"] |
||
1368 | |||
1369 | par_array = table_param[0]["VALUE"] |
||
1370 | par_delta = np.diff(par_array) * 0.5 |
||
1371 | |||
1372 | param_lo, param_hi = par_array, par_array # initialisation |
||
1373 | param_lo[0] = par_min - par_delta[0] |
||
1374 | param_lo[1:] -= par_delta |
||
1375 | param_hi[:-1] += par_delta |
||
1376 | param_hi[-1] = par_max |
||
1377 | |||
1378 | # Get energy values |
||
1379 | table_energy = Table.read(filename, hdu="ENERGIES") |
||
1380 | energy_lo = u.Quantity(table_energy["ENERG_LO"], "keV", copy=False) # unit not stored in file |
||
1381 | energy_hi = u.Quantity(table_energy["ENERG_HI"], "keV", copy=False) # unit not stored in file |
||
1382 | |||
1383 | # Get spectrum values |
||
1384 | table_spectra = Table.read(filename, hdu="SPECTRA") |
||
1385 | data = table_spectra["INTPSPEC"].data |
||
1386 | |||
1387 | return cls( |
||
1388 | energy_lo=energy_lo, |
||
1389 | energy_hi=energy_hi, |
||
1390 | param_lo=param_lo, |
||
1391 | param_hi=param_hi, |
||
1392 | data=data, |
||
1393 | ) |
||
1394 | |||
1395 | @classmethod |
||
1396 | def read_builtin(cls, name): |
||
1397 | """Read one of the built-in absorption models. |
||
1398 | |||
1399 | Parameters |
||
1400 | ---------- |
||
1401 | name : {'franceschini', 'dominguez', 'finke'} |
||
1402 | name of one of the available model in gammapy-extra |
||
1403 | |||
1404 | References |
||
1405 | ---------- |
||
1406 | .. [1] Franceschini et al., "Extragalactic optical-infrared background radiation, its time evolution and the cosmic photon-photon opacity", |
||
1407 | `Link <http://adsabs.harvard.edu/abs/2008A%26A...487..837F>`__ |
||
1408 | .. [2] Dominguez et al., " Extragalactic background light inferred from AEGIS galaxy-SED-type fractions" |
||
1409 | `Link <http://adsabs.harvard.edu/abs/2011MNRAS.410.2556D>`__ |
||
1410 | .. [3] Finke et al., "Modeling the Extragalactic Background Light from Stars and Dust" |
||
1411 | `Link <http://adsabs.harvard.edu/abs/2010ApJ...712..238F>`__ |
||
1412 | """ |
||
1413 | models = dict() |
||
1414 | models["franceschini"] = "$GAMMAPY_DATA/ebl/ebl_franceschini.fits.gz" |
||
1415 | models["dominguez"] = "$GAMMAPY_DATA/ebl/ebl_dominguez11.fits.gz" |
||
1416 | models["finke"] = "$GAMMAPY_DATA/ebl/frd_abs.fits.gz" |
||
1417 | |||
1418 | return cls.read(models[name]) |
||
1419 | |||
1420 | def table_model(self, parameter, unit="TeV"): |
||
1421 | """Table model for a given parameter (`~gammapy.spectrum.models.TableModel`). |
||
1422 | |||
1423 | Parameters |
||
1424 | ---------- |
||
1425 | parameter : float |
||
1426 | Parameter value. |
||
1427 | unit : str, (optional) |
||
1428 | desired value for energy axis |
||
1429 | """ |
||
1430 | energy_axis = self.data.axes[1] |
||
1431 | energy = (energy_axis.log_center()).to(unit) |
||
1432 | |||
1433 | values = self.evaluate(energy=energy, parameter=parameter) |
||
1434 | |||
1435 | return TableModel(energy=energy, values=values, values_scale="lin") |
||
1436 | |||
1437 | def evaluate(self, energy, parameter): |
||
1438 | """Evaluate model for energy and parameter value.""" |
||
1439 | return self.data.evaluate(energy=energy, parameter=parameter) |
||
1440 | |||
1441 | |||
1442 | class AbsorbedSpectralModel(SpectralModel): |
||
1443 | """Spectral model with EBL absorption. |
||
1444 | |||
1445 | Parameters |
||
1446 | ---------- |
||
1447 | spectral_model : `~gammapy.spectrum.models.SpectralModel` |
||
1448 | Spectral model. |
||
1449 | absorption : `~gammapy.spectrum.models.Absorption` |
||
1450 | Absorption model. |
||
1451 | parameter : float |
||
1452 | parameter value for absorption model |
||
1453 | parameter_name : str, optional |
||
1454 | parameter name |
||
1455 | """ |
||
1456 | |||
1457 | def __init__( |
||
1458 | self, spectral_model, absorption, parameter, parameter_name="redshift" |
||
1459 | ): |
||
1460 | self.spectral_model = spectral_model |
||
1461 | self.absorption = absorption |
||
1462 | self.parameter = parameter |
||
1463 | self.parameter_name = parameter_name |
||
1464 | |||
1465 | # initialise list parameters from spectral model |
||
1466 | param_list = [] |
||
1467 | for param in spectral_model.parameters.parameters: |
||
1468 | param_list.append(param) |
||
1469 | |||
1470 | # Add parameter to the list |
||
1471 | min_ = self.absorption.data.axes[0].lo[0] |
||
1472 | max_ = self.absorption.data.axes[0].lo[-1] |
||
1473 | par = Parameter(parameter_name, parameter, min=min_, max=max_, frozen=True) |
||
1474 | param_list.append(par) |
||
1475 | |||
1476 | self.parameters = Parameters(param_list) |
||
1477 | |||
1478 | def evaluate(self, energy, **kwargs): |
||
1479 | """Evaluate the model at a given energy.""" |
||
1480 | # assign redshift value and remove it from dictionnary |
||
1481 | # since it does not belong to the spectral model |
||
1482 | parameter = kwargs[self.parameter_name] |
||
1483 | del kwargs[self.parameter_name] |
||
1484 | |||
1485 | flux = self.spectral_model.evaluate(energy=energy, **kwargs) |
||
1486 | absorption = self.absorption.evaluate(energy=energy, parameter=parameter) |
||
1487 | return flux * absorption |
||
1488 |