{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# COSI Pulsar Analysis Tutorial: Phase Assignment and Filtering\n",
"\n",
"## Overview\n",
"This tutorial demonstrates a complete workflow for performing phase-resolved pulsar analysis using `cosipy` and custom time-domain tools. Using Data Challenge 2 (DC2) FITS data for the Crab pulsar and its accompanying albedo background, this notebook guides you through the process of preparing unbinned photon events, assigning rotational phases, and extracting specific pulse intervals.\n",
"\n",
"## Objectives\n",
"By the end of this notebook, you will learn how to:\n",
"* **Combine and Time-Slice Data:** Merge unbinned source and background datasets, and filter them down to specific Mission Elapsed Time (MET) windows for faster processing.\n",
"* **Assign Pulsar Phases:** Use a pulsar ephemeris parameter file (`.par`) to calculate and append the precise rotational phase (0.0 to 1.0) to individual photon events.\n",
"* **Visualize the Pulse Profile:** Generate a folded light curve (pulse profile) to visually inspect the on-pulse peaks and off-pulse baseline.\n",
"* **Perform Phase Selection:** Filter the phase-assigned data to isolate specific windows, such as the `on-pulse` and `off-pulse` regions, outputting the results into new FITS files for downstream spatial or spectral analysis.\n",
"\n",
"## Prerequisites\n",
"To run this notebook, ensure you have the following in your working directory:\n",
"* A valid ephemeris file (e.g., `crab.par`)\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Welcome to JupyROOT 6.28/12\n"
]
},
{
"data": {
"text/html": [
"14:28:06 WARNING The naima package is not available. Models that depend on it will not be functions.py : 43 available \n",
" WARNING The GSL library or the pygsl wrapper cannot be loaded. Models that depend on it functions.py : 65 will not be available. \n",
"14:28:07 WARNING The ebltable package is not available. Models that depend on it will not be absorption.py : 33 available \n",
"14:28:07 INFO Starting 3ML! __init__.py : 44 WARNING WARNINGs here are NOT errors __init__.py : 45 WARNING but are inform you about optional packages that can be installed __init__.py : 46 WARNING to disable these messages, turn off start_warning in your config file __init__.py : 47 WARNING Multinest minimizer not available minimization.py : 1218 WARNING PyGMO is not available minimization.py : 1228 WARNING Could not import plugin FermiLATLike.py. Do you have the relative instrument __init__.py : 126 software installed and configured? \n",
"14:28:08 WARNING No fermitools installed lat_transient_builder.py : 44 WARNING Env. variable OMP_NUM_THREADS is not set. Please set it to 1 for optimal __init__.py : 335 performances in 3ML \n",
" WARNING Env. variable MKL_NUM_THREADS is not set. Please set it to 1 for optimal __init__.py : 335 performances in 3ML \n",
" WARNING Env. variable NUMEXPR_NUM_THREADS is not set. Please set it to 1 for optimal __init__.py : 335 performances in 3ML \n",
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"filename = \"crab_with_albedo_10d_with_pulse_phase.fits\"\n",
"with fits.open(filename) as hdul:\n",
" pp = PlotPulseProfile(hdul[1].data)\n",
" pp.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The PhaseSelector is designed to isolate specific \"windows\" of a pulsar's rotation (e.g., the Main Pulse or the Interpulse). Unlike basic selectors that only allow one range, this class accepts a list of tuples. This allows you to select multiple distinct regions of the pulse profile in a single execution."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"# input_fits = \"crab_with_albedo_with_pulse_phase_column.fits\"\n",
"# output_fits = \"crab_with_albedo_with_pulse_phase_column_onpulse.fits\"\n",
"\n",
"input_fits = \"crab_with_albedo_10d_with_pulse_phase.fits\"\n",
"output_fits = \"crab_with_albedo_10d_with_pulse_phase_onpulse.fits\"\n",
"\n",
"# Each tuple represents (pstart, pstop).\n",
"# Values must be between 0.0 and 1.0.\n",
"intervals = [\n",
" (0.25, 0.40), \n",
" (0.55, 0.75), \n",
"]\n",
"\n",
"# Initialize the selector with your list of intervals\n",
"phase_sel = PhaseSelector(intervals=intervals)\n",
"\n",
"# Run the filter. This returns a structured NumPy array for immediate analysis\n",
"# and writes the filtered data to a new FITS file.\n",
"on_pulse_events = phase_sel.filter_events(\n",
" input_fits, \n",
" output_fits=output_fits\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Loading SpacecraftHistory from 20280301_3_month_with_orbital_info.ori...\n",
"\n",
"--- BEFORE CORRECTION ---\n",
"Original Livetime Array (first 5 bins): \n",
"[1. 1. 1. 1. 1.] s\n",
"\n",
"Phase Cut Total Width: 35%\n",
"Applying phase exposure correction...\n",
"\n",
"--- AFTER CORRECTION ---\n",
"Corrected Livetime Array (first 5 bins): \n",
"[0.35 0.35 0.35 0.35 0.35] s\n",
"\n",
"Exposure is now perfectly scaled to 35%!\n"
]
}
],
"source": [
"# --- CORRECT MISSION EXPOSURE ---\n",
"\n",
"# 1. Load the original spacecraft pointing history\n",
"ori_file = '20280301_3_month_with_orbital_info.ori'\n",
"print(f\"Loading SpacecraftHistory from {ori_file}...\")\n",
"sc_history = SpacecraftHistory.open(ori_file)\n",
"\n",
"# Let's look at the first few bins of the original livetime\n",
"print(\"\\n--- BEFORE CORRECTION ---\")\n",
"print(f\"Original Livetime Array (first 5 bins): \\n{sc_history._livetime_hist.contents[:5]}\")\n",
"\n",
"# Define the phase cut\n",
"intervals = [\n",
" (0.25, 0.40), \n",
" (0.55, 0.75), \n",
"]\n",
"\n",
"expected_fraction = sum([stop - start for start, stop in intervals])\n",
"print(f\"\\nPhase Cut Total Width: {expected_fraction * 100:.0f}%\")\n",
"\n",
"# 2. Apply the correction using the protocol\n",
"print(\"Applying phase exposure correction...\")\n",
"sc_history.update_ephemeris(timing_model, intervals)\n",
"\n",
"print(\"\\n--- AFTER CORRECTION ---\")\n",
"print(f\"Corrected Livetime Array (first 5 bins): \\n{sc_history._livetime_hist.contents[:5]}\")\n",
"print(f\"\\nExposure is now perfectly scaled to {expected_fraction * 100:.0f}%!\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Conclusion\n",
"\n",
"In this tutorial, we successfully processed raw, unbinned COSI data to assign rotational phases and extract specific phase intervals for the Crab pulsar and its corresponding albedo background. \n",
"\n",
"By isolating the `on-pulse` and `off-pulse` regions into separate FITS files, the data is now fully prepared for downstream analysis. Specifically, these newly created files can be plugged directly into `cosipy` to perform a **Spectral Fit** for the Crab pulsar. This will allow you to characterize its emission spectrum and analyze the background-subtracted source components."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "cosipy",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.13"
}
},
"nbformat": 4,
"nbformat_minor": 2
}