Blind search on Chandra/HETG observation of Vela X-1¶

plt.figure(figsize=(9, 4.5))
plt.errorbar(spectrum.energy, spectrum.values, yerr=spectrum.uncertainties,fmt="b.",alpha=0.2)
plt.errorbar(spectrum.energy, spectrum.values,fmt="b.")
plt.xlabel("Energy (keV)")
plt.ylabel("counts s$^{-1}$ keV$^{-1}$")

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

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"
)

plt.figure(figsize=(9, 4.5))
plt.errorbar(spectrum_bins.energy, spectrum_bins.values, yerr=spectrum_bins.uncertainties,fmt="b.",alpha=0.2)
plt.errorbar(spectrum_bins.energy, spectrum_bins.values,fmt="b.")
plt.xlabel("Energy (keV)")
plt.ylabel("counts s$^{-1}$ keV$^{-1}$")

baseline_window = 0.4 #Default

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


# Figure ##########################################
fig, ax = plt.subplots(figsize=(9, 4.5))

# Spectrum
ax.errorbar(
    spectrum_bins.energy,
    spectrum_bins.values,
    yerr = spectrum_bins.uncertainties,
    fmt=".",
    color="blue",
    linewidth=0.1,
    alpha=0.2,
    label="Spectrum"
)

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

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

ax.legend()
fig.tight_layout()
plt.show()

en1, en2 = 0.5, 8

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)")

# Mark quick-search centers on a separately prepared three-bin spectrum.
fig, ax = plt.subplots(figsize=(9, 4.5))

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.5, 8),
    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)




ax.legend()
fig.tight_layout()
plt.show()

N_SIM = 100
SEED = 1
output_dir = "bliss_velax1_chandra"

config = BlindLineSearchConfig(
    en1=0.5,
    en2=8.0,
    energy_pad=0.0,
    rebin_method="bins",
    rebin_scale=2,  # combina 2 canales consecutivos
    baseline_window=0.4,
    max_sigma_line=0.1,
    num_synthetic_simulations=N_SIM,
    synthetic_seed=SEED,
    noise_model="poisson",
)

inicio = perf_counter()

pipeline = BlindLineSearchPipeline(config=config)
candidates_evaluated = pipeline.run(
    pha=pha,
    rmf=rmf,
    arf=arf,
    output_dir=output_dir,
)

tiempo = perf_counter() - inicio

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

plot_bliss_score(candidates_evaluated)
candidates_evaluated

plot_bliss_score(candidates_evaluated)

pvalue_columns = [
    "p_global_area",
    "p_global_snr_area",
    "p_global_snr_peak",
    "p_global_ew",
    "p_global_relative_power",
]

selected_candidates = candidates_evaluated.loc[
    candidates_evaluated[pvalue_columns].lt(0.01).any(axis=1)
    | candidates_evaluated["bliss_score"].gt(0.95)
].copy()

selected_candidates.head()

selected_candidates_global_fit = fit_global(
    selected_candidates,
    spectrum_bins,
    energy_min=en1, energy_max=en2,
    base=baseline_bins,
)

v_dop=800
add_most_probable_ion(selected_candidates_global_fit, v_doppler_kms=v_dop )

selected_candidates_global_fit_low = fit_global(
    selected_candidates,
    spectrum_bins,
    energy_min=1, energy_max=2,
    base=baseline_bins,
)

selected_candidates_global_fit_low.head()