"""
IXPE Data Visualization Module.
This module provides plotting functions for IXPE spectral and polarimetric data,
including Stokes I/Q/U spectra and Q-U plane visualizations with confidence contours.
"""
import numpy
import matplotlib.pyplot as plt
from threeML.io.plotting.post_process_data_plots import display_spectrum_model_counts
from ixpepy.utils.stokes import get_polarization_angle, get_polarization_degree
from scipy.stats import gaussian_kde
[docs]
def plot_I_spectrum(like, emin=2., emax=8., ymin=None, ymax=None, basename='ixpe'):
"""
Plot the Stokes I (total intensity) spectrum for all three IXPE detector units.
Parameters
----------
like : JointLikelihood or similar
Likelihood object containing the fitted model and data
emin : float
Minimum energy for x-axis (keV)
emax : float
Maximum energy for x-axis (keV)
ymin : float, optional
Minimum count rate for y-axis. Default is None (auto-scale)
ymax : float, optional
Maximum count rate for y-axis. Default is None (auto-scale)
basename : str, optional
Base name for the IXPE plugins. Default is 'ixpe'
Returns
-------
matplotlib.figure.Figure
Figure object containing the Stokes I spectrum plot
"""
figI = display_spectrum_model_counts(like, data=[f'{basename}_DU1_I',f'{basename}_DU2_I',f'{basename}_DU3_I'], show_background=False)
figI.axes[0].set_ylim(ymin, ymax)
figI.axes[0].set_xlim(emin, emax)
plt.tight_layout()
return figI
[docs]
def plot_QU_spectra(like, emin=2., emax=8., ymin=None, ymax=None, basename='ixpe'):
"""
Plot the Stokes Q and U spectra for all three IXPE detector units.
Parameters
----------
like : JointLikelihood or similar
Likelihood object containing the fitted model and data
emin : float
Minimum energy for x-axis (keV)
emax : float
Maximum energy for x-axis (keV)
ymin : float, optional
Minimum count rate for y-axis. Default is None (auto-scale)
ymax : float, optional
Maximum count rate for y-axis. Default is None (auto-scale)
basename : str, optional
Base name for the IXPE plugins. Default is 'ixpe'
Returns
-------
tuple of matplotlib.figure.Figure
(figQ, figU) - Figure objects for Stokes Q and U spectra
"""
figQ = display_spectrum_model_counts(like, data=[f'{basename}_DU1_Q', f'{basename}_DU2_Q', f'{basename}_DU3_Q'], min_rate=-10)
figQ.axes[0].set_yscale('linear')
figQ.axes[0].set_ylim(ymin, ymax)
figQ.axes[0].set_xlim(emin, emax)
plt.tight_layout()
figU = display_spectrum_model_counts(like, data=[f'{basename}_DU1_U', f'{basename}_DU2_U', f'{basename}_DU3_U'], min_rate=-10)
figU.axes[0].set_yscale('linear')
figU.axes[0].set_ylim(ymin, ymax)
figU.axes[0].set_xlim(emin, emax)
plt.tight_layout()
return figQ, figU
[docs]
def cdf_2d(sigma):
"""
Convert sigma confidence level to 2D cumulative distribution function value.
Parameters
----------
sigma : float
Confidence level in sigmas (e.g., 1, 2, 3)
Returns
-------
float
CDF value corresponding to the sigma level
"""
return numpy.exp(-(sigma**2)/2)
[docs]
def fmt(x):
"""
Format contour level value as sigma string.
Parameters
----------
x : float
CDF value
Returns
-------
str
Formatted string with sigma symbol
"""
s = numpy.sqrt(-2.0 * numpy.log(x))
ss = f"{s:.1f}"
if ss.endswith("0"):
ss=f"{s:.0f}"
return rf"{ss} $\sigma$"
[docs]
def fmt_s(x):
"""
Format contour level value as sigma string (alternative).
Parameters
----------
x : float
CDF value
Returns
-------
str
Formatted string with sigma symbol
"""
s = numpy.sqrt(-2.0 * numpy.log(x))
ss = f"{s:.1f}"
if ss.endswith("0"):
ss=f"{s:.0f}"
return rf"{ss} $\sigma$"
[docs]
def fmt_p(x):
"""
Format contour level value as percentage string.
Parameters
----------
x : float
CDF value (between 0 and 1)
Returns
-------
str
Formatted percentage string
"""
s=100*(1.0-x)
ss = f"{s:.1f}"
if ss.endswith("0"):
ss=f"{s:.0f}"
return rf"{ss} $\%$"
[docs]
def plot_pdf_2d(x, y, sigma=False, ax=None):
"""
Plot 2D probability density function contours for Q-U samples.
Uses kernel density estimation to create smooth confidence contours
from discrete samples.
Parameters
----------
x : array-like
X-axis values (e.g., Q samples)
y : array-like
Y-axis values (e.g., U samples)
sigma : bool, optional
If True, use sigma levels (1, 2, 3σ). If False, use percentage
levels (50%, 90%, 99%). Default is False.
ax : matplotlib.axes.Axes, optional
Axes object to plot on. If None, uses current axes. Default is None.
Returns
-------
matplotlib.contour.QuadContourSet
Contour set object
"""
if ax is None:
ax = plt.gca()
k = gaussian_kde(numpy.vstack([x, y]))
xi, yi = numpy.mgrid[x.min():x.max():x.size ** 0.5 * 1j, y.min():y.max():y.size ** 0.5 * 1j]
zi = k(numpy.vstack([xi.flatten(), yi.flatten()]))
zi = (zi - zi.min()) / (zi.max() - zi.min())
levels=[]
if sigma:
levels = (cdf_2d(3), cdf_2d(2), cdf_2d(1))
fmt = fmt_s
else:
levels = (1-0.99, 1-0.90, 1-0.50)
fmt = fmt_p
CS = ax.contour(xi, yi, zi.reshape(xi.shape), colors='k', levels=levels)
return CS
[docs]
def plot_QU_polar(q, u, rmax=None, nbins=50, mdp_99=0, sigma=False, manual_locations=False):
"""
Create a polar plot of polarization degree vs. polarization angle.
Visualizes the posterior distribution in polar coordinates with
2D histogram and confidence contours.
Parameters
----------
q : RandomVariate or similar
Stokes Q parameter samples (normalized)
u : RandomVariate or similar
Stokes U parameter samples (normalized)
rmax : float, optional
Maximum radius (polarization degree) for plot. If None, auto-scales.
Default is None.
nbins : int, optional
Number of bins for 2D histogram. Default is 50.
mdp_99 : float, optional
Minimum Detectable Polarization at 99% confidence to overlay.
Default is 0 (no MDP line).
sigma : bool, optional
If True, label contours with sigma levels. If False, use percentages.
Default is False.
manual_locations : bool, optional
If True, allows manual positioning of contour labels. Default is False.
Returns
-------
matplotlib.figure.Figure
Figure object containing the polar plot
"""
pol_ang = get_polarization_angle(q, u)
pol_deg = get_polarization_degree(q, u)
fig, ax = plt.subplots(subplot_kw=dict(projection='polar'))
if mdp_99 > 0:
_x = numpy.linspace(0, 2*numpy.pi, 100)
_y = mdp_99 * numpy.ones_like(_x)
plt.plot(_x, _y, 'r:', label='MDP$_{99}$')
x = numpy.radians(pol_ang.samples)
y = pol_deg.samples
if rmax is None:
rmax = max(y)*1.1
_ = plt.hist2d(x, y, bins=nbins,range=[[numpy.radians(0),numpy.radians(90)],[0,rmax]], cmap=plt.cm.Reds)
CS = plot_pdf_2d(x, y, sigma=sigma, ax=ax)
plt.clabel(CS, fmt=fmt, colors='k', fontsize=14,manual=manual_locations)
plt.grid()
plt.draw()
return fig
[docs]
def plot_QU(q, u, nbins=40, mdp_99=0, sigma=False, pcube=False, manual_locations=False):
"""
Create a Cartesian Q-U plane plot with confidence contours.
Visualizes the posterior distribution of Stokes Q and U parameters
with 2D histogram, confidence contours, and polarization degree circles.
Parameters
----------
q : RandomVariate or similar
Stokes Q parameter samples (normalized by I)
u : RandomVariate or similar
Stokes U parameter samples (normalized by I)
nbins : int, optional
Number of bins for 2D histogram. Default is 40.
mdp_99 : float, optional
Minimum Detectable Polarization at 99% confidence to overlay
as a circle. Default is 0 (no MDP circle).
sigma : bool, optional
If True, label contours with sigma levels. If False, use percentages.
Default is False.
pcube : dict or False, optional
If dict, overlay PCUBE results with keys 'i', 'q', 'u', 'err'.
If False, no PCUBE overlay. Default is False.
manual_locations : bool, optional
If True, allows manual positioning of contour labels. Default is False.
Returns
-------
matplotlib.figure.Figure
Figure object containing the Q-U plane plot
Notes
-----
The plot includes:
- 2D histogram of Q-U samples (orange colormap)
- Best-fit value marked with a cross
- Confidence contours (50%, 90%, 99% or 1σ, 2σ, 3σ)
- Polarization degree circles at 5%, 10%, 15%, ... intervals
- Polarization angle lines at 30° intervals
- Optional MDP circle and PCUBE comparison
"""
fig, ax = plt.subplots(figsize=(6.4, 6.2))
ax.spines['top'].set_visible(True)
ax.spines['right'].set_visible(True)
if mdp_99 > 0:
_x = numpy.linspace(0, 2*numpy.pi, 100)
_y = mdp_99 * numpy.ones_like(_x)
plt.plot(_x,_y,'r:',label='MDP$_{99}$')
x = q.samples
y = u.samples
qmin, qmax = x.min(), x.max()
umin, umax = y.min(), y.max()
xmin = min(0.8 * qmin, -0.02)
xmax = max(1.2 * qmax, 0.02)
ymin = min(0.8 * umin, -0.02)
ymax = max(1.2 * umax, 0.02)
plt.hist2d(x, y, bins=nbins, range=[[xmin, xmax],[ymin, ymax]], cmap=plt.cm.Oranges)
plt.plot(q.value, u.value, '+', lw=2, color='k', label='3ML [ixpepy]')
if pcube:
pcube_i = pcube['i']
pcube_q, pcube_u = pcube['q'] / pcube_i, pcube['u'] / pcube_i
pcube_err = pcube['err'] / pcube_i
plt.plot(pcube_q, pcube_u, 'x', color='green', linewidth=2, label='PCUBE')
c = plt.Circle((pcube_q, pcube_u), radius=pcube_err, facecolor='none', edgecolor='green')
plt.gca().add_artist(c)
qumax = max(qmax, umax)
if qumax < 0.2:
step = 0.05
elif qumax < 0.5:
step = 0.1
else:
step = 0.2
for i in numpy.arange(1, numpy.ceil(qumax/step)):
r = i*step
c = plt.Circle((0,0), radius=r, facecolor='none', edgecolor='k', linewidth=1, linestyle='--', alpha=0.5)
plt.gca().add_artist(c)
_x = r * numpy.cos(numpy.sign(x.average) * numpy.pi / 4.)
_y = r * numpy.sin(numpy.sign(y.average) * numpy.pi / 4.)
plt.text(_x, _y, '%d%%' % int(r*100), fontsize=12, color='k', alpha=1)
plt.hlines(0, -1, 1, linewidth=1, color='k', linestyle='--', alpha=0.5)
plt.vlines(0, -1, 1, linewidth=1, color='k', linestyle='--', alpha=0.5)
r = 0.15
for phi in numpy.arange(0., 360., 30):
x_ = r * numpy.cos(numpy.radians(phi))
y_ = r * numpy.sin(numpy.radians(phi))
plt.plot([0, 2. * x_], [0, 2. * y_], color='gray', ls='dashed', lw=1)
ax.set_xlabel('Q/I', size=18)
ax.set_ylabel('U/I', size=18)
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
CS = plot_pdf_2d(x, y, sigma=sigma, ax=ax)
plt.clabel(CS, fmt=fmt, colors='k', fontsize=14, manual=manual_locations)
plt.legend(loc=8, fontsize=16)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.tight_layout()
ax = plt.gca()
ax.set_aspect('equal', adjustable='box')
plt.draw()
return fig