diff --git a/src/datasets/ogip.jl b/src/datasets/ogip.jl index 934cbb7b..dd64946a 100644 --- a/src/datasets/ogip.jl +++ b/src/datasets/ogip.jl @@ -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 @@ -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, @@ -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( @@ -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)) @@ -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, @@ -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) @@ -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 diff --git a/src/datasets/ogipdataset.jl b/src/datasets/ogipdataset.jl index d98d8d01..882c073a 100644 --- a/src/datasets/ogipdataset.jl +++ b/src/datasets/ogipdataset.jl @@ -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 diff --git a/src/datasets/spectraldata.jl b/src/datasets/spectraldata.jl index 90a1af89..63aff79e 100644 --- a/src/datasets/spectraldata.jl +++ b/src/datasets/spectraldata.jl @@ -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, @@ -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 diff --git a/test/datasets/test-ogip.jl b/test/datasets/test-ogip.jl index 175a0429..d20b2511 100644 --- a/test/datasets/test-ogip.jl +++ b/test/datasets/test-ogip.jl @@ -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 @@ -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 @@ -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" @@ -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