Exposure Time Calculations

Overview

Exposure time calculations determine how long a gamma-ray observatory needs to observe a transient source to achieve a specified detection significance. This is one of the core functionalities of sensipy and is particularly useful for rapidly-evolving time-varying sources.

Given:

• A time-evolving source spectrum
• Observatory sensitivity
• An observation start time (delay from event onset)
• A target significance level $\sigma$

sensipy calculates the required observation time to reach that significance $\sigma$, or reports that the source is not detectable.

The observe Method

The primary method for exposure calculations is Source.observe():

from sensipy.source import Source
from sensipy.sensitivity import Sensitivity
from sensipy.ctaoirf import IRFHouse
from sensipy.util import get_data_path
import astropy.units as u

# Load IRFs
house = IRFHouse(base_directory="./IRFs/CTAO")
irf = house.get_irf(
    site="south",
    configuration="alpha",
    zenith=20,
    duration=1800,
    azimuth="average",
    version="prod5-v0.1",
)

# Load source
mock_data_path = get_data_path("mock_data/GRB_42_mock.csv")

# Define energy range (typical for CTA)
min_energy = 20 * u.GeV
max_energy = 10 * u.TeV

source = Source(
    filepath=str(mock_data_path),
    min_energy=min_energy,
    max_energy=max_energy,
    ebl="franceschini"
)

# Create sensitivity calculator
sens = Sensitivity(
    irf=irf,
    observatory="ctao_south",
    min_energy=min_energy,
    max_energy=max_energy,
    radius=3.0 * u.deg,
)
sens.get_sensitivity_curve(source=source)

# Simulate observation starting 30 minutes after event
result = source.observe(
    sensitivity=sens,
    start_time=30 * u.min,
    min_energy=min_energy,
    max_energy=max_energy,
)

print(f"Required observation time: {result['obs_time']}")
print(f"Detection achieved: {result['seen']}")
print(f"Start time: {result['start_time']}")
print(f"End time: {result['end_time']}")

Key Parameters

Required Parameters

Parameter	Description	Example
sensitivity	Sensitivity object with computed curves	sens
start_time	Delay from event onset to start observing	30 * u.min
min_energy	Minimum energy for detection	20 * u.GeV
max_energy	Maximum energy for detection	10 * u.TeV

Optional Parameters

Parameter	Description	Default
target_precision	Precision for binary search	1 * u.s
max_time	Maximum allowed observation time	12 * u.h
sensitivity_mode	"sensitivity" or "photon_flux"	"sensitivity"
target_significance	Target sigma for detection	5.0
n_time_steps	Number of time steps for integration	10

Return Values

The observe() method returns a dictionary with the following keys:

{
    'obs_time': <Quantity>,      # Required observation time
    'seen': bool,                 # True if detection is possible
    'start_time': <Quantity>,     # Observation start time
    'end_time': <Quantity>,       # Observation end time (start + obs_time)
    'min_energy': <Quantity>,     # Energy range minimum
    'max_energy': <Quantity>,     # Energy range maximum
    'ebl_model': str,             # EBL model used (or None)
    'error_message': str,         # Error message if not detected

    # Source metadata (if available)
    'long': <Quantity>,           # RA/longitude
    'lat': <Quantity>,            # Dec/latitude
    'dist': <Quantity>,           # Distance
    'id': int,                    # Event ID
}

Usage Examples

Basic Detection Simulation

Simple Case

import astropy.units as u
from sensipy.source import Source
from sensipy.sensitivity import Sensitivity
from sensipy.ctaoirf import IRFHouse

# Setup
house = IRFHouse(base_directory="./IRFs/CTAO")
irf = house.get_irf(
    site="south",
    configuration="alpha",
    zenith=20,
    duration=1800,
    azimuth="average",
    version="prod5-v0.1",
)

from sensipy.util import get_data_path
mock_data_path = get_data_path("mock_data/GRB_42_mock.csv")

# Define energy range (typical for CTA)
min_energy = 20 * u.GeV
max_energy = 10 * u.TeV

source = Source(
    filepath=str(mock_data_path),
    min_energy=min_energy,
    max_energy=max_energy,
    ebl="franceschini"
)

sens = Sensitivity(
    irf=irf,
    observatory="ctao_south",
    min_energy=min_energy,
    max_energy=max_energy,
    radius=3.0 * u.deg,
)
sens.get_sensitivity_curve(source=source)

# Observe starting 30 minutes after event
result = source.observe(
    sensitivity=sens,
    start_time=30 * u.min,
    min_energy=min_energy,
    max_energy=max_energy,
)

if result['seen']:
    print(f"✓ Detectable in {result['obs_time']}")
else:
    print(f"✗ Not detectable: {result['error_message']}")

Multiple Delays

import numpy as np
import astropy.units as u
import matplotlib.pyplot as plt

# Test detection at different delay times
delays = np.logspace(1, 4, 20) * u.s  # 10s to 10,000s
obs_times = []

for delay in delays:
    result = source.observe(
        sensitivity=sens,
        start_time=delay,
        min_energy=min_energy,
        max_energy=max_energy,
    )

    if result['seen']:
        obs_times.append(result['obs_time'].to(u.s).value)
    else:
        obs_times.append(np.nan)  # Not detectable

# Plot results
plt.figure(figsize=(10, 6))
plt.loglog(delays, obs_times, 'o-')
plt.xlabel("Delay from Event [s]")
plt.ylabel("Required Observation Time [s]")
plt.title("Detectability vs. Observation Delay")
plt.grid(True, alpha=0.3)
plt.show()

With EBL Comparison

import astropy.units as u

delay = 30 * u.min

from sensipy.util import get_data_path
mock_data_path = get_data_path("mock_data/GRB_42_mock.csv")

# Define energy range
min_energy = 20 * u.GeV
max_energy = 10 * u.TeV

for ebl_model in [None, "franceschini", "dominguez"]:
    source = Source(
        filepath=str(mock_data_path),
        min_energy=min_energy,
        max_energy=max_energy,
        ebl=ebl_model
    )

    sens.get_sensitivity_curve(source=source)
    result = source.observe(
        sensitivity=sens,
        start_time=delay,
        min_energy=min_energy,
        max_energy=max_energy,
    )

    ebl_label = ebl_model if ebl_model else "No EBL"
    if result['seen']:
        print(f"{ebl_label:20s}: {result['obs_time']}")
    else:
        print(f"{ebl_label:20s}: Not detectable")

Understanding the Algorithm

The exposure calculation uses a binary search algorithm to find the minimum observation time:

1. Set search bounds

   • Lower bound: Short observation time (e.g., 10s)
   • Upper bound: Maximum allowed time (e.g., 12 hours)

2. Compare source flux to sensitivity

   • Calculate the source flux at the given start time
   • Compare to the sensitivity curve at different exposure times

3. Binary search

   • Test the midpoint exposure time
   • If source is detectable → try shorter time
   • If not detectable → try longer time
   • Repeat until precision target is reached

4. Return result

   • If detectable within max_time → return observation time
   • Otherwise → return seen=False with error message

Sensitivity Modes

The algorithm can use two different metrics:

1. Integral Sensitivity Mode (default: sensitivity_mode="sensitivity")

• Integrates source flux and sensitivity over energy
• Compares integrated energy flux
• Used for most applications

2. Photon Flux Mode (sensitivity_mode="photon_flux")

• Uses photon flux instead of energy flux

# Use photon flux mode
# Assume source and sens are already set up with min_energy = 20 * u.GeV
result = source.observe(
    sensitivity=sens,
    start_time=30 * u.min,
    min_energy=min_energy,
    max_energy=max_energy,
    sensitivity_mode="photon_flux",  # Use photon flux
)

Advanced Parameters

Precision Control

The target_precision parameter controls how accurately the observation time is determined:

# High precision (slower)
result_precise = source.observe(
    sensitivity=sens,
    start_time=30 * u.min,
    target_precision=0.1 * u.s,  # Precise to 0.1 seconds
    # ... other params
)

# Lower precision (faster)
result_fast = source.observe(
    sensitivity=sens,
    start_time=30 * u.min,
    target_precision=10 * u.s,  # Precise to 10 seconds
    # ... other params
)

Note

For most applications, 1 * u.s precision (default) is sufficient and provides a good balance between accuracy and speed.

Maximum Time Limit

The max_time parameter sets the longest observation time to consider:

# Only consider observations up to 2 hours
result = source.observe(
    sensitivity=sens,
    start_time=30 * u.min,
    max_time=2 * u.h,  # Max 2 hours
    # ... other params
)

if not result['seen']:
    print(f"Not detectable within {2 * u.h}")

Time Step Integration

The n_time_steps parameter controls how finely the time-integrated spectrum is sampled:

result = source.observe(
    sensitivity=sens,
    start_time=30 * u.min,
    n_time_steps=20,  # More steps = better accuracy
    # ... other params
)

Batch Processing Multiple Events

For large-scale studies (e.g. catalogs), you can batch process many events:

import pandas as pd
from pathlib import Path

# Catalog of event files
source_files = Path("./event_catalog/").glob("event_*.csv")
results_list = []

# Define energy range
min_energy = 20 * u.GeV
max_energy = 10 * u.TeV

for source_file in source_files:
    source = Source(
        filepath=str(source_file),
        min_energy=min_energy,
        max_energy=max_energy,
        ebl="franceschini"
    )

    sens.get_sensitivity_curve(source=source)
    result = source.observe(
        sensitivity=sens,
        start_time=30 * u.min,
        min_energy=min_energy,
        max_energy=max_energy,
    )

    results_list.append(result)

# Convert to DataFrame for analysis
df = pd.DataFrame(results_list)
print(f"Detection fraction: {df['seen'].sum() / len(df):.2%}")
print(f"Median obs time: {df[df['seen']]['obs_time'].median()}")

With a parallelization package like ray or dask, you can speed up the previous loop significantly.

Problem: All sources report as not detectable

Solution:

• Check that min_energy and max_energy match between source, sensitivity, and observe() call
• Verify the spectral model has sufficient flux
• Try a longer max_time
• Verify EBL absorption isn’t too strong (try without EBL first)

---

Problem: Observation times are unrealistically short/long

Solution:

• Check the units of your spectral model (should be cm⁻² s⁻¹ GeV⁻¹)
• Verify the IRF loading is correct
• Ensure the sensitivity curve was calculated for the correct source

---

Problem: Inconsistent results between runs

Solution:

• This may indicate numerical instability - try adjusting target_precision or n_time_steps
• Check for edge cases (very weak or very strong sources)