Lightcurve recipes :)

Project.toml

@@ -16,8 +16,10 @@ JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 ProgressBars = "49802e3a-d2f1-5c88-81d8-b72133a6f568"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 ResumableFunctions = "c5292f4c-5179-55e1-98c5-05642aab7184"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
@@ -35,8 +37,10 @@ JuliaFormatter = "1.0.62"
 LinearAlgebra = "1.11.0"
 Logging = "1.11.0"
 NaNMath = "0.3, 1"
+Plots = "1.40.17"
 ProgressBars = "1.4"
 Random = "1.11.0"
+RecipesBase = "1.3.4"
 ResumableFunctions = "0.6"
 StatsBase = "0.33"
 julia = "1.11"

src/Stingray.jl

@@ -6,6 +6,7 @@ using ProgressBars: tqdm as show_progress
 using DocStringExtensions
 using LinearAlgebra
 using Random
+using RecipesBase
 
 include("fourier.jl")
 export positive_fft_bins
@@ -76,4 +77,8 @@ export create_filtered_lightcurve
 export check_gtis
 export split_by_gtis
 
+include("plotting/plots_recipes_lightcurve.jl")
+export create_segments
+include("plotting/plots_recipes_gti.jl")
+export BTIAnalysisPlot
 end

src/events.jl

@@ -290,7 +290,7 @@ gti_info(ev)
 
 # Notes
 - Issues a warning if no GTI information is available
-- Information is logged using the `@info` macro for structured output
+- Information is logged using the `@debug` macro for structured output
 """
 function gti_info(ev::EventList)
     if !has_gti(ev)
@@ -806,10 +806,9 @@ function readevents(
             gti_hdu_candidates, gti_hdu_indices, combine_gtis)
 
         if !isnothing(gti_data)
-            println("Found GTI data: $(size(gti_data, 1)) intervals")
-            println("GTI time range: $(minimum(gti_data)) to $(maximum(gti_data))")
+            @debug "Found GTI data" n_intervals=size(gti_data, 1) time_range=(minimum(gti_data), maximum(gti_data))
         else
-            println("No GTI data found")
+            @debug "No GTI data found"
         end
     end
 

src/lightcurve.jl

@@ -495,6 +495,77 @@ function apply_filters(
 
     return filtered_times, filtered_energies, start_t, stop_t
 end
+"""
+    apply_filters(times, energies, tstart, tstop, energy_filter)
+
+Basic event filtering without GTI consideration.
+
+# Arguments
+- `times::AbstractVector{T}`: Event arrival times
+- `energies::Union{Nothing,AbstractVector{T}}`: Event energies (optional)
+- `tstart::Union{Nothing,Real}`: Start time filter (inclusive)
+- `tstop::Union{Nothing,Real}`: Stop time filter (inclusive)
+- `energy_filter::Union{Nothing,Tuple{Real,Real}}`: Energy range (emin, emax)
+
+# Returns
+`Tuple{Vector, Union{Nothing,Vector}, Real, Real}`: 
+- Filtered times
+- Filtered energies (if provided)
+- Actual start time
+- Actual stop time
+
+# Examples
+```julia
+# Time filtering only
+filtered_times, _, start_t, stop_t = apply_filters(times, nothing, 1000.0, 2000.0, nothing)
+
+# Energy and time filtering
+filtered_times, filtered_energies, start_t, stop_t = apply_filters(
+    times, energies, 1000.0, 2000.0, (0.5, 10.0)
+)
+Notes
+- Energy filter is applied as: emin ≤ energy < emax
+- Time filter is applied as: tstart ≤ time ≤ tstop
+"""
+function apply_filters(
+    times::AbstractVector{T},
+    energies::Union{Nothing,AbstractVector{T}},
+    eventlist::EventList,
+    tstart::Union{Nothing,Real},
+    tstop::Union{Nothing,Real},
+    energy_filter::Union{Nothing,Tuple{Real,Real}},
+    binsize::Real
+) where T
+    mask = trues(length(times))
+
+    # Apply energy filter
+    if !isnothing(energy_filter) && !isnothing(energies)
+        emin, emax = energy_filter
+        mask = mask .& (energies .>= emin) .& (energies .< emax)
+    end
+
+    # Apply time filters
+    if !isnothing(tstart)
+        mask = mask .& (times .>= tstart)
+    end
+    if !isnothing(tstop)
+        mask = mask .& (times .<= tstop)
+    end
+    # If GTI is present, apply GTI mask
+    if has_gti(eventlist)
+        gti_mask, _ = create_gti_mask(times, eventlist.meta.gti, dt=binsize)
+        mask = mask .& gti_mask
+    end
+
+    !any(mask) && throw(ArgumentError("No events remain after applying filters"))
+
+    filtered_times = times[mask]
+    filtered_energies = isnothing(energies) ? nothing : energies[mask]
+    start_t = isnothing(tstart) ? minimum(filtered_times) : tstart
+    stop_t = isnothing(tstop) ? maximum(filtered_times) : tstop
+
+    return filtered_times, filtered_energies, start_t, stop_t
+end
 
 """
     calculate_event_properties(times, energies, dt, bin_centers) -> Vector{EventProperty}
@@ -605,17 +676,24 @@ function extract_metadata(eventlist::EventList, start_time, stop_time, binsize,
         n_filtered_events  # Assume it's already a count
     end
 
+    # Create extra metadata with GTI information[storing purpose]
     extra_metadata = Dict{String,Any}(
         "filtered_nevents" => n_events,
         "total_nevents" => length(eventlist.times),
         "energy_filter" => energy_filter,
-        "binning_method" => "histogram"
+        "binning_method" => "histogram",
+        "gti" => eventlist.meta.gti  #GTI information[this 'eventlist.meta.gti' need to bw rembered since it is the one where u will all gti information:)]
     )
 
     if hasfield(typeof(eventlist.meta), :extra)
         merge!(extra_metadata, eventlist.meta.extra)
     end
 
+    # Add GTI source information if available
+    if !isnothing(eventlist.meta.gti_source)
+        extra_metadata["gti_source"] = eventlist.meta.gti_source
+    end
+
     return LightCurveMetadata(
         telescope, instrument, object, Float64(mjdref),
         (Float64(start_time), Float64(stop_time)), Float64(binsize),
@@ -691,22 +769,20 @@ function create_lightcurve(
     !(err_method in [:poisson, :gaussian]) && throw(ArgumentError(
         "Unsupported error method: $err_method. Use :poisson or :gaussian"
     ))
-
     binsize_t = convert(T, binsize)
 
-    # Apply filters to get filtered times and energies
+
     filtered_times, filtered_energies, start_t, stop_t = apply_filters(
         eventlist.times,
         eventlist.energies,
+        eventlist,
         tstart,
         tstop,
-        energy_filter
+        energy_filter,
+        binsize_t
     )
 
-    # Check if we have any events left after filtering
-    if isempty(filtered_times)
-        throw(ArgumentError("No events remain after filtering"))
-    end
+    isempty(filtered_times) && throw(ArgumentError("No events remain after filtering"))
 
     # Determine time range
     start_time = minimum(filtered_times)
@@ -716,17 +792,17 @@ function create_lightcurve(
     dt, bin_centers = create_time_bins(start_time, stop_time, binsize_t)
     counts = bin_events(filtered_times, dt)
 
-    @info "Created light curve: $(length(bin_centers)) bins, bin size = $(binsize_t) s"
+    @debug "Created light curve: $(length(bin_centers)) bins, bin size = $(binsize_t) s"
 
     # Calculate exposure and properties
     exposure = fill(binsize_t, length(bin_centers))
     properties = calculate_event_properties(filtered_times, filtered_energies, dt, bin_centers)
 
-    # Extract metadata - use filtered time range if time filtering was applied
+    # Extract metadata with GTI information
     actual_start = !isnothing(tstart) ? T(tstart) : start_time
     actual_stop = !isnothing(tstop) ? T(tstop) : stop_time
     metadata = extract_metadata(eventlist, actual_start, actual_stop, binsize_t, 
-                               length(filtered_times), energy_filter)
+                              length(filtered_times), energy_filter)
 
     # Create light curve (errors will be calculated when needed)
     lc = LightCurve{T}(
@@ -737,6 +813,11 @@ function create_lightcurve(
     # Calculate initial errors
     calculate_errors!(lc)
 
+    # Add debug info about GTI
+    if has_gti(eventlist)
+        @debug "GTI information preserved" n_intervals=size(eventlist.meta.gti, 1) time_range=extrema(eventlist.meta.gti)
+    end
+
     return lc
 end
 """