Remove OGIP config struct

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

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
 

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

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