Mission Support - follow up

src/Stingray.jl

@@ -77,4 +77,13 @@ export create_filtered_lightcurve
 export check_gtis
 export split_by_gtis
 
+include("missionSupport.jl")
+export MissionSupport,
+       get_mission_support,
+       apply_calibration,
+       patch_mission_info,
+       SIMPLE_CALIBRATION_FUNCS,
+       interpret_fits_data!,
+       AbstractMissionSupport
+
 end

src/events.jl

@@ -1031,6 +1031,9 @@ Times are ready for direct use in timing analysis without additional corrections
 """
 function readevents(
     path::AbstractString;
+    mission::Union{String,Nothing} = nothing,
+    instrument::Union{String,Nothing} = nothing,
+    epoch::Union{Float64,Nothing} = nothing,
     hdu::Int = 2,
     T::Type = Float64,
     sort::Bool = false,
@@ -1045,6 +1048,16 @@ function readevents(
     kwargs...,
 )::EventList{Vector{T},FITSMetadata{FITSIO.FITSHeader}}
 
+    # Get mission support if specified
+    mission_support = if !isnothing(mission)
+        ms = get_mission_support(mission, instrument, epoch)
+        # Use mission-specific energy alternatives if available
+        energy_alternatives = ms.energy_alternatives
+        ms
+    else
+        nothing
+    end
+
     # Read GTI first if requested
     gti_data, gti_source = nothing, nothing
     if load_gti
@@ -1072,6 +1085,12 @@ function readevents(
     time_del::Union{Nothing,Float64} = FITS(path, "r") do f
 
         selected_hdu = f[hdu]
+
+        # Apply mission-specific FITS interpretation if available
+        if !isnothing(mission_support) && !isnothing(mission_support.interpretation_func)
+            interpret_fits_data!(f, mission_support)
+        end
+
         header = read_header(selected_hdu)
         all_cols = FITSIO.colnames(selected_hdu)
         time = convert(Vector{T}, read(selected_hdu, "TIME", case_sensitive = false))
@@ -1083,6 +1102,14 @@ function readevents(
             energy_alternatives = energy_alternatives,
         )
 
+        # Apply mission-specific calibration if we have PI data and mission support
+        if !isnothing(energy) && !isnothing(mission_support) && 
+           !isnothing(energy_column) && uppercase(energy_column) == "PI"
+            energy = convert(Vector{T}, apply_calibration(mission_support, energy))
+            # Update the energy column name to reflect that it's now calibrated
+            energy_column = "ENERGY"
+        end
+
         # Read extra columns
         extra_data = Dict{String,Vector}()
         for col_name in extra_columns
@@ -1101,6 +1128,23 @@ function readevents(
         (time, energy, energy_column, header, extra_data, mjd_ref, time_zero, time_unit, time_sys, time_pixr, time_del)
     end
 
+    # Apply mission-specific header patches if available
+    if !isnothing(mission_support)
+        # Convert header to dictionary for patching
+        header_dict = Dict{String,Any}()
+
+        # Use the proper way to access FITSHeader keys and values
+        for key in keys(header)
+            header_dict[key] = header[key]
+        end
+
+        # Apply mission patches
+        patched_header_dict = patch_mission_info(header_dict, mission)
+
+        # Note: We keep the original header structure but could extend this
+        # to update the header with patched information if needed
+    end
+
     # Apply time corrections and optionally convert to MJD
     if !isnothing(mjd_ref)
         effective_timezero = isnothing(time_zero) ? 0.0 : time_zero
@@ -1206,14 +1250,8 @@ function Base.summary(ev::EventList)
     if has_energies(ev)
         energy_range = extrema(ev.energies)
         summary_str *= ", energies: $(energy_range[1]) - $(energy_range[2])"
-        # Map energy units to display units
         if !isnothing(ev.meta.energy_units)
-            display_units = if ev.meta.energy_units in ["ENERGY", "PI", "PHA"]
-                "keV"
-            else
-                ev.meta.energy_units
-            end
-            summary_str *= " ($display_units)"
+            summary_str *= " ($(ev.meta.energy_units))"
         end
     end
     if has_gti(ev)

src/gti.jl

@@ -97,13 +97,13 @@ function check_gtis(gti::AbstractMatrix)
     gti_start = @view gti[:, 1]
     gti_end = @view gti[:, 2]
 
-    if any(gti_end < gti_start)
+    if any(gti_end .< gti_start)
         throw(ArgumentError(
             "The GTI end times must be larger than the GTI start times."
         )) 
     end
 
-    if any(@view(gti_start[begin+1:end]) < @view(gti_end[begin:end-1]))
+    if any(@view(gti_start[begin+1:end]) .< @view(gti_end[begin:end-1]))
         throw(ArgumentError(
             "This GTI has overlaps"
         ))

src/missionSupport.jl

@@ -0,0 +1,226 @@
+"""
+Dictionary of simple conversion functions for different missions.
+
+This dictionary provides PI (Pulse Invariant) to energy conversion functions
+for various X-ray astronomy missions. Each function takes a PI channel value
+and returns the corresponding energy in keV.
+
+Supported missions:
+- NuSTAR: Nuclear Spectroscopic Telescope Array
+- XMM: X-ray Multi-Mirror Mission  
+- NICER: Neutron star Interior Composition Explorer
+- IXPE: Imaging X-ray Polarimetry Explorer
+- AXAF/Chandra: Advanced X-ray Astrophysics Facility
+- XTE/RXTE: Rossi X-ray Timing Explorer
+"""
+const SIMPLE_CALIBRATION_FUNCS = Dict{String, Function}(
+    "nustar" => (pi) -> pi * 0.04 + 1.62,
+    "xmm" => (pi) -> pi * 0.001,
+    "nicer" => (pi) -> pi * 0.01,
+    "ixpe" => (pi) -> pi / 375 * 15,
+    "axaf" => (pi) -> (pi - 1) * 14.6e-3,  # Chandra/AXAF
+    "chandra" => (pi) -> (pi - 1) * 14.6e-3,  # Explicit chandra entry
+    "xte" => (pi) -> pi * 0.025  # RXTE/XTE
+)
+
+"""
+Abstract type for mission-specific calibration and interpretation.
+
+This serves as the base type for all mission support implementations,
+allowing for extensibility and type safety in mission-specific operations.
+"""
+abstract type AbstractMissionSupport end
+
+"""
+    MissionSupport{T} <: AbstractMissionSupport
+
+Structure containing mission-specific calibration and interpretation information.
+
+This structure encapsulates all the necessary information for handling
+data from a specific X-ray astronomy mission, including calibration
+functions, energy column alternatives, and GTI extension preferences.
+
+# Fields
+- `name::String`: Mission name (normalized to lowercase)
+- `instrument::Union{String, Nothing}`: Instrument identifier
+- `epoch::Union{T, Nothing}`: Observation epoch in MJD (for time-dependent calibrations)
+- `calibration_func::Function`: PI to energy conversion function
+- `interpretation_func::Union{Function, Nothing}`: Mission-specific FITS interpretation function
+- `energy_alternatives::Vector{String}`: Preferred energy column names in order of preference
+- `gti_extensions::Vector{String}`: GTI extension names in order of preference
+
+# Type Parameters
+- `T`: Type of the epoch parameter (typically Float64)
+"""
+struct MissionSupport{T} <: AbstractMissionSupport
+    name::String
+    instrument::Union{String, Nothing}
+    epoch::Union{T, Nothing}
+    calibration_func::Function
+    interpretation_func::Union{Function, Nothing}
+    energy_alternatives::Vector{String}
+    gti_extensions::Vector{String}
+end
+
+"""
+    get_mission_support(mission::String, instrument=nothing, epoch=nothing) -> MissionSupport
+
+Create mission support object with mission-specific parameters.
+
+This function creates a MissionSupport object containing all the necessary
+information for processing data from a specified X-ray astronomy mission.
+It handles mission aliases (e.g., Chandra/AXAF) and provides appropriate
+defaults for each mission.
+
+# Arguments
+- `mission::String`: Mission name (case-insensitive)
+- `instrument::Union{String, Nothing}=nothing`: Instrument identifier
+- `epoch::Union{Float64, Nothing}=nothing`: Observation epoch in MJD
+
+# Returns
+- `MissionSupport{Float64}`: Mission support object
+
+# Throws
+- `ArgumentError`: If mission name is empty
+
+# Examples
+```julia
+# Basic usage
+ms = get_mission_support("nustar")
+
+# With instrument specification
+ms = get_mission_support("nustar", "FPM_A")
+
+# With epoch for time-dependent calibrations
+ms = get_mission_support("xte", "PCA", 50000.0)
+```
+"""
+function get_mission_support(mission::String, 
+                           instrument::Union{String, Nothing}=nothing,
+                           epoch::Union{Float64, Nothing}=nothing)
+
+    # Check for empty mission string
+    if isempty(mission)
+        throw(ArgumentError("Mission name cannot be empty"))
+    end
+
+    mission_lower = lowercase(mission)
+
+    if mission_lower in ["chandra", "axaf"]
+        mission_lower = "chandra"
+    end
+
+    calib_func = if haskey(SIMPLE_CALIBRATION_FUNCS, mission_lower)
+        SIMPLE_CALIBRATION_FUNCS[mission_lower]
+    else
+        @warn "Mission $mission not recognized, using identity function"
+        identity
+    end
+
+    energy_alts = if mission_lower in ["chandra", "axaf"]
+        ["ENERGY", "PI", "PHA"]
+    elseif mission_lower == "xte"
+        ["PHA", "PI", "ENERGY"]
+    elseif mission_lower == "nustar"
+        ["PI", "ENERGY", "PHA"]
+    else
+        ["ENERGY", "PI", "PHA"]
+    end
+
+    # Mission-specific GTI extensions
+    gti_exts = if mission_lower == "xmm"
+        ["GTI", "GTI0", "STDGTI"]
+    elseif mission_lower in ["chandra", "axaf"]
+        ["GTI", "GTI0", "GTI1", "GTI2", "GTI3"]
+    else
+        ["GTI", "STDGTI"]
+    end
+
+    MissionSupport{Float64}(mission_lower, instrument, epoch, calib_func, nothing, energy_alts, gti_exts)
+end
+
+"""
+    apply_calibration(mission_support::MissionSupport, pi_channels::AbstractArray) -> Vector{Float64}
+
+Apply calibration function to PI channels.
+
+Converts PI (Pulse Invariant) channel values to energies in keV using
+the mission-specific calibration function stored in the MissionSupport object.
+
+# Arguments
+- `mission_support::MissionSupport`: Mission support object containing calibration function
+- `pi_channels::AbstractArray{T}`: Array of PI channel values
+
+# Returns
+- `Vector{Float64}`: Array of energy values in keV
+
+# Examples
+```julia
+ms = get_mission_support("nustar")
+pi_values = [100, 500, 1000]
+energies = apply_calibration(ms, pi_values)
+```
+"""
+function apply_calibration(mission_support::MissionSupport, pi_channels::AbstractArray{T}) where T
+    if isempty(pi_channels)
+        return similar(pi_channels, Float64)
+    end
+    return mission_support.calibration_func.(pi_channels)
+end
+
+"""
+    patch_mission_info(info::Dict{String,Any}, mission=nothing) -> Dict{String,Any}
+
+Apply mission-specific patches to header information.
+
+This function applies mission-specific modifications to FITS header information
+to handle mission-specific quirks and conventions. It's based on the Python
+implementation in Stingray's mission interpretation module.
+
+# Arguments
+- `info::Dict{String,Any}`: Dictionary containing header information
+- `mission::Union{String,Nothing}=nothing`: Mission name
+
+# Returns
+- `Dict{String,Any}`: Patched header information dictionary
+
+# Examples
+```julia
+info = Dict("gti" => "STDGTI", "ecol" => "PHA")
+patched = patch_mission_info(info, "xmm")  # Adds GTI0 to gti field
+```
+"""
+function patch_mission_info(info::Dict{String,Any}, mission::Union{String,Nothing}=nothing)
+    if isnothing(mission)
+        return info
+    end
+
+    mission_lower = lowercase(mission)
+    patched_info = copy(info)
+
+    # Normalize chandra/axaf
+    if mission_lower in ["chandra", "axaf"]
+        mission_lower = "chandra"
+    end
+
+    if mission_lower == "xmm" && haskey(patched_info, "gti")
+        patched_info["gti"] = string(patched_info["gti"], ",GTI0")
+    elseif mission_lower == "xte" && haskey(patched_info, "ecol")
+        patched_info["ecol"] = "PHA"
+        patched_info["ccol"] = "PCUID"
+    elseif mission_lower == "chandra"
+        # Chandra-specific patches
+        if haskey(patched_info, "DETNAM")
+            patched_info["detector"] = patched_info["DETNAM"]
+        end
+        if haskey(patched_info, "TIMESYS")
+            patched_info["time_system"] = patched_info["TIMESYS"]
+        end
+    end
+
+    return patched_info
+end
+function interpret_fits_data!(f::FITS, mission_support::MissionSupport)
+    # Placeholder for mission-specific interpretation
+    return nothing
+end

test/runtests.jl

@@ -15,4 +15,8 @@ end
 
 @testset "lightcurve" begin
     include("test_lightcurve.jl")
+end
+
+@testset "Synthetic Events Tests" begin
+    include("test_missionSupport.jl")
 end

test/test_events.jl

@@ -669,7 +669,7 @@ let
     @test occursin("4 events", summary_str)
     @test occursin("over 3.0 time units", summary_str)
     @test occursin("energies: 100.0 - 400.0", summary_str)
-    @test occursin("(keV)", summary_str)
+    @test occursin("(ENERGY)", summary_str)
     @test occursin("GTI: 2 intervals (3.0 s exposure)", summary_str)
     @test occursin("1 extra columns", summary_str)