Fix nan/inf or negative values (if not log scale) for evaluate function in TableModel
@@ -27,6 +27,7 @@ __all__ = [
'ExponentialCutoffPowerLaw3FGL',
'LogParabola',
'TableModel',
+ 'AbsorbedSpectralModel',
]
@@ -588,6 +589,9 @@ class TableModel(SpectralModel):
self.values = values
self.scale_logy = scale_logy
+ self.lo_threshold = energy[0]
+ self.hi_threshold = energy[-1]
+
loge = np.log10(self.energy.to('eV').value)
try:
self.unit = self.values.unit
@@ -660,10 +664,38 @@ class TableModel(SpectralModel):
return cls(energy=energy, values=values, scale_logy=False)
def evaluate(self, energy, amplitude):
- interpy = self.interpy(np.log10(energy.to('eV').value))
+
+ is_array = True
+ try:
+ len(energy)
+ except:
+ is_array = False
+
+ # Not working for astropy.units.quantity.Quantity
+ # if isinstance(energy, (np.ndarray, np.generic)):
+ if is_array: # Test if array
+ # initialise array value to zero (dim energy)
+ values = np.zeros(len(energy), dtype=float)
+ # mask for energy range
+ mask = (energy >= self.lo_threshold) & (
+ energy <= self.hi_threshold)
+ # apply interpolation for masked values
+ values[mask] = self.interpy(np.log10(energy[mask].to('eV').value))
+ # Get rid of negative values (due to interpolation)
+ # Needed because of the rand.poisson used in SpectrumSimulation class
+ # Should be fixed in the class itself ?
+ if self.scale_logy is False:
+ values[values < 0] = 0.
+ else: # if not array
+ # test if energy is in range
+ if (energy >= self.lo_threshold or energy <= self.hi_threshold):
+ values = self.interpy(np.log10(energy.to('eV').value))
+ else:
+ values = 0
+
if self.scale_logy:
- interpy = np.power(10, interpy)
- return amplitude * interpy * self.unit
+ values = np.power(10, values)
+ return amplitude * values * self.unit
def plot(self, energy_range, ax=None, energy_unit='TeV',
n_points=100, **kwargs):
@@ -708,3 +740,33 @@ class TableModel(SpectralModel):
if self.scale_logy:
ax.set_yscale("log", nonposy='clip')
return ax
+
+
+class AbsorbedSpectralModel(SpectralModel):
+
+ def __init__(self, spectral_model, table_model):
+ """Absorbed spectral model
+
+ Parameters
+ ---------
+ spectral_model : `~gammapy.spectrum.models.SpectralModel`
+ spectral model
+ table_model : `~gammapy.spectrum.models.TableModel`
+ table model
+ """
+ self.spectral_model = spectral_model
+ self.table_model = table_model
+ # Will be implemented later for sherpa fit
+ self.parameters = {}
+
+ def evaluate(self, energy):
+ flux = self.spectral_model.__call__(energy)
+ absorption = self.table_model.__call__(energy)
+ return flux * absorption
+
+ def to_sherpa(self, name='default'):
+ """Convert to Sherpa model
+
+ To be implemented by subclasses
+ """
+ raise NotImplementedError('{}'.format(self.__class__.__name__))
