Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

make RhessiLoader compatible with legacy fitter args #181

Draft
wants to merge 5 commits into
base: main
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
114 changes: 81 additions & 33 deletions sunkit_spex/extern/rhessi.py
Original file line number Diff line number Diff line change
Expand Up @@ -28,16 +28,15 @@ class RhessiLoader(instruments.InstrumentBlueprint):
[1] https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/XSappendixStatistics.html
"""

def __init__(self, spectrum_fn, srm_fn, **kwargs):
def __init__(self, spectrum_fn, srm_fn):
"""
Spectrum and SRM files are both required: attenuator state change times
are in the spectrum file,
and the state determines which SRM will be used.
"""
self._construction_string = f"RhessiLoader(spectrum_fn={spectrum_fn}, " f"srm_fn={srm_fn}," f"**{kwargs})"
self._construction_string = f"RhessiLoader(spectrum_fn={spectrum_fn}, " f"srm_fn={srm_fn})"
self._systematic_error = 0
self.load_prepare_spectrum_srm(spectrum_fn, srm_fn)
self._start_background_time, self._end_background_time = None, None

@property
def systematic_error(self):
Expand Down Expand Up @@ -96,9 +95,11 @@ def load_prepare_spectrum_srm(self, spectrum_fn, srm_fn):
}
self._update_event_data_with_times()
self._original_data = copy.deepcopy(self._loaded_spec_data)
# Set these times to sensible defaults
self.update_background_times(self._start_event_time, self._end_event_time)

@property
def subtract_background(self):
def subtract_background(self) -> bool:
"""
States whether the the data is event minus background or not.
"""
Expand Down Expand Up @@ -445,10 +446,9 @@ def lightcurve(self, energy_ranges=None, axes=None, rebin_time=1):
elif len(np.shape(energy_ranges)) == 1:
energy_ranges = [energy_ranges]
elif len(np.shape(energy_ranges)) == 0:
print(
raise ValueError(
"The `energy_ranges` input should be a range of two energy values (e.g., [4,8] meaning 4-8 keV inclusive) or a list of these ranges."
)
return

ax = axes or plt.gca()
font_size = plt.rcParams["font.size"]
Expand All @@ -462,11 +462,11 @@ def lightcurve(self, energy_ranges=None, axes=None, rebin_time=1):

# plot each energy range
time_binning = (flat_times[1:] - flat_times[:-1]).to_value("s")
left_edges, right_edges = self._loaded_spec_data["count_channel_bins"].T
for er in energy_ranges:
i = np.where(
(self._loaded_spec_data["count_channel_bins"][:, 0] >= er[0])
& (self._loaded_spec_data["count_channel_bins"][:, -1] <= er[-1])
)
# Restrict the energy range to the data energy bins
ea, eb = max(er[0], left_edges[0]), min(er[1], right_edges[-1])
i = np.where((left_edges >= ea) & (right_edges <= eb))
e_range_cts = np.sum(self._spectrum["counts"][:, i].reshape((len(time_binning), -1)), axis=1)

if rebin_time > 1:
Expand All @@ -475,7 +475,7 @@ def lightcurve(self, energy_ranges=None, axes=None, rebin_time=1):

e_range_ctr, e_range_ctr_err = e_range_cts / time_binning, np.sqrt(e_range_cts) / time_binning
lc = e_range_ctr
p = ax.stairs(lc, flat_times.datetime, label=f"{er[0]}$-${er[-1]} keV")
p = ax.stairs(lc, flat_times.datetime, label=f"{ea}$-${eb} keV")
ax.errorbar(
x=time_mids.datetime, y=e_range_ctr, yerr=e_range_ctr_err, c=p.get_edgecolor(), ls=""
) # error bar in middle of the bin
Expand Down Expand Up @@ -604,32 +604,65 @@ def spectrogram(self, axes=None, rebin_time=1, rebin_energy=1, **kwargs):

return ax

def select_time(self, start=None, end=None, background=False):
"""Provides method to set start and end time of the event or background in one line.
# "retroactive" properties for setting the times like in older versions,
# for legacy compatibility.
@property
def start_event_time(self):
time_depr_warn('event')
return self._start_event_time

Parameters
----------
start, end : str, `astropy.Time`, None
String to be given to astropy's Time, `astropy.Time` is used directly, None sets the
start/end event time to be the first time of the data. None doesn't add, or will remove,
any background data in `_loaded_spec_data["extras"]` if background is set to True.
Default: None
@start_event_time.setter
def start_event_time(self, new_time: atime.Time):
time_depr_warn('event')
self.update_event_times(new_time, self._end_event_time)

background : bool
Determines whether the start and end times are for the event (False) or background
time (True).
Default: False
@property
def end_event_time(self):
time_depr_warn('event')
return self._end_event_time

Returns
-------
None.
"""
self.__warn = False if (start is not None) and (end is not None) else True
@end_event_time.setter
def end_event_time(self, new_time: atime.Time):
time_depr_warn('event')
self.update_event_times(self._start_event_time, new_time)

if background:
self.start_background_time, self.end_background_time = start, end
else:
self.start_event_time, self.end_event_time = start, end
@property
def start_background_time(self):
time_depr_warn('background')
return self._start_background_time

@start_background_time.setter
def start_background_time(self, new_time: atime.Time):
time_depr_warn('background')
self.update_background_times(new_time, self._end_background_time)

@property
def end_background_time(self):
time_depr_warn('background')
return self._end_background_time

@end_background_time.setter
def end_background_time(self, new_time: atime.Time):
time_depr_warn('background')
self.update_background_times(self._start_background_time, new_time)

# "retroactive" properties for specifying background subtraction,
# like in older versions.
@property
def data2data_minus_background(self) -> bool:
warnings.warn(
"Please use the `subtract_background` property instead.",
stacklevel=2
)
return self.subtract_background

@data2data_minus_background.setter
def data2data_minus_background(self, do_subtract: bool):
warnings.warn(
"Please use the `subtract_background` property instead.",
stacklevel=2
)
self.subtract_background = do_subtract


def load_spectrum(spec_fn: str):
Expand Down Expand Up @@ -801,3 +834,18 @@ def load_srm(srm_file: str):
}

return dict(channel_bins=channel_bins, photon_bins=photon_bins, srm_options=ret_srms)


class LegacyRhessiLoader(RhessiLoader):
'''Allow legacy `Fitter` to load in RHESSI data with the expected args format'''
def __init__(self, spec_fn, **kw):
srm_fn = kw.pop('srm_file')
super().__init__(spec_fn, srm_fn)


def time_depr_warn(partial_name: str):
warnings.warn(
"This time-related property is deprecated. "
f"In the future, please use `RhessiLoader.update_{partial_name}_times`",
stacklevel=2
)
2 changes: 1 addition & 1 deletion sunkit_spex/legacy/fitting/data_loader.py
Original file line number Diff line number Diff line change
Expand Up @@ -136,7 +136,7 @@ def __init__(
f"srm_file={srm_file},srm_custom={srm_custom},custom_channel_bins={custom_channel_bins}, **{kwargs})"
)

self.instrument_loaders = {"NuSTAR": inst.NustarLoader, "RHESSI": rhessi.RhessiLoader, "Solar Orbiter": stix.STIXLoader}
self.instrument_loaders = {"NuSTAR": inst.NustarLoader, "RHESSI": rhessi.LegacyRhessiLoader, "Solar Orbiter": stix.STIXLoader}

pha_file, arf_file, rmf_file, srm_file, srm_custom, custom_channel_bins, instruments = self._sort_files(
pha_file=pha_file,
Expand Down
Loading