Skip to content

Commit

Permalink
Merge pull request #93 from fjebaker/fergus/rm-ogip-config
Browse files Browse the repository at this point in the history 
Remove OGIP config struct
  • Loading branch information
@fjebaker
fjebaker authored May 8, 2024
2 parents b38142d + 147f5d1 commit 3b8e65c
Show file tree
Hide file tree
Showing 4 changed files with 61 additions and 68 deletions.
85 changes: 41 additions & 44 deletions src/datasets/ogip.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,29 +5,6 @@ import SpectralFitting: SpectralUnits
using FITSIO
using SparseArrays

export AbstractOGIPConfig, StandardOGIPConfig

abstract type AbstractOGIPConfig{T} end
rmf_matrix_index(c::AbstractOGIPConfig) = error("Not implemented for $(typeof(c))")
rmf_energy_index(c::AbstractOGIPConfig) = error("Not implemented for $(typeof(c))")

struct StandardOGIPConfig{T} <: AbstractOGIPConfig{T}
rmf_matrix_index::Int
rmf_energy_index::Int
rmf_matrix_key::String
function StandardOGIPConfig(;
rmf_matrix_index = 3,
rmf_energy_index = 2,
rmf_matrix_key = "SPECRESP",
T::Type = Float64,
)
@assert rmf_matrix_index != rmf_energy_index
new{T}(rmf_matrix_index, rmf_energy_index, rmf_matrix_key)
end
end
rmf_matrix_index(c::StandardOGIPConfig) = c.rmf_matrix_index
rmf_energy_index(c::StandardOGIPConfig) = c.rmf_energy_index

struct MissingHeader <: Exception
header::String
end
Expand Down Expand Up @@ -153,9 +130,6 @@ function read_rmf_matrix(table::TableHDU, header::RMFHeader, T::Type)
f_chan::Vector{Vector{Int}} = _translate_channel_array(f_chan_raw)
n_chan::Vector{Vector{Int}} = _translate_channel_array(n_chan_raw)

@show typeof(matrix_raw)
@show size(f_chan), size(matrix_raw)

RMFMatrix(
f_chan,
n_chan,
Expand All @@ -178,29 +152,52 @@ function _read_fits_and_close(f, path)
res
end

function read_rmf(path::String, ogip_config::AbstractOGIPConfig{T}) where {T}
i = rmf_matrix_index(ogip_config)
j = rmf_energy_index(ogip_config)

function read_rmf(path::String; T::Type = Float64)
(header, rmf, channels::RMFChannels{T}) = _read_fits_and_close(path) do fits
hdr = parse_rmf_header(fits[i])
_rmf = read_rmf_matrix(fits[i], hdr, T)
_channels = read_rmf_channels(fits[j], T)
rmf_i = find_extension(fits, ["RESP", "MATRIX"])
energy_i = find_extension(fits, "EBOUND")

hdr = parse_rmf_header(fits[rmf_i])
_rmf = read_rmf_matrix(fits[rmf_i], hdr, T)
_channels = read_rmf_channels(fits[energy_i], T)
(hdr, _rmf, _channels)
end

_build_reponse_matrix(header, rmf, channels, T)
end

function read_ancillary_response(path::String, ogip_config::AbstractOGIPConfig{T}) where {T}
function find_extension(fits, extension::T) where {T <: Union{<:AbstractString, <:AbstractVector}}
# find the correct extensions
i::Int = 1
for hdu in fits
header = read_header(hdu)
extname = get(header, "EXTNAME", "")
if T <: AbstractString
if contains(extname, extension)
return i
end
elseif T <: AbstractVector
for ext in extension
if contains(extname, ext)
return i
end
end
end
i += 1
end
return nothing
end

function read_ancillary_response(path::String; T::Type = Float64)
fits = FITS(path)
(bins_low, bins_high, effective_area) = _read_fits_and_close(path) do fits
area::Vector{T} = convert.(T, read(fits[2], ogip_config.rmf_matrix_key))
lo::Vector{T} = convert.(T, read(fits[2], "ENERG_LO"))
hi::Vector{T} = convert.(T, read(fits[2], "ENERG_HI"))
i = find_extension(fits, "RESP")
area::Vector{T} = convert.(T, read(fits[i], "SPECRESP"))
lo::Vector{T} = convert.(T, read(fits[i], "ENERG_LO"))
hi::Vector{T} = convert.(T, read(fits[i], "ENERG_HI"))
(lo, hi, area)
end
SpectralFitting.AncillaryResponse(bins_low, bins_high, effective_area)
SpectralFitting.AncillaryResponse{T}(bins_low, bins_high, effective_area)
end

function _build_reponse_matrix(
Expand Down Expand Up @@ -261,8 +258,8 @@ function _get_stable(::Type{T}, header, name, default)::T where {T}
get(header, name, T(default))
end

function read_spectrum(path, config::AbstractOGIPConfig{T}) where {T}
info = _read_fits_and_close(path) do fits
function read_spectrum(path; T::Type = Float64)
info::SpectralFitting.Spectrum{T} = _read_fits_and_close(path) do fits
header = read_header(fits[2])
# if not set, assume not poisson errors
is_poisson = _string_boolean(get(header, "POISSERR", false))
Expand Down Expand Up @@ -309,7 +306,7 @@ function read_spectrum(path, config::AbstractOGIPConfig{T}) where {T}
SpectralFitting.ErrorStatistics.Unknown, T[0 for _ in values]
end

SpectralFitting.Spectrum(
SpectralFitting.Spectrum{T}(
channels,
quality,
grouping,
Expand All @@ -328,8 +325,8 @@ function read_spectrum(path, config::AbstractOGIPConfig{T}) where {T}
info
end

function read_background(path::String, config::AbstractOGIPConfig{T}) where {T}
read_spectrum(path, config)
function read_background(path::String)
read_spectrum(path)
end

function read_paths_from_spectrum(path::String)
Expand Down Expand Up @@ -383,7 +380,7 @@ end
end # module

using .OGIP
export OGIP, StandardOGIPConfig
export OGIP

function read_fits_header(path; hdu = 2)
OGIP._read_fits_and_close(path) do f
Expand Down
4 changes: 1 addition & 3 deletions src/datasets/ogipdataset.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,15 +21,13 @@ function load_ogip_dataset(
response = response,
ancillary = ancillary,
)
config = StandardOGIPConfig(; kwargs...)

header = read_fits_header(paths.spectrum; hdu = hdu)

obs_id = haskey(header, "OBS_ID") ? header["OBS_ID"] : "[no observation id]"
exposure_id = haskey(header, "EXP_ID") ? header["EXP_ID"] : "[no exposure id]"
object = haskey(header, "OBJECT") ? header["OBJECT"] : "[no object]"

data = SpectralData(paths, config)
data = SpectralData(paths)
(data, paths, obs_id, exposure_id, object, header)
end

Expand Down
14 changes: 7 additions & 7 deletions src/datasets/spectraldata.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,8 +57,8 @@ end

# constructor

SpectralData(paths::SpectralDataPaths, config::OGIP.AbstractOGIPConfig; kwargs...) =
_dataset_from_ogip(paths, config; kwargs...)
SpectralData(paths::SpectralDataPaths; kwargs...) =
_dataset_from_ogip(paths; kwargs...)

function SpectralData(
spectrum::Spectrum,
Expand Down Expand Up @@ -249,23 +249,23 @@ function rebin_if_different_domains!(output, data_domain, model_domain, input)
output
end

function _dataset_from_ogip(paths::SpectralDataPaths, config::OGIP.AbstractOGIPConfig)
spec = OGIP.read_spectrum(paths.spectrum, config)
function _dataset_from_ogip(paths::SpectralDataPaths)
spec = OGIP.read_spectrum(paths.spectrum)
back = if !ismissing(paths.background)
OGIP.read_background(paths.background, config)
OGIP.read_background(paths.background)
else
@warn "No background file found."
missing
end
resp = if !ismissing(paths.response)
OGIP.read_rmf(paths.response, config)
OGIP.read_rmf(paths.response)
else
throw(
"No response file found in the header. Response must be specified with the keyword `response=PATH`.",
)
end
ancillary = if !ismissing(paths.ancillary)
OGIP.read_ancillary_response(paths.ancillary, config)
OGIP.read_ancillary_response(paths.ancillary)
else
@warn "No ancillary file found."
missing
Expand Down
26 changes: 12 additions & 14 deletions test/datasets/test-ogip.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,16 +3,14 @@ using SpectralFitting

@show testdir

xmm_config = StandardOGIPConfig(rmf_matrix_index = 2, rmf_energy_index = 3)
nustar_config = StandardOGIPConfig(rmf_matrix_index = 3, rmf_energy_index = 2)

xmm_rmf_path = joinpath(testdir, "xmm/pn.rmf")
nustar_rmf_path = joinpath(testdir, "nustar/nu60001047002A01_sr.rmf")

# test we can read in the different response matrices correctly
rmf_xmm = OGIP.read_rmf(xmm_rmf_path, xmm_config)
@inferred OGIP.read_rmf(xmm_rmf_path, xmm_config)
rmf_nustar = OGIP.read_rmf(nustar_rmf_path, nustar_config)
rmf_xmm = OGIP.read_rmf(xmm_rmf_path)
@inferred OGIP.read_rmf(xmm_rmf_path)
rmf_nustar = OGIP.read_rmf(nustar_rmf_path)

# make sure the correct number of values have been written in
@test count(!=(0), rmf_xmm.matrix) == 928061
Expand All @@ -22,9 +20,9 @@ rmf_nustar = OGIP.read_rmf(nustar_rmf_path, nustar_config)
xmm_arf_path = joinpath(testdir, "xmm/pn.arf")
nustar_arf_path = joinpath(testdir, "nustar/nu60001047002A01_sr.arf")

arf_xmm = OGIP.read_ancillary_response(xmm_arf_path, xmm_config)
@inferred OGIP.read_ancillary_response(xmm_arf_path, xmm_config)
arf_nustar = OGIP.read_ancillary_response(nustar_arf_path, nustar_config)
arf_xmm = OGIP.read_ancillary_response(xmm_arf_path)
@inferred OGIP.read_ancillary_response(xmm_arf_path)
arf_nustar = OGIP.read_ancillary_response(nustar_arf_path)

@test arf_xmm.bins_low ≈ rmf_xmm.bins_low atol = 1e-7
@test arf_nustar.bins_low ≈ rmf_nustar.bins_low atol = 1e-3
Expand All @@ -33,9 +31,9 @@ arf_nustar = OGIP.read_ancillary_response(nustar_arf_path, nustar_config)
xmm_spec_path = joinpath(testdir, "xmm/pn_spec_grp.fits")
nustar_spec_path = joinpath(testdir, "nustar/nu60001047002A01_sr_grp_simple.pha")

spec_xmm = OGIP.read_spectrum(xmm_spec_path, xmm_config)
@inferred OGIP.read_spectrum(xmm_spec_path, xmm_config)
spec_nustar = OGIP.read_spectrum(nustar_spec_path, nustar_config)
spec_xmm = OGIP.read_spectrum(xmm_spec_path)
@inferred OGIP.read_spectrum(xmm_spec_path)
spec_nustar = OGIP.read_spectrum(nustar_spec_path)

@test spec_xmm.telescope_name == "XMM"
@test spec_nustar.telescope_name == "NuSTAR"
Expand All @@ -45,11 +43,11 @@ xmm_backgroud_path = joinpath(testdir, "xmm/pn_bg_spec.fits")
nustar_background_path = joinpath(testdir, "nustar/nu60001047002A01_bk.pha")


bg_xmm = OGIP.read_background(xmm_backgroud_path, xmm_config)
bg_xmm = OGIP.read_background(xmm_backgroud_path)
@test length(bg_xmm.channels) == 4096

bg_nustar = OGIP.read_background(nustar_background_path, nustar_config)
@inferred OGIP.read_background(nustar_background_path, nustar_config)
bg_nustar = OGIP.read_background(nustar_background_path)
@inferred OGIP.read_background(nustar_background_path)
@test length(bg_xmm.channels) == 4096

# test reading the associated paths
Expand Down

0 comments on commit 3b8e65c

Please sign in to comment.