search

USRA-STI / gdt-core

Gamma-ray Data Tools - Core Components
Apache License 2.0
7 stars 11 forks source link

Cannot set energy or channel range of simulated PHA data when converting to gdt.core.pha.PHA #66

Open rachel-hamburg opened 1 month ago

rachel-hamburg commented 1 month ago

When converting simulated PHA data from a gdt.core.simulate.pha.PhaSimulator object to a gdt.core.pha.PHA object using gdt.core.simulate.pha.PhaSimulator.to_pha, setting the energy_range attribute does not update the pha.energy_range or the pha.channel_mask.

Code to produce the issue:

from gdt.missions.fermi.gbm.tte import GbmTte
from gdt.core.binning.unbinned import bin_by_time
from gdt.missions.fermi.gbm.response import GbmRsp
from gdt.core.background.binned import Polynomial
from gdt.core.background.fitter import BackgroundFitter
from gdt.core.spectra.functions import Band
from gdt.core.simulate.pha import PhaSimulator

# Open TTE file
tte = GbmTte.open('data/150506630/glg_tte_n9_bn150506630_v00.fit')

# Open RSP file
rsp = GbmRsp.open('data/150506630/glg_cspec_n9_bn150506630_v02.rsp')

# Bin the TTE data by time and convert to PHAII
phaii = tte.to_phaii(bin_by_time, 1.024, time_ref=0.0)

# Fit the background
bkgd_fitter = BackgroundFitter.from_phaii(phaii, Polynomial, time_ranges=[(-22.17, -4.),( 2., 15.17)])
bkgd_fitter.fit(order=1)
bkgd = bkgd_fitter.interpolate_bins(phaii.data.tstart, phaii.data.tstop)

# Integrate background model over some time
model_bkgd_time = (11.,13.)
spec_bkgd_model = bkgd.integrate_time(*model_bkgd_time)

# Define source selection
exposure = 1.024

# Draw from parameters in the spectral catalog
model_params = (0.1, 300.0, -1., -2.)

# Simulate the source and background
pha_sim = PhaSimulator(rsp, Band(), model_params, exposure, spec_bkgd_model, 'Gaussian')

# Convert to PHA files
pha = pha_sim.to_pha(1, energy_range=(8.,900.))[0]  # apply energy range
print (pha.energy_range)                            # energy range is not applied
print (pha.channel_mask)                            # channel mask is not updated
print (len(pha.channel_mask))                       # there are still 128 channels so no channels have been cut

Output:

(4.3979997634887695, 2000.0)
[ True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True False  True  True  True]
128

If I slice the PHA object in energy and then try to fit the spectra, I receive an Index Error due to the number of channels removed from the full 128 channel range:

from gdt.core.spectra.fitting import SpectralFitterPstat

# Slice the PHA in energy
pha = pha.slice_energy((8.,900.))
print (pha.energy_range)
print (pha.channel_mask)
print (len(pha.channel_mask))

# Fit the spectrum
specfitter = SpectralFitterPstat([pha], [spec_bkgd_model], [rsp], method='TNC')
specfitter.fit(Band(), options={'maxiter': 1000})

Output:

(7.146999835968018, 911.9829711914062)
[ True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True False]
122

Traceback (most recent call last):
  File "/Users/hamburg/Dropbox/Work/GBM/batools_project/spec-catalog/catfiles/spectral_analysis/test2.py", line 50, in <module>
    specfitter = SpectralFitterPstat([pha], [spec_bkgd_model], [rsp], method='TNC')
  File "/Users/hamburg/venv/gdt-devel/lib/python3.10/site-packages/gdt/core/spectra/fitting.py", line 1072, in __init__
    super().__init__(pha_list, bkgd_list, rsp_list, pstat, **kwargs)
  File "/Users/hamburg/venv/gdt-devel/lib/python3.10/site-packages/gdt/core/spectra/fitting.py", line 140, in __init__
    self._back_rates = self._apply_masks(self._back_rates)
  File "/Users/hamburg/venv/gdt-devel/lib/python3.10/site-packages/gdt/core/spectra/fitting.py", line 725, in _apply_masks
    return [np.asarray(one_list)[one_mask] for one_list, one_mask in zip(a_list, self._chan_masks)]
  File "/Users/hamburg/venv/gdt-devel/lib/python3.10/site-packages/gdt/core/spectra/fitting.py", line 725, in <listcomp>
    return [np.asarray(one_list)[one_mask] for one_list, one_mask in zip(a_list, self._chan_masks)]
IndexError: boolean index did not match indexed array along dimension 0; dimension is 128 but corresponding boolean dimension is 122
rachel-hamburg commented 1 month ago

Ah I see related to Issue #50