Finding emission-line candidates in the Vela X-1 RGS spectrum¶

This worked example uses RGS1, first order, observation 0841890201, with a search interval of 0.62–0.92 keV. It illustrates loading a spectrum, exploring rebinning, inspecting an empirical baseline, detecting candidates, evaluating them with synthetic spectra, and fitting and identifying selected features.

In [1]:

Copied!





# Import the public BLiSS API and the utilities used throughout this notebook.
from bliss import *
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
from time import perf_counter

from bliss.line_search.blind_line_search import (
    BlindLineSearchConfig,)
# Import the public BLiSS API and the utilities used throughout this notebook.
from bliss import *
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
from time import perf_counter

from bliss.line_search.blind_line_search import (
    BlindLineSearchConfig,)

1.1 Define the configuration¶

BlindLineSearchConfig holds the search settings. Commented-out arguments use the library defaults; they do not inherit values from earlier cells. Here the active overrides include a 0.1 keV energy margin, a maximum Gaussian sigma of 0.005 keV (5 eV), and Poisson simulations.

Execution dependency: en1 and en2 are first assigned later in the saved notebook. For a fresh session, define en1, en2 = 0.62, 0.92 before this cell.

Creating pipeline does not execute a search. The later calls to find_candidate_lines construct their own pipeline. Only the evaluated search explicitly receives this config; keyword arguments in that call can override it.

In [2]:

Copied!





en1 = 0.5
en2 = 2.2
# Define en1 and en2 before this cell in a fresh kernel.
# Commented arguments use defaults; this pipeline object is not run below.
config = BlindLineSearchConfig(
    # Search interval (keV)
    en1=en1,                           # Default: 0.2
    en2=en2,                           # Default: 10.0
    energy_pad=0.1,                    # Extra margin on each side, in keV

    # Rebinning
    #rebin_method=rebin_method,         # Default: "none"
    #rebin_scale=rebin_scale,           # Default: None
    #rebin_min_bins=1,                  # Minimum input bins per group for "snr"

    # Baseline
    #baseline_window=window_fn,         # Default: 0.4 keV
    #max_range_fraction=0.2,            # Maximum window / input energy range
    #min_points=3,                     # Minimum points for the local median

    # Line fitting
    max_sigma_line=0.005,                # Maximum Gaussian sigma, in keV
    #final_fit_maxfev=100000,           # Maximum global-fit function evaluations

    # Simulations
    #num_synthetic_simulations=100,     # Default: 10; 0 disables simulations
    #synthetic_seed=None,              # Use an integer for reproducibility
    noise_model="poisson",             # "gaussian" for validation/comparison

    # Evaluation and flags
    #min_area_snr=None,                # No preliminary area-S/N cut for the GMM
    #response_feature_threshold=5.0,   # Threshold for flagging ARF structure
)

pipeline = BlindLineSearchPipeline(config=config)
en1 = 0.5
en2 = 2.2
# Define en1 and en2 before this cell in a fresh kernel.
# Commented arguments use defaults; this pipeline object is not run below.
config = BlindLineSearchConfig(
    # Search interval (keV)
    en1=en1,                           # Default: 0.2
    en2=en2,                           # Default: 10.0
    energy_pad=0.1,                    # Extra margin on each side, in keV

    # Rebinning
    #rebin_method=rebin_method,         # Default: "none"
    #rebin_scale=rebin_scale,           # Default: None
    #rebin_min_bins=1,                  # Minimum input bins per group for "snr"

    # Baseline
    #baseline_window=window_fn,         # Default: 0.4 keV
    #max_range_fraction=0.2,            # Maximum window / input energy range
    #min_points=3,                     # Minimum points for the local median

    # Line fitting
    max_sigma_line=0.005,                # Maximum Gaussian sigma, in keV
    #final_fit_maxfev=100000,           # Maximum global-fit function evaluations

    # Simulations
    #num_synthetic_simulations=100,     # Default: 10; 0 disables simulations
    #synthetic_seed=None,              # Use an integer for reproducibility
    noise_model="poisson",             # "gaussian" for validation/comparison

    # Evaluation and flags
    #min_area_snr=None,                # No preliminary area-S/N cut for the GMM
    #response_feature_threshold=5.0,   # Threshold for flagging ARF structure
)

pipeline = BlindLineSearchPipeline(config=config)

2. Locate and load the RGS files and load the spectrum¶

File	Role
`vela_RGS1_order1_grp.fits`	Source spectrum with grouping flags previously written by SAS
`P0841890201R1S004RSPMAT1001.FIT`	Corresponding RGS response
`P0841890201R1S004BGSPEC1001.FIT`	Background spectrum

In [34]:

Copied!





# Load a background-subtracted spectrum; SAS GROUPING flags are not applied by this loader.
rgs_path = Path("./0841890201/odf/")  # Directory containing the FITS files

pha = rgs_path / "vela_RGS1_order1_grp.fits"
rmf = rgs_path / "P0841890201R1S004RSPMAT1001.FIT"
bkg = rgs_path / "P0841890201R1S004BGSPEC1001.FIT"

spectrum = prepare_spectrum(
    pha_path=pha,
    rmf_path=rmf,
    background_path=bkg,
)
# Load a background-subtracted spectrum; SAS GROUPING flags are not applied by this loader.
rgs_path = Path("./0841890201/odf/")  # Directory containing the FITS files

pha = rgs_path / "vela_RGS1_order1_grp.fits"
rmf = rgs_path / "P0841890201R1S004RSPMAT1001.FIT"
bkg = rgs_path / "P0841890201R1S004BGSPEC1001.FIT"

spectrum = prepare_spectrum(
    pha_path=pha,
    rmf_path=rmf,
    background_path=bkg,
)

3. Inspect the native spectrum and compare rebinning choices¶

The plotting helper restricts the plotted arrays to the requested energy interval. The printed median bin width is computed over the entire loaded spectrum, not just the displayed interval.

3.1 Native channels¶

rebin_method="none" preserves the retained input channels. Unlike Section 2, this call does not pass a background file, so native contains the source-region spectrum without background subtraction. It is a separate object from spectrum.

In [4]:

Copied!





def plot_spectrum(spectra, labels, xlim, title=None):
    plt.figure(figsize=(9, 4.5))
    for i, (s, lab) in enumerate(zip(spectra, labels)):
        m = (s.energy >= xlim[0]) & (s.energy <= xlim[1])
        plt.errorbar(s.energy[m], s.values[m], yerr=s.uncertainties[m],
                     fmt=".", ms=3 if i == 0 else 4, alpha=0.3 if i == 0 else 0.9, label=lab)
    plt.xlabel("Energy (keV)")
    plt.ylabel("counts s$^{-1}$ keV$^{-1}$")
    if title:
        plt.title(title)
    plt.legend()
    plt.tight_layout()
    plt.show()


native = prepare_spectrum(pha, rmf, rebin_method="none")
print(f"{len(native.energy)} channels, median width {1e3*np.median(native.bin_width):.2f} eV, "
      f"exposure {native.native.exposure:.0f} s")

plot_spectrum([native], ["RGS Vela X-1"], xlim=(0.62, 0.92))
def plot_spectrum(spectra, labels, xlim, title=None):
    plt.figure(figsize=(9, 4.5))
    for i, (s, lab) in enumerate(zip(spectra, labels)):
        m = (s.energy >= xlim[0]) & (s.energy <= xlim[1])
        plt.errorbar(s.energy[m], s.values[m], yerr=s.uncertainties[m],
                     fmt=".", ms=3 if i == 0 else 4, alpha=0.3 if i == 0 else 0.9, label=lab)
    plt.xlabel("Energy (keV)")
    plt.ylabel("counts s$^{-1}$ keV$^{-1}$")
    if title:
        plt.title(title)
    plt.legend()
    plt.tight_layout()
    plt.show()


native = prepare_spectrum(pha, rmf, rebin_method="none")
print(f"{len(native.energy)} channels, median width {1e3*np.median(native.bin_width):.2f} eV, "
      f"exposure {native.native.exposure:.0f} s")

plot_spectrum([native], ["RGS Vela X-1"], xlim=(0.62, 0.92))

2303 channels, median width 0.26 eV, exposure 113470 s

No description has been provided for this image

3.2 Group by signal-to-noise ratio¶

rebin_method="snr", rebin_scale=3.0 accumulates adjacent input bins until their combined signal-to-noise ratio reaches 3. This can produce wider bins in weak parts of the spectrum; a trailing group can fall below the requested threshold.

In [5]:

Copied!





# S/N-based grouping of the source-region spectrum (no background passed here).
spectrum_snr = prepare_spectrum(
    pha,
    rmf,
    rebin_method="snr",
    rebin_scale=3.0,
)

print(
    f"{len(spectrum_snr.energy)} bins, "
    f"median width {1e3 * np.median(spectrum_snr.bin_width):.2f} eV, "
    f"exposure {spectrum_snr.native.exposure:.0f} s"
)

plot_spectrum(
    [spectrum_snr],
    ["RGS Vela X-1 — S/N ≥ 3"],
    xlim=(0.62, 0.92),
)
# S/N-based grouping of the source-region spectrum (no background passed here).
spectrum_snr = prepare_spectrum(
    pha,
    rmf,
    rebin_method="snr",
    rebin_scale=3.0,
)

print(
    f"{len(spectrum_snr.energy)} bins, "
    f"median width {1e3 * np.median(spectrum_snr.bin_width):.2f} eV, "
    f"exposure {spectrum_snr.native.exposure:.0f} s"
)

plot_spectrum(
    [spectrum_snr],
    ["RGS Vela X-1 — S/N ≥ 3"],
    xlim=(0.62, 0.92),
)

981 bins, median width 1.38 eV, exposure 113470 s

3.3 Group a fixed number of channels¶

rebin_method="bins", rebin_scale=3 combines groups of three retained input bins, sums their values, and propagates errors in quadrature. The implementation merges any incomplete trailing group into the preceding group. This is a channel-count criterion, not a constant energy width.

Here background_path=bkg is supplied, so spectrum_bins is background-subtracted. To compare grouping methods alone, all examples would need the same treatment of the background.

In [6]:

Copied!





# Combine three input bins per group after background subtraction.
spectrum_bins = prepare_spectrum(
    pha_path=pha,
    rmf_path=rmf,
    background_path=bkg,
    rebin_method="bins",
    rebin_scale=3,
)

print(
    f"{len(spectrum_bins.energy)} bins, "
    f"median width {1e3 * np.median(spectrum_bins.bin_width):.2f} eV, "
    f"exposure {spectrum_bins.native.exposure:.0f} s"
)

plot_spectrum(
    [spectrum_bins],
    ["RGS Vela X-1 — agrupado de 3 en 3"],
    xlim=(0.62, 0.92),
)
# Combine three input bins per group after background subtraction.
spectrum_bins = prepare_spectrum(
    pha_path=pha,
    rmf_path=rmf,
    background_path=bkg,
    rebin_method="bins",
    rebin_scale=3,
)

print(
    f"{len(spectrum_bins.energy)} bins, "
    f"median width {1e3 * np.median(spectrum_bins.bin_width):.2f} eV, "
    f"exposure {spectrum_bins.native.exposure:.0f} s"
)

plot_spectrum(
    [spectrum_bins],
    ["RGS Vela X-1 — agrupado de 3 en 3"],
    xlim=(0.62, 0.92),
)

891 bins, median width 0.71 eV, exposure 113470 s

3.4 Repeated fixed-channel grouping cell¶

For a separate fixed-energy-width experiment, BLiSS accepts rebin_method="resolution", rebin_scale=0.002 (2 eV). That method uses fixed energy intervals; it does not automatically follow the RGS instrumental resolution.

In [35]:

Copied!





# This repeats the preceding fixed-channel example; it is not resolution-based grouping.
spectrum_bins = prepare_spectrum(
    pha_path=pha,
    rmf_path=rmf,
    background_path=bkg,
    rebin_method="bins",
    rebin_scale=3,
)

print(
    f"{len(spectrum_bins.energy)} bins, "
    f"median width {1e3 * np.median(spectrum_bins.bin_width):.2f} eV, "
    f"exposure {spectrum_bins.native.exposure:.0f} s"
)

plot_spectrum(
    [spectrum_bins],
    ["RGS Vela X-1 — agrupado de 3 en 3"],
    xlim=(0.62, 0.92),
)
# This repeats the preceding fixed-channel example; it is not resolution-based grouping.
spectrum_bins = prepare_spectrum(
    pha_path=pha,
    rmf_path=rmf,
    background_path=bkg,
    rebin_method="bins",
    rebin_scale=3,
)

print(
    f"{len(spectrum_bins.energy)} bins, "
    f"median width {1e3 * np.median(spectrum_bins.bin_width):.2f} eV, "
    f"exposure {spectrum_bins.native.exposure:.0f} s"
)

plot_spectrum(
    [spectrum_bins],
    ["RGS Vela X-1 — agrupado de 3 en 3"],
    xlim=(0.62, 0.92),
)

891 bins, median width 0.71 eV, exposure 113470 s

4. Inspect the empirical baseline¶

BLiSS estimates the baseline using a running median in energy. This provides an empirical description of the observed spectral shape rather than a physical continuum model. The method exploits the difference in scale between the smoothly varying continuum and narrow emission-line features. A window that is too narrow may absorb part of the line structure into the baseline, whereas a window that is too broad may fail to follow relatively rapid changes in the underlying continuum. If different spectral regions require different smoothing scales, an array of window sizes can also be provided.

In [8]:

Copied!





baseline_window = 0.4 #Default

baseline_bins, info_bins = base_calculator(
    spectrum_bins.energy,
    spectrum_bins.values,
    baseline_window=baseline_window,
    return_info=True)


fig, ax = plt.subplots(figsize=(6, 4))

ax.errorbar(
    spectrum_bins.energy,
    spectrum_bins.values,
    yerr = spectrum_bins.uncertainties,
    fmt=".",
    color="black",
    linewidth=0.8,
   
    label="RGS Vela X-1", 
)

ax.plot(
    spectrum_bins.energy,
    baseline_bins,
    color="tab:orange",
    linewidth=2,
    label="Baseline",
)

ax.set(
    xlim=(0.62, 0.92),
    xlabel="Energía (keV)",
    ylabel=r"Cuentas s$^{-1}$ keV$^{-1}$",
)

plt.ylim(-0.05,0.25)
ax.legend()
fig.tight_layout()
plt.show()
baseline_window = 0.4 #Default

baseline_bins, info_bins = base_calculator(
    spectrum_bins.energy,
    spectrum_bins.values,
    baseline_window=baseline_window,
    return_info=True)


fig, ax = plt.subplots(figsize=(6, 4))

ax.errorbar(
    spectrum_bins.energy,
    spectrum_bins.values,
    yerr = spectrum_bins.uncertainties,
    fmt=".",
    color="black",
    linewidth=0.8,
   
    label="RGS Vela X-1", 
)

ax.plot(
    spectrum_bins.energy,
    baseline_bins,
    color="tab:orange",
    linewidth=2,
    label="Baseline",
)

ax.set(
    xlim=(0.62, 0.92),
    xlabel="Energía (keV)",
    ylabel=r"Cuentas s$^{-1}$ keV$^{-1}$",
)

plt.ylim(-0.05,0.25)
ax.legend()
fig.tight_layout()
plt.show()

4.2 Compare a smaller requested window: 0.01 keV¶

This cell recomputes baseline_bins with a 10 eV window. It overwrites the previous baseline but does not update config.baseline_window.

In [18]:

Copied!





# Compare a smaller baseline window on the same spectrum.
# TODO: use info_bins in the print and baseline_bins in the plotted curve.
baseline_window = 0.01


baseline_bins, info_bins = base_calculator(
    spectrum_bins.energy,
    spectrum_bins.values,
    baseline_window=baseline_window,
    return_info=True,
)


fig, ax = plt.subplots(figsize=(6, 4))

ax.errorbar(
    spectrum_bins.energy,
    spectrum_bins.values,
    yerr = spectrum_bins.uncertainties,
        fmt=".",
    color="black",
    linewidth=0.8,
    label="RGS Vela X-1", 
)

ax.plot(
    spectrum_bins.energy,
    baseline_bins,
    color="tab:orange",
    linewidth=2,
    label="Baseline",
)

ax.set(
    xlim=(0.62, 0.92),
    xlabel="Energía (keV)",
    ylabel=r"Cuentas s$^{-1}$ keV$^{-1}$",
)

plt.ylim(-0.05,0.25)
ax.legend()
fig.tight_layout()
plt.show()
# Compare a smaller baseline window on the same spectrum.
# TODO: use info_bins in the print and baseline_bins in the plotted curve.
baseline_window = 0.01


baseline_bins, info_bins = base_calculator(
    spectrum_bins.energy,
    spectrum_bins.values,
    baseline_window=baseline_window,
    return_info=True,
)


fig, ax = plt.subplots(figsize=(6, 4))

ax.errorbar(
    spectrum_bins.energy,
    spectrum_bins.values,
    yerr = spectrum_bins.uncertainties,
        fmt=".",
    color="black",
    linewidth=0.8,
    label="RGS Vela X-1", 
)

ax.plot(
    spectrum_bins.energy,
    baseline_bins,
    color="tab:orange",
    linewidth=2,
    label="Baseline",
)

ax.set(
    xlim=(0.62, 0.92),
    xlabel="Energía (keV)",
    ylabel=r"Cuentas s$^{-1}$ keV$^{-1}$",
)

plt.ylim(-0.05,0.25)
ax.legend()
fig.tight_layout()
plt.show()

5. Run a quick candidate search without simulations¶

Setting num_synthetic_simulations=0 detects and locally fits candidates without assigning simulation-based significance. This is useful for an initial inspection; the resulting candidates are not confirmed detections.

This call uses two-bin grouping, omits the background, and does not receive config. It therefore uses the default baseline and line-width settings rather than the custom configuration above. The timer measures the entire function call, including loading and preparation.

In [25]:

Copied!





# Quick search: no simulations, no background, and no custom config.
# The two-bin grouping differs from the evaluated search below.
en1, en2 = 0.62, 0.92

inicio = perf_counter()

candidate_lines_not_evaluated = find_candidate_lines(
    pha=pha,
    rmf=rmf,
    rebin_method="bins",
    rebin_scale=2,
    en1=en1,
    en2=en2,
    num_synthetic_simulations=0,
    synthetic_seed=1,
)

tiempo = perf_counter() - inicio

print(f"Tiempo total: {tiempo:.2f} s ({tiempo / 60:.2f} min)")
# Quick search: no simulations, no background, and no custom config.
# The two-bin grouping differs from the evaluated search below.
en1, en2 = 0.62, 0.92

inicio = perf_counter()

candidate_lines_not_evaluated = find_candidate_lines(
    pha=pha,
    rmf=rmf,
    rebin_method="bins",
    rebin_scale=2,
    en1=en1,
    en2=en2,
    num_synthetic_simulations=0,
    synthetic_seed=1,
)

tiempo = perf_counter() - inicio

print(f"Tiempo total: {tiempo:.2f} s ({tiempo / 60:.2f} min)")

Tiempo total: 4.92 s (0.08 min)

5.1 Overlay candidate energies¶

The vertical markers show the centers returned by the quick search.

In [26]:

Copied!





fig, ax = plt.subplots(figsize=(6, 4))

ax.errorbar(
    spectrum_bins.energy,
    spectrum_bins.values,
    yerr = spectrum_bins.uncertainties,
        fmt=".",
    color="black",
    linewidth=0.8,
    label="RGS Vela X-1", 
)

ax.set(
    xlim=(0.62, 0.92),
    xlabel="Energía (keV)",
    ylabel=r"Cuentas s$^{-1}$ keV$^{-1}$",
)

for i in range(len(candidate_lines_not_evaluated )):
    plt.vlines(x=candidate_lines_not_evaluated .center[i], ymin=-1, ymax=2)



plt.ylim(-0.05,0.35)
ax.legend()
fig.tight_layout()
plt.show()
fig, ax = plt.subplots(figsize=(6, 4))

ax.errorbar(
    spectrum_bins.energy,
    spectrum_bins.values,
    yerr = spectrum_bins.uncertainties,
        fmt=".",
    color="black",
    linewidth=0.8,
    label="RGS Vela X-1", 
)

ax.set(
    xlim=(0.62, 0.92),
    xlabel="Energía (keV)",
    ylabel=r"Cuentas s$^{-1}$ keV$^{-1}$",
)

for i in range(len(candidate_lines_not_evaluated )):
    plt.vlines(x=candidate_lines_not_evaluated .center[i], ymin=-1, ymax=2)



plt.ylim(-0.05,0.35)
ax.legend()
fig.tight_layout()
plt.show()

6. Evaluate candidates with 100 synthetic spectra¶

The evaluated search uses three-bin grouping and a fixed random seed. Poisson null spectra provide a reference for the BLiSS score and the five p_global_* statistics.

With 100 simulations, the smallest attainable Monte Carlo p-value is 1 / (100 + 1) ≈ 0.0099. Each global p-value compares one statistic with the maximum found across the search interval in each null realization.

In [28]:

Copied!





# Evaluated search: these keyword arguments override the supplied config.
# No background is passed; the wrapper also defaults energy_pad to zero.
en1, en2 = 0.62, 0.92

inicio = perf_counter()

candidate_lines_evaluated = find_candidate_lines(
    pha=pha,
    rmf=rmf,
    rebin_method="bins",
    rebin_scale=2,
    en1=en1,
    en2=en2,
    num_synthetic_simulations=100,
    synthetic_seed=1,
    config=config
)

tiempo = perf_counter() - inicio

print(f"Tiempo total: {tiempo:.2f} s ({tiempo / 60:.2f} min)")
plot_bliss_score(candidate_lines_evaluated)
# Evaluated search: these keyword arguments override the supplied config.
# No background is passed; the wrapper also defaults energy_pad to zero.
en1, en2 = 0.62, 0.92

inicio = perf_counter()

candidate_lines_evaluated = find_candidate_lines(
    pha=pha,
    rmf=rmf,
    rebin_method="bins",
    rebin_scale=2,
    en1=en1,
    en2=en2,
    num_synthetic_simulations=100,
    synthetic_seed=1,
    config=config
)

tiempo = perf_counter() - inicio

print(f"Tiempo total: {tiempo:.2f} s ({tiempo / 60:.2f} min)")
plot_bliss_score(candidate_lines_evaluated)

Tiempo total: 24.30 s (0.41 min)

6.1 Inspect the evaluated candidate catalogue¶

Review fitted centers and widths, their uncertainties, areas, equivalent widths, BLiSS scores and global p-values.

In [29]:

Copied!

candidate_lines_evaluated
candidate_lines_evaluated

Out[29]:

	center	ecenter	sigma	esigma	amplitude	eamplitude	area	earea	ew	relative_power	snr_peak	snr_area	bliss_score	response_feature	p_global_area	p_global_snr_area	p_global_snr_peak	p_global_ew	p_global_relative_power
0	0.654858	0.000389	0.001489	0.000346	0.168458	0.069991	0.000629	0.000184	9.731621	0.495620	4.514946	3.411284	0.97	False	0.039604	0.019802	0.158416	0.039604	0.009901
1	0.666276	0.000206	0.000425	0.000140	0.120340	0.048330	0.000128	0.000049	1.996888	0.484043	3.566491	2.590678	0.66	False	0.792079	0.079208	0.257426	0.712871	0.009901
2	0.674441	0.000719	0.000938	0.000698	0.039789	0.027878	0.000094	0.000062	1.442680	0.194666	1.274415	1.518102	0.76	False	0.950495	0.851485	1.000000	0.950495	1.000000
3	0.728888	0.000659	0.000481	0.000392	0.038487	0.032518	0.000046	0.000045	0.722312	0.230611	1.500361	1.038997	0.58	False	1.000000	1.000000	0.990099	1.000000	0.940594
4	0.740764	0.000266	0.001144	0.000230	0.154728	0.034270	0.000444	0.000088	6.911610	0.512523	4.297775	5.062115	0.97	False	0.059406	0.009901	0.158416	0.059406	0.009901
5	0.745201	0.000620	0.001207	0.000591	0.055658	0.025987	0.000168	0.000074	2.619034	0.254779	1.856703	2.277022	0.76	False	0.475248	0.227723	0.960396	0.435644	0.871287
6	0.775804	0.001912	0.000643	0.000709	0.157047	0.300922	0.000253	0.000708	3.968915	0.499244	5.302232	0.357303	0.88	False	0.227723	1.000000	0.108911	0.217822	0.009901
7	0.819364	0.000409	0.000602	0.000258	0.077397	0.035748	0.000117	0.000057	1.670122	0.353894	2.995098	2.059655	0.66	False	0.851485	0.396040	0.425743	0.851485	0.178218
8	0.827235	0.000891	0.000875	0.000786	0.041519	0.029876	0.000091	0.000079	1.284386	0.223651	1.613572	1.150588	0.66	False	0.960396	1.000000	0.990099	0.960396	0.970297
9	0.830111	0.001452	0.000921	0.001541	0.024969	0.027354	0.000058	0.000085	0.804222	0.097591	0.970380	0.678712	0.58	False	1.000000	1.000000	1.000000	1.000000	1.000000
10	0.853213	0.000824	0.001653	0.000793	0.050685	0.022812	0.000210	0.000090	2.821594	0.171288	1.905282	2.343529	0.76	False	0.297030	0.178218	0.940594	0.356436	1.000000
11	0.881865	0.000448	0.001585	0.000421	0.149449	0.042294	0.000594	0.000192	8.248092	0.463636	4.615811	3.099036	0.97	False	0.039604	0.019802	0.148515	0.049505	0.009901
12	0.891074	0.000876	0.001953	0.001056	0.116100	0.026894	0.000568	0.000283	7.637624	0.412184	3.585798	2.005262	0.97	False	0.039604	0.465347	0.257426	0.059406	0.019802
13	0.873987	0.001621	0.003267	0.001353	0.135719	0.023407	0.001112	0.000468	15.478022	0.510917	4.191749	2.377517	1.00	False	0.009901	0.168317	0.158416	0.009901	0.009901
14	0.863480	0.001100	0.001381	0.001214	0.061536	0.027973	0.000213	0.000156	2.886813	0.281927	1.900580	1.363811	0.97	False	0.297030	0.950495	0.940594	0.336634	0.673267
15	0.887117	0.000735	0.000798	0.000957	0.088362	0.047472	0.000177	0.000202	2.394520	0.402884	2.729088	0.876848	0.76	False	0.415842	1.000000	0.534653	0.554455	0.019802
16	0.905548	0.017899	0.004559	0.006921	0.221981	0.533353	0.002536	0.009860	33.070452	0.528455	6.855978	0.257256	1.00	False	0.009901	1.000000	0.059406	0.009901	0.009901

7. Select candidates and perform a global fit¶

The example keeps candidates with bliss_score > 0.5 and jointly fits their line profiles

In [41]:

Copied!





candidate_lines_evaluated_filtered = candidate_lines_evaluated[candidate_lines_evaluated.bliss_score>0.5]

lines_res_global = fit_global(
    candidate_lines_evaluated_filtered ,
    spectrum_snr,
    baseline_window=0.025,
    energy_min=0.62,
    energy_max=0.92)
candidate_lines_evaluated_filtered = candidate_lines_evaluated[candidate_lines_evaluated.bliss_score>0.5]

lines_res_global = fit_global(
    candidate_lines_evaluated_filtered ,
    spectrum_snr,
    baseline_window=0.025,
    energy_min=0.62,
    energy_max=0.92)

8. Assign tentative atomic identifications¶

add_most_probable_ion compares each fitted center with the bundled atomic line table, using a matching tolerance corresponding to 600 km s⁻¹.

In [42]:

Copied!

global_lines = add_most_probable_ion(lines_res_global, v_doppler_kms=600)
global_lines
global_lines = add_most_probable_ion(lines_res_global, v_doppler_kms=600)
global_lines

Out[42]:

	ion	doppler_kms	center	ecenter	sigma	esigma	amplitude	eamplitude	area	earea	...	relative_power	snr_peak	snr_area	bliss_score	response_feature	p_global_area	p_global_snr_area	p_global_snr_peak	p_global_ew	p_global_relative_power
0	o_viii	184.270376	0.653879	0.000445	0.001049	0.000325	0.118537	5.221817e-02	0.000312	0.000115	...	0.495620	3.176983	2.706087e+00	0.97	False	0.039604	0.019802	0.158416	0.039604	0.009901
1	o_vii	274.400682	0.666225	0.016019	0.000210	0.009376	0.224257	2.157240e+01	0.000118	0.006109	...	0.484043	6.646253	1.936383e-02	0.66	False	0.792079	0.079208	0.257426	0.712871	0.009901
2	ca_xix	-273.929155	0.674605	0.016097	0.000321	0.019077	0.046813	5.184791e+00	0.000038	0.001933	...	0.194666	1.499390	1.948460e-02	0.76	False	0.950495	0.851485	1.000000	0.950495	1.000000
3	o_vii	-332.861199	0.728938	0.022937	0.000426	0.103217	0.049639	1.201808e+01	0.000053	0.000042	...	0.230611	1.935118	1.272300e+00	0.58	False	1.000000	1.000000	0.990099	1.000000	0.940594
4	o_vii	593.392120	0.740211	0.000214	0.000619	0.000189	0.161857	5.387056e-02	0.000251	0.000073	...	0.512523	4.495800	3.416931e+00	0.97	False	0.059406	0.009901	0.158416	0.059406	0.009901
5	fe_viii	392.271865	0.744372	0.082832	0.000291	0.042923	0.034029	1.384725e+01	0.000025	0.006432	...	0.254779	1.135193	3.855514e-03	0.76	False	0.475248	0.227723	0.960396	0.435644	0.871287
6	o_viii	348.042966	0.775510	0.000407	0.001290	0.000387	0.095527	3.042818e-02	0.000309	0.000085	...	0.499244	3.225167	3.642996e+00	0.88	False	0.227723	1.000000	0.108911	0.217822	0.009901
7	ne_v	-401.648785	0.819390	0.005868	0.000280	0.012427	0.169272	1.541875e+01	0.000119	0.005539	...	0.353894	6.550498	2.142975e-02	0.66	False	0.851485	0.396040	0.425743	0.851485	0.178218
8	ne_v	376.580752	0.827600	0.006319	0.000289	0.012603	0.124321	1.080221e+01	0.000090	0.003895	...	0.223651	4.831511	2.311219e-02	0.66	False	0.960396	1.000000	0.990099	0.960396	0.970297
9	ne_vi	-219.317135	0.831891	593.009951	0.000170	105.218482	0.079800	9.345799e+04	0.000034	60.929040	...	0.097591	3.101262	5.589518e-07	0.58	False	1.000000	1.000000	1.000000	1.000000	1.000000
10	ne_iii	-473.687519	0.853653	0.001082	0.001093	0.001005	0.032650	2.918316e-02	0.000089	0.000075	...	0.171288	1.227339	1.188204e+00	0.76	False	0.297030	0.178218	0.940594	0.356436	1.000000
11	ne_vi	117.407447	0.881227	0.003902	0.000223	0.013684	0.187099	2.249822e+01	0.000104	0.006142	...	0.463636	5.778647	1.700594e-02	0.97	False	0.039604	0.019802	0.148515	0.049505	0.009901
12	ne_vii	-398.424552	0.892006	0.008568	0.000228	0.028999	0.088966	2.218657e+01	0.000051	0.006208	...	0.412184	2.747770	8.187055e-03	0.97	False	0.039604	0.465347	0.257426	0.059406	0.019802
13	ne_v	554.300097	0.874311	0.003020	0.000216	0.006740	0.314242	1.953341e+01	0.000170	0.005258	...	0.510917	9.705519	3.233121e-02	1.00	False	0.009901	0.168317	0.158416	0.009901	0.009901
14	ne_iv	425.989901	0.863364	0.009155	0.000272	0.034651	0.099100	1.517076e+01	0.000068	0.001743	...	0.281927	3.060756	3.878468e-02	0.97	False	0.297030	0.950495	0.940594	0.336634	0.673267
15	ne_vi	238.758474	0.887257	73681.001915	0.000181	29046.007763	0.017489	1.482360e+07	0.000008	5463.966070	...	0.402884	0.540149	1.454731e-09	0.76	False	0.415842	1.000000	0.534653	0.554455	0.019802
16	ne_vii	334.434293	0.901725	0.007456	0.000306	0.051880	0.147998	2.874242e+01	0.000113	0.002785	...	0.528455	4.570977	4.073182e-02	1.00	False	0.009901	1.000000	0.059406	0.009901	0.009901

17 rows × 21 columns

In [ ]: