Merge pull request #81 from fjebaker/fergus/ogip-flexibility

More flexibility with OGIP datasets

Project.toml

@@ -1,7 +1,7 @@
 name = "SpectralFitting"
 uuid = "f2c56810-742e-4b72-8bf4-27af3bb81a12"
 authors = ["Fergus Baker <[email protected]>"]
-version = "0.5.4"
+version = "0.5.5"
 
 [deps]
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"

src/datasets/datasets.jl

@@ -169,8 +169,5 @@ export make_model_domain,
 include("spectrum.jl")
 include("response.jl")
 include("spectraldata.jl")
-
-# mission specifics
 include("ogipdataset.jl")
-include("xmm-newton.jl")
-include("nustart.jl")
+include("mission-specifics.jl")

src/datasets/mission-specifics.jl

@@ -0,0 +1,87 @@
+
+struct NuStarData{T,H} <: AbstractDataset
+    data::SpectralData{T}
+    paths::SpectralDataPaths
+    observation_id::String
+    exposure_id::String
+    object::String
+    header::H
+end
+
+NuStarData(spec_path; rmf_matrix_index = 3, rmf_energy_index = 2, kwargs...) = NuStarData(
+    load_ogip_dataset(
+        spec_path;
+        rmf_matrix_index = rmf_matrix_index,
+        rmf_energy_index = rmf_energy_index,
+        kwargs...,
+    )...,
+)
+
+make_label(data::NuStarData) = data.observation_id
+
+@_forward_SpectralData_api NuStarData.data
+
+function Base.show(io::IO, @nospecialize(data::NuStarData{T})) where {T}
+    print(io, "NuStarData[obs_id=$(data.observation_id)]")
+end
+
+function _printinfo(io, data::NuStarData{T}) where {T}
+    descr = """NuStarData:
+      . Object              : $(data.object)
+      . Observation ID      : $(data.observation_id)
+      . Exposure ID         : $(data.exposure_id)
+    """
+    print(io, descr)
+    _printinfo(io, data.data)
+end
+
+
+abstract type AbstractXmmNewtonDevice end
+struct XmmEPIC <: AbstractXmmNewtonDevice end
+
+struct XmmData{T,H,D} <: AbstractDataset
+    device::D
+    data::SpectralData{T}
+    paths::SpectralDataPaths
+    observation_id::String
+    exposure_id::String
+    object::String
+    header::H
+end
+
+XmmData(
+    device::AbstractXmmNewtonDevice,
+    spec_path;
+    rmf_matrix_index = 2,
+    rmf_energy_index = 3,
+    kwargs...,
+) = XmmData(
+    device,
+    load_ogip_dataset(
+        spec_path;
+        rmf_matrix_index = rmf_matrix_index,
+        rmf_energy_index = rmf_energy_index,
+        kwargs...,
+    )...,
+)
+
+make_label(data::XmmData) = data.observation_id
+
+@_forward_SpectralData_api XmmData.data
+
+function Base.show(io::IO, @nospecialize(data::XmmData{T})) where {T}
+    print(io, "XmmData[dev=$(data.device),obs_id=$(data.observation_id)]")
+end
+
+function _printinfo(io, data::XmmData{T}) where {T}
+    descr = """XmmData for $(Base.typename(typeof(data.device)).name):
+      . Object              : $(data.object)
+      . Observation ID      : $(data.observation_id)
+      . Exposure ID         : $(data.exposure_id)
+    """
+    print(io, descr)
+    _printinfo(io, data.data)
+end
+
+
+export NuStarData, XmmData, XmmEPIC

src/datasets/ogip.jl

@@ -14,13 +14,15 @@ 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)
+        new{T}(rmf_matrix_index, rmf_energy_index, rmf_matrix_key)
     end
 end
 rmf_matrix_index(c::StandardOGIPConfig) = c.rmf_matrix_index
@@ -60,6 +62,21 @@ function _parse_any(::Type{T}, @nospecialize(value::V))::T where {T,V}
     end
 end
 
+function _string_boolean(@nospecialize(value::V))::Bool where {V}
+    if V <: AbstractString
+        if value == "F"
+            false
+        elseif value == "T"
+            true
+        else
+            @warn("Unknown boolean string: $(value)")
+            false
+        end
+    else
+        value
+    end
+end
+
 # functions
 function parse_rmf_header(table::TableHDU)
     header = FITSIO.read_header(table)
@@ -107,6 +124,24 @@ function _chan_to_vectors(chan::Matrix)
     end
 end
 
+function _translate_channel_array(channel)
+    if channel isa Matrix
+        _chan_to_vectors(channel)
+    elseif eltype(channel) <: AbstractVector
+        channel
+    else
+        map(i -> [i], channel)
+    end
+end
+
+function _adapt_matrix_type(T::Type, mat::M) where {M}
+    if eltype(M) <: AbstractVector
+        map(row -> convert.(T, row), mat)
+    elseif M <: AbstractMatrix
+        map(row -> convert.(T, row), eachcol(mat))
+    end
+end
+
 function read_rmf_matrix(table::TableHDU, header::RMFHeader, T::Type)
     energy_low = convert.(T, read(table, "ENERG_LO"))
     energy_high = convert.(T, read(table, "ENERG_HI"))
@@ -115,17 +150,18 @@ function read_rmf_matrix(table::TableHDU, header::RMFHeader, T::Type)
     matrix_raw = read(table, "MATRIX")
 
     # type stable: convert to common vector of vector format
-    f_chan::Vector{Vector{Int}} =
-        f_chan_raw isa Matrix ? _chan_to_vectors(f_chan_raw) : f_chan_raw
-    n_chan::Vector{Vector{Int}} =
-        n_chan_raw isa Matrix ? _chan_to_vectors(n_chan_raw) : n_chan_raw
+    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,
         energy_low,
         energy_high,
-        map(row -> convert.(T, row), matrix_raw),
+        _adapt_matrix_type(T, matrix_raw),
         header,
     )
 end
@@ -159,7 +195,7 @@ end
 function read_ancillary_response(path::String, ogip_config::AbstractOGIPConfig{T}) where {T}
     fits = FITS(path)
     (bins_low, bins_high, effective_area) = _read_fits_and_close(path) do fits
-        area::Vector{T} = convert.(T, read(fits[2], "SPECRESP"))
+        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"))
         (lo, hi, area)
@@ -229,7 +265,7 @@ function read_spectrum(path, config::AbstractOGIPConfig{T}) where {T}
     info = _read_fits_and_close(path) do fits
         header = read_header(fits[2])
         # if not set, assume not poisson errors
-        is_poisson = get(header, "POISSERR", false)
+        is_poisson = _string_boolean(get(header, "POISSERR", false))
         # read general infos
         instrument = header["INSTRUME"]
         telescope = header["TELESCOP"]
@@ -314,6 +350,9 @@ function read_filename(header, entry, parent, exts...)
     parent_name = basename(parent)
     if haskey(header, entry)
         path::String = strip(header[entry])
+        if path == "NONE"
+            return missing
+        end
         name = find_file(data_directory, path, parent_name, exts)
         if !ismissing(name)
             return name

src/datasets/ogipdataset.jl

@@ -7,9 +7,8 @@ struct OGIPDataset{T,H} <: AbstractDataset
     header::H
 end
 
-function OGIPDataset(
+function load_ogip_dataset(
     spec_path;
-    T::Type = Float64,
     hdu = 2,
     background = missing,
     response = missing,
@@ -22,20 +21,24 @@ function OGIPDataset(
         response = response,
         ancillary = ancillary,
     )
-    config = StandardOGIPConfig(rmf_matrix_index = 2, rmf_energy_index = 3, T = T)
+    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; kwargs...)
-    OGIPDataset(data, paths, obs_id, exposure_id, object, header)
+    data = SpectralData(paths, config)
+    (data, paths, obs_id, exposure_id, object, header)
 end
 
+OGIPDataset(spec_path; kwargs...) = OGIPDataset(load_ogip_dataset(spec_path; kwargs...)...)
+
 @_forward_SpectralData_api OGIPDataset.data
 
+make_label(d::OGIPDataset) = d.observation_id
+
 function _printinfo(io, data::OGIPDataset{T}) where {T}
     descr = """OGIPDataset:
       . Object              : $(data.object)

src/datasets/spectraldata.jl

@@ -260,8 +260,9 @@ function _dataset_from_ogip(paths::SpectralDataPaths, config::OGIP.AbstractOGIPC
     resp = if !ismissing(paths.response)
         OGIP.read_rmf(paths.response, config)
     else
-        @warn "No response file found."
-        missing
+        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)