{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Fitting for a point source polarization with `LeakageLib`\n", "\n", "_Weights used:_\n", "* Spatial\n", "* Spectral\n", "\n", "This jupyter notebook fits for the polarization of a single simulated point source, with a photon and particle background. See other notebooks for how to use the more advanced weights." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[93m>>> PyXSPEC is not installed, you will no be able to use it.\u001b[0m\n" ] } ], "source": [ "import leakagelib" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, we need to load the data. See the next notebook for how to download and prepare real IXPE data. This notebook uses simulated data I made, which is located in `/docs/examples/data/ps`. This notebook also comes downloaded with LeakageLib, located at `/docs/examples/fit-point-source.ipynb`. If you run that notebook, the following load command should work." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ">>> Reading (in memory) /opt/homebrew/anaconda3/lib/python3.12/site-packages/ixpeobssim/caldb/ixpe/xrt/bcf/vign/ixpe_d1_obssim20240101_vign_v013.fits...\n" ] } ], "source": [ "datas = leakagelib.IXPEData.load_all_detectors_with_path(\"data\", \"ps\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "**Alternatives to** `IXPEData.load_all_detectors_with_path`\n", "\n", "- `IXPEData.load_all_detectors`. Searches all directories listed in the DATA_DIRECTORIES variable defined in settings.py for the given obs id, and load all detectors.\n", "- `IXPEData.load_all_detectors_with_path`. Provide a path to the observation and the function will load all detectors. This will not use the DATA_DIRECTORIES variable.\n", "- `IXPEData`. You must pass specific file names.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we center the data and cut to 2-8 keV." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "for data in datas:\n", " data.iterative_centroid_center()\n", " data.retain(data.evt_energies > 2)\n", " data.retain(data.evt_energies < 8)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "**Alternatives to** `IXPEData.iterative_centroid_center`\n", "\n", "- `IXPEData.explicit_center`: Center on a specific point, in pixels.\n", "- `IXPEData.iterative_centroid_center`: Center by iteratively zooming in on the average event. This finds central PSs rather well.\n", "
\n", "\n", "Now we have to tell the fitter how we want to do the fit. All this information is stored in the FitSettings object. First, we tell it to only consider events within 280 arcseconds of the center. More distant events will automatically be removed from the data sets." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "6282 events were cut for being outside the region of interest.\n" ] } ], "source": [ "settings = leakagelib.FitSettings(datas)\n", "settings.apply_circular_roi(280)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "**Alternatives to** `FitSettings.apply_circular_roi`\n", "\n", "- `IXPEData.apply_circular_roi`: Creates a circular ROI centered on the origin\n", "- `IXPEData.apply_roi_cut`: Adds a region file to the ROI (Or excludes it if you set `exclude=True`)\n", "- `IXPEData.apply_roi`: Apply an arbitrary ROI, which you pass in as an image. The image must have the same pixel size and dimensions as the fit sources.\n", "
\n", "\n", "You might be curious how many events are left. You can check this here:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "15894" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sum([len(data) for data in datas])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For this simplified scenario there are two components: a point source and a photon background. To create the point source," ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "settings.add_point_source() # Named \"src\" by default\n", "settings.fix_flux(\"src\", 1)\n", "settings.set_initial_qu(\"src\", (0.5,0))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Without loss of generality, we have defined the flux to be in units of the source flux.\n", "\n", "
\n", "\n", "You can use `FitSettings.set_initial_qu` to set the start point of the fitter. It's important to set the start point as close to the true value as possible to avoid falling into the wrong maximum of the likelihood.\n", "\n", "
\n", "\n", "Now we create the background models. These background sources are assumed to be spatially uniform and already have spatial weights applied. " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "settings.add_background() # Named \"bkg\" by default\n", "settings.fix_qu(\"bkg\", (0, 0)) # Assume an unpolarized background" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To use energy weights, we must assign a spectrum. We create functions which take in the energy in keV and return the photon spectrum (dN/dE). Normalization doesn't matter. The ARF and RMF are automatically included, but you can disable them with the `use_arf` and `use_rmf` arguments." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ">>> Reading (in memory) /opt/homebrew/anaconda3/lib/python3.12/site-packages/ixpeobssim/caldb/ixpe/gpd/cpf/arf/ixpe_d1_obssim20240101_v013.arf...\n", ">>> Using cached xEffectiveArea object at /opt/homebrew/anaconda3/lib/python3.12/site-packages/ixpeobssim/caldb/ixpe/gpd/cpf/arf/ixpe_d1_obssim20240101_v013.arf...\n" ] } ], "source": [ "settings.set_spectrum(\"bkg\", lambda e: e**-2.5) # Gamma=2.5, unabsorbed powerlaw\n", "settings.set_spectrum(\"src\", lambda e: e**-1.5) # Gamma=1.5, unabsorbed powerlaw" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have assigned all the weights. Now we are ready to create the fitter..." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Data set ps DU 1 had no exposure map loaded. Please load an exposure map if you are fitting to events in the vignetted portion.\n" ] } ], "source": [ "fitter = leakagelib.Fitter(settings)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "... and run the fit" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "FitResult:\n", "\tq (src) = 0.4697 +/- 0.0725\n", "\tu (src) = -0.0260 +/- 0.0721\n", "\tf (bkg) = 2.1813 +/- 0.0423\n", "\n", "Polarization:\n", "\tPD (src): 0.4704 +/- 0.0725\n", "\tPA (src): -1.5836 deg +/- 4.3936\n", "Likelihood 17886.60119278403, dof 15891\n", "Optimization terminated successfully." ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "result = fitter.fit()\n", "result" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The true polarization of this simulated source was q = 0.5, u = 0.0. The fit result agrees." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "**Alternatives to** `Fitter.fit`\n", "\n", "- `Fitter.fit`: Fit by numerically maximizing the likelihood and get Gaussian uncertainties by taking the Hessian and using Laplace's approximation.\n", "- `IXPEData.fit_mcmc`: Do a full MCMC fit.\n", "\n", "
\n", "\n", "You can access `result.params` and `result.sigmas` to get the best-fit values and uncertainties. Also `result.cov` to get the full covariance matrix. The covariance matrix entries are stored in the same order as the `result.parameter_names` array." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "===Results as human-readable dictionaries===\n", "Parameters {('q', 'src'): 0.46965223111815946, ('u', 'src'): -0.025987502310235725, ('f', 'bkg'): 2.1813398855067225}\n", "Uncertainties {('q', 'src'): 0.07247836727152544, ('u', 'src'): 0.07213767036302343, ('f', 'bkg'): 0.0423476862090614}\n", "\n", "===Results as arrays===\n", "Parameter names [('q', 'src'), ('u', 'src'), ('f', 'bkg')]\n", "Best fit values [ 0.46965223 -0.0259875 2.18133989]\n", "Covariance [[ 5.25311372e-03 -1.73106799e-05 0.00000000e+00]\n", " [-1.73106799e-05 5.20384349e-03 0.00000000e+00]\n", " [ 0.00000000e+00 0.00000000e+00 1.79332653e-03]]\n" ] } ], "source": [ "print(\"===Results as human-readable dictionaries===\")\n", "print(\"Parameters\", result.params)\n", "print(\"Uncertainties\", result.sigmas)\n", "print()\n", "print(\"===Results as arrays===\")\n", "print(\"Parameter names\", result.parameter_names)\n", "print(\"Best fit values\", result.best_params)\n", "print(\"Covariance\", result.cov)" ] } ], "metadata": { "kernelspec": { "display_name": "base", "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.2" } }, "nbformat": 4, "nbformat_minor": 2 }