304.1. PSF analysis for cosmology

304.1. Cosmological Applications of Point Spread Function Analysis¶

For the Rubin Science Platform at data.lsst.cloud.
Data Release: Data Preview 1
Container Size: large
LSST Science Pipelines version: r29.2.0
Last verified to run: 2025-09-02
Repository: github.com/lsst/tutorial-notebooks

Learning objective: To use single-epoch and coadded PSFs in cosmological analyses.

LSST data products: source, object, visit_images, coadd

Packages: treecorr, lsst.daf.butler

Credit: Originally developed by Andrés A. Plazas Malagón and the Rubin Community Science team. Sections 3, 4, and 5 on cluster lensing are based on a notebook by Céline Combet (plus others in the cluster commissioning team). Section 6 on "PSF residuals in focal plane" is based on a notebook by Pierre-François Leget. Please consider acknowledging them if this notebook is used for the preparation of journal articles, software releases, or other notebooks.

Get Support: Everyone is encouraged to ask questions or raise issues in the Support Category of the Rubin Community Forum. Rubin staff will respond to all questions posted there.

1. Introduction¶

Accurate modeling of the Point Spread Function (PSF) is a fundamental component of modern cosmological analyses. In particular, weak gravitational lensing—both in the form of cosmic shear measurements and galaxy cluster mass calibration—relies on precise knowledge of the PSF to infer subtle distortions in the shapes of background galaxies. Imperfections in PSF modeling can lead to significant biases in shear measurements and, consequently, in cosmological parameter estimation.

This notebook explores several diagnostic tools and validation techniques for PSF modeling using early commissioning data from LSSTComCam, focusing on the galaxy cluster Abell 360.

Related tutorials: The 200-level tutorials on the PSF in visit_image and deep_coadd images. The 100-level tutorials on how to use the Butler.

1.1. Import packages¶

In [1]:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
from scipy.stats import binned_statistic
from astropy.coordinates import SkyCoord
import treecorr

from lsst.daf.butler import Butler
from lsst.rsp import get_tap_service

1.2. Define parameters and functions¶

Function: get_ellipticity

Compute the star and PSF model ellipticity components (e1, e2) from second moments.

In [2]:

def get_ellipticity(catalog, band, obj_type='star'):
    """Compute the ellipticity components (e1, e2) from second moments
    for stars or PSF models.

    Parameters
    ----------
    catalog : `astropy.table.Table`
        Catalog containing second moment columns for the
        specified band and object type.
    band : `str`
        Photometric band (e.g., 'i', 'r') used to construct the column names.
    obj_type : `str`, optional
        Type of object: 'star' (default) or 'PSF'. Determines whether to read
        columns like '{band}_ixx' or '{band}_ixxPSF'.

    Returns
    -------
    e1 : `np.array`
        First ellipticity component, defined as (Ixx - Iyy) / (Ixx + Iyy).
    e2 : `np.array`
        Second ellipticity component, defined as 2 * Ixy / (Ixx + Iyy).

    Raises
    ------
    ValueError
        If `obj_type` is not one of 'star' or 'PSF'.
    """
    if obj_type == 'star':
        suffix = ''
    elif obj_type == 'PSF':
        suffix = 'PSF'
    else:
        raise ValueError("obj_type must be either 'star' or 'PSF'.")

    ixx = catalog[f'{band}_ixx{suffix}']
    iyy = catalog[f'{band}_iyy{suffix}']
    ixy = catalog[f'{band}_ixy{suffix}']
    denom = ixx + iyy
    e1 = (ixx - iyy) / denom
    e2 = 2. * ixy / denom
    return e1, e2

Function: compute_et

Computes the tangential ellipticity component around a point.

In [3]:

def compute_et(cat, e1_star, e2_star, e1_psf, e2_psf,
               center_ra, center_dec, r_bins):
    """Compute the tangential and cross ellipticity components
    relative to a central coordinate.

    Parameters
    ----------
    cat : `astropy.table.Table`
        Catalog containing object positions with 'coord_ra' and
        'coord_dec' columns (in degrees).
    e1_star : `np.ndarray`
        Measured e1 ellipticity component of the stars.
    e2_star : `np.ndarray`
        Measured e2 ellipticity component of the stars.
    e1_psf : `np.ndarray`
        PSF model e1 ellipticity component at the star positions.
    e2_psf : `np.ndarray`
        PSF model e2 ellipticity component at the star positions.
    center_ra : `float`
        Right ascension of the center point, in degrees.
    center_dec : `float`
        Declination of the center point, in degrees.
    r_bins : `array_like`
        Bin edges (in arcminutes) for computing radial statistics.

    Returns
    -------
    mid : `np.ndarray`
        Midpoints of the radial bins (arcminutes).
    e_t_mean : `np.ndarray`
        Mean tangential ellipticity residual in each radial bin.
    e_t_err : `np.ndarray`
        Standard error on the tangential ellipticity residual in
        each bin.
    count : `np.ndarray`
        Number of objects in each radial bin.

    Notes
    -----
    `e_t` is the tangential ellipticity, `e_x` is the cross ellipticity,
    `pos_angle` is the position angle toward the cluster center
    (with North up and East to the left),
    and `e_1` and `e_2` are the original ellipticity components.

    The position angle `pos_angle` is measured from the negative RA direction
    (i.e., West) and increases counterclockwise. However, in `astropy`,
    the position angle is measured from North and increases toward the East.
    To match the convention used in the ellipticity decomposition,
    add `pi/2` to the position angle computed by astropy.
    """

    coord0 = SkyCoord(center_ra, center_dec, frame='icrs', unit='deg')
    coord1 = SkyCoord(cat['coord_ra'], cat['coord_dec'],
                      frame='icrs', unit='deg')

    pos_angle = coord0.position_angle(coord1).rad + np.pi / 2.
    r = coord0.separation(coord1).arcmin

    e1 = e1_star - e1_psf
    e2 = e2_star - e2_psf

    e_t = -e1 * np.cos(2. * pos_angle) - e2 * np.sin(2. * pos_angle)
    mid = 0.5 * (r_bins[1:] + r_bins[:-1])
    count, _, _ = binned_statistic(r, e_t, bins=r_bins, statistic='count')
    e_t_mean, _, _ = binned_statistic(r, e_t, bins=r_bins)
    e_t_std, _, _ = binned_statistic(r, e_t, bins=r_bins, statistic='std')
    e_t_err = e_t_std / np.sqrt(count)

    return mid, e_t_mean, e_t_err, count

Instantiate the Butler with the appropiate DP1 repository and collection.

In [4]:

butler = Butler('dp1', collections='LSSTComCam/DP1')

Instantiate the Table Access Protocol (TAP) service.

In [5]:

rsp_tap = get_tap_service("tap")
assert rsp_tap is not None

2. PSF size correlation function: source table.¶

Retrieve information from the source table for a particular visit.

First, retrieve the reference object around a chosen coordinate on the sky, for one particular visit. Use the center location (in degrees) of cluster Abell 360. Use the coordinates of the brightest cluster galaxy (BCG) as the cluster center (in degrees).

In [6]:

ra_bcg, dec_bcg = 37.862, 6.98

In [7]:

my_band = 'i'

In [8]:

visitimage_refs = butler.query_datasets(
    'visit_image',
    where=("visit.region OVERLAPS POINT(ra, dec) "
           f"and BAND='{my_band}'"),
    bind={'ra': ra_bcg, 'dec': dec_bcg},
    order_by=['visit', 'detector'],
    with_dimension_records=True,
    limit=1
)

In [9]:

visit_id = visitimage_refs[0].dataId['visit']

Use the TAP service to choose bright objects that were used to produce the PSF model.

In [10]:

query = (
    "SELECT sourceId, x, y, "
    "coord_ra, coord_dec, "
    "ixxPSF, iyyPSF, "
    "psfFlux, "
    "calib_psf_used, "
    "calibFlux, "
    "extendedness, "
    "pixelFlags_suspectCenter "
    "FROM dp1.Source "
    f"WHERE visit = {visit_id} "
    "AND calib_psf_used = 1 "
    "AND calibFlux > 360 "
    "AND extendedness = 0.0 "
    "AND pixelFlags_suspectCenter = 0 "
    "ORDER BY sourceId "
)

In [11]:

job = rsp_tap.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)
if job.phase == 'ERROR':
    job.raise_if_error()
assert job.phase == 'COMPLETED'
psf_used_table_visit = job.fetch_result().to_table()

Job phase is COMPLETED

Calculate the trace radius of the PSF, in pixels, from the PSF model second moments.

In [12]:

psf_used_table_visit['size_PSF'] = np.sqrt(
    (psf_used_table_visit['ixxPSF'] + psf_used_table_visit['iyyPSF']) / 2
)

Temporarily store the sources right ascention and declination coordinates in the src_ra and src_dec variables.

In [13]:

src_ra = psf_used_table_visit['coord_ra']
src_dec = psf_used_table_visit['coord_dec']

Use the sofwatre treecorr to calculate the correlation functions.

Define treecorr's catalog and adjust a few of its parameters, calculate the two-point correlation function of the size, and retrieve it, along with the angular seprations.

In [14]:

cat = treecorr.Catalog(ra=src_ra, dec=src_dec,
                       k=psf_used_table_visit['size_PSF']
                       - np.mean(psf_used_table_visit['size_PSF']),
                       ra_units='deg', dec_units='deg')
kk_config = {'max_sep': .06, 'min_sep': .0001, 'nbins': 12}
kk = treecorr.KKCorrelation(kk_config)
kk.process(cat)
xi_source = kk.xi
bins_source = kk.rnom

Plot the PSF size as a function of sky coordinates, and the correlation function.

In [15]:

_, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4.5),
                             gridspec_kw={'wspace': .3})

ax1.set_title(f'PSF size. Visit: {visit_id}')
scatter_plot = ax1.scatter(src_ra, src_dec,
                           c=psf_used_table_visit['size_PSF'],
                           s=1, cmap='cividis')
ax1.set_xlabel('RA [deg]')
ax1.set_ylabel('DEC [deg]')
plt.colorbar(scatter_plot, ax=ax1, label='[pixels]')

ax2.set_title(f'PSF size correlation. Visit: {visit_id}')
ax2.plot(np.degrees(bins_source), xi_source*1e4, 'o-', color='darkblue')
ax2.axhline(0, linestyle='--', color='lightgrey')
ax2.set_xscale('log')
ax2.set_ylabel(r'$\xi \times 10^{4}$', labelpad=-12)
ax2.set_xlabel(r'$\theta$ [degree]')

plt.show()

Figure 1: Left: spatial distribution of PSF size (pixels) for a region in a particular visit, using the source table .Right: PSF size correlation function in the same field as the panel on the left.

3. PSF whisker plots in Abell 360 with object table.¶

Whisker plots are visualizations that represent the shape and orientation of point spread functions (PSFs) or galaxy ellipticities as small line segments ("whiskers") at their sky positions. Each whisker’s direction encodes the ellipticity angle, and its length reflects the ellipticity magnitude.

In weak lensing, whisker plots are crucial because they help visualize spatial patterns and systematics in the PSF across the field of view. Since accurate PSF modeling is essential for measuring weak gravitational lensing shear, these plots are an intuitive way to spot residual anisotropies or discontinuities that could bias shear measurements.

Select stars in the object table in a radius of 0.75 degrees around the cluster BCG.

Define a template query to select stars that were used in PSF modeling, as well as "reserved" stars, which were excluded from the modeling process to test the fit.

In [16]:

query_template = """
SELECT
    objectId, x, y,
    refExtendedness,
    {my_band}_calib_psf_used,
    {my_band}_calib_psf_reserved,
    {my_band}_pixelFlags_inexact_psfCenter,
    {my_band}_calibFlux,
    {my_band}_centroid_x, {my_band}_centroid_y,
    coord_ra, coord_dec,
    {my_band}_ixx, {my_band}_ixxPSF,
    {my_band}_ixy, {my_band}_ixyPSF,
    {my_band}_iyy, {my_band}_iyyPSF,
    {my_band}_psfFlux, {my_band}_psfFluxErr,
    {my_band}_extendedness
FROM dp1.Object
WHERE
    CONTAINS(POINT('ICRS', coord_ra, coord_dec),
             CIRCLE('ICRS', {ra_bcg}, {dec_bcg}, 0.75)) = 1
    AND refExtendedness = 0.0
    AND {psf_flag} = 1
    AND {my_band}_pixelFlags_inexact_psfCenter = 0
    {flux_cut}
ORDER BY objectId
"""

Define a list with the charcateristics for the used and reserved stars queries.

In [17]:

query_types = [
    ("psf_used_table", f"{my_band}_calib_psf_used", f"AND {my_band}_calibFlux > 360"),
    ("psf_reserved_table", f"{my_band}_calib_psf_reserved", "")
]

Loop over the query types list and use the TAP service to obtain the two tables.

In [18]:

results = {}
for name, psf_flag, flux_cut in query_types:
    query = query_template.format(
        my_band=my_band,
        ra_bcg=ra_bcg,
        dec_bcg=dec_bcg,
        psf_flag=psf_flag,
        flux_cut=flux_cut
    )
    job = rsp_tap.submit_job(query)
    job.run()
    job.wait(phases=['COMPLETED', 'ERROR'])
    print('Job phase is', job.phase)
    if job.phase == 'ERROR':
        job.raise_if_error()
    assert job.phase == 'COMPLETED'
    results[name] = job.fetch_result().to_table()

Job phase is COMPLETED

Job phase is COMPLETED

In [19]:

psf_used_table = results["psf_used_table"]
psf_reserved_table = results["psf_reserved_table"]

The reserved stars are about 10 percent of the total stars in the region.

In [20]:

len(psf_used_table)

Out[20]:

2617

In [21]:

len(psf_reserved_table)

Out[21]:

302

Calculate ellipticities (magnitude and orientations), and their Cartesian components for the whisker plots for the used and reserved PSF models and observed stars.

1. For the PSF model, at the location of the used stars.

In [22]:

e1_psf_used, e2_psf_used = get_ellipticity(psf_used_table,
                                           my_band,
                                           obj_type='PSF')
e_psf_used = np.sqrt(e1_psf_used * e1_psf_used + e2_psf_used * e2_psf_used)
theta_psf_used = 0.5 * np.arctan2(e2_psf_used, e1_psf_used)

cx_psf_used = e_psf_used * np.cos(theta_psf_used)
cy_psf_used = e_psf_used * np.sin(theta_psf_used)

2. For the PSF model, at the location of the reserved stars.

In [23]:

e1_psf_reserved, e2_psf_reserved = get_ellipticity(psf_reserved_table,
                                                   my_band,
                                                   obj_type='PSF')
e_psf_reserved = np.sqrt(e1_psf_reserved * e1_psf_reserved + e2_psf_reserved * e2_psf_reserved)
theta_psf_reserved = 0.5 * np.arctan2(e2_psf_reserved, e1_psf_reserved)

cx_psf_reserved = e_psf_reserved * np.cos(theta_psf_reserved)
cy_psf_reserved = e_psf_reserved * np.sin(theta_psf_reserved)

3. For the observed stars that were used in PSF modeling.

In [24]:

e1_star_used, e2_star_used = get_ellipticity(psf_used_table,
                                             my_band,
                                             obj_type='star')
e_star_used = np.sqrt(
    e1_star_used * e1_star_used + e2_star_used * e2_star_used
)
theta_star_used = 0.5 * np.arctan2(e2_star_used, e1_star_used)

cx_star_used = e_star_used * np.cos(theta_star_used)
cy_star_used = e_star_used * np.sin(theta_star_used)

4. For the observed stars that were reserved from the PSF modeling process.

In [25]:

e1_star_reserved, e2_star_reserved = get_ellipticity(psf_reserved_table,
                                                     my_band,
                                                     obj_type='star')
e_star_reserved = np.sqrt(
    e1_star_reserved * e1_star_reserved + e2_star_reserved * e2_star_reserved
)
theta_star_reserved = 0.5 * np.arctan2(e2_star_reserved, e1_star_reserved)

cx_star_reserved = e_star_reserved * np.cos(theta_star_reserved)
cy_star_reserved = e_star_reserved * np.sin(theta_star_reserved)

Define ellipticity residuals for used stars.

In [26]:

e1_residual_used = e1_star_used - e1_psf_used
e2_residual_used = e2_star_used - e2_psf_used

theta_residual_used = 0.5 * np.arctan2(e2_residual_used,
                                       e1_residual_used)
e_residual_used = np.sqrt(
    e1_residual_used * e1_residual_used + e2_residual_used * e2_residual_used
)

cx_residual_used = e_residual_used * np.cos(theta_residual_used)
cy_residual_used = e_residual_used * np.sin(theta_residual_used)

In [27]:

cx_star_used = np.asarray(cx_star_used)
cy_star_used = np.asarray(cy_star_used)
cx_psf_used = np.asarray(cx_psf_used)
cy_psf_used = np.asarray(cy_psf_used)
cx_residual_used = np.asarray(cx_residual_used)
cy_residual_used = np.asarray(cy_residual_used)

Whisker plots for measured stars and PSF models, and their residuals.

In [28]:

scale = 0.6
ref_norm = scale * 0.1
ref_x, ref_y = 38.4, 6.27
ra_limits = [38.7, 37.3]
dec_limits = [6.25, 7.6]
x_centroid = np.asarray(psf_used_table['coord_ra'])
y_centroid = np.asarray(psf_used_table['coord_dec'])

Prepare whisker data. Switch the sign of the abscisa component for the plots to display RA increasing to the left.

In [29]:

whiskers = [
    (scale * -cx_star_used, scale * cy_star_used,
     'Star ellipticity (from i_ixx, i_ixy, i_iyy)'),
    (scale * -cx_psf_used, scale * cy_psf_used,
     'PSF model ellipticity (from i_ixxPSF, i_ixyPSF, i_iyyPSF)'),
    (scale * -cx_residual_used, scale * cy_residual_used,
     'Star - PSF ellipticity residuals')
]

In [30]:

fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(15, 5))

for i, ax in enumerate(axes):
    cx, cy, title = whiskers[i]
    ax.quiver(x_centroid, y_centroid, cx, cy, angles='xy', color='black',
              scale_units='xy', scale=1, headlength=0, headwidth=0,
              headaxislength=0)
    ax.scatter([ra_bcg], [dec_bcg], marker='+', s=100, c='darkviolet')
    ax.add_patch(Circle((ra_bcg, dec_bcg), 0.5, color='darkviolet',
                        fill=False, linewidth=1))

    ax.quiver(ref_x, ref_y, ref_norm, 0, angles='xy', color='red',
              scale_units='xy', scale=1, headlength=0, headwidth=0,
              headaxislength=0,
              label=f'Ref: {ref_norm}')
    ax.text(ref_x + ref_norm / 2, ref_y + 0.02, 'e=0.1', color='red', ha='center')

    ax.invert_xaxis()
    ax.set_title(title)
    ax.set_xlabel('RA (deg)')
    ax.set_xlim(ra_limits)
    ax.set_ylim(dec_limits)
    if i == 0:
        ax.set_ylabel('Dec (deg)')

fig.tight_layout()

Figure 2: Whisker plot for measured stars used in PSF modeling (left), PSF models with the Piff software (Jarvis et al. 2015, center), and residuals (right) in a region with a half-degree diameter around Abell 360's BCG. Coherent structures in ellipticity are visible, but the modeling captures them well for most stars, as shown by the residuals. These structures may originate from weather patterns or from the nature of the engineering data taken by LSSTComCam, which emphasized Active Optics System optimization (SITCOMTN-149).

Calculate the histograms for the ellipticity residuals.

In [31]:

fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(10, 5))

ax[0].hist((e1_star_used - e1_psf_used), bins=30, range=[-0.04, 0.04],
           linestyle='-', histtype='step', density=True, label='used')
ax[0].hist((e1_star_reserved-e1_psf_reserved), bins=30, range=[-0.04, 0.04],
           linestyle='--', histtype='step', density=True, label='reserved')
ax[0].set_xlabel(r'$\delta e_1 = $ Star e1 - PSF e1')
ax[0].legend()

ax[1].hist((e2_star_used - e2_psf_used), bins=30, range=[-0.04, 0.04],
           linestyle='-', histtype='step', density=True,
           label='used')
ax[1].hist((e2_star_reserved - e2_psf_reserved), bins=30, range=[-0.04, 0.04],
           linestyle='--', histtype='step', density=True,
           label='reserved')
ax[1].set_xlabel(r'$\delta e_2 = $ Star e2 - PSF e2')
ax[1].legend()

plt.show()

Figure 3: Ellipticity residuals (e1 on the left; e2 on the right).

4. PSF size correlation function: object table.¶

Calculate the correlation function of the size of the PSF in the i band, using treecorr and the psf_used_table object table from the previous section.

In [32]:

psf_used_table['size_'+my_band+'_PSF'] = np.sqrt((psf_used_table[my_band + '_ixxPSF']
                                                 + psf_used_table[my_band + '_iyyPSF']) / 2)

In [33]:

obj_ra = psf_used_table['coord_ra']
obj_dec = psf_used_table['coord_dec']

In [34]:

cat = treecorr.Catalog(ra=obj_ra, dec=obj_dec,
                       k=psf_used_table[f'size_{my_band}_PSF']
                       - np.mean(psf_used_table[f'size_{my_band}_PSF']),
                       ra_units='deg', dec_units='deg')

kk_config = {'max_sep': .06, 'min_sep': .0001, 'nbins': 12}

In [35]:

kk = treecorr.KKCorrelation(kk_config)
kk.process(cat)

In [36]:

xi = kk.xi
bins = kk.rnom

In [37]:

_, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4.5),
                             gridspec_kw={'wspace': .3})

ax1.set_title(f'PSF size. Band: {my_band}')
scatter_plot = ax1.scatter(obj_ra, obj_dec, c=psf_used_table[f'size_{my_band}_PSF'],
                           s=1, cmap='cividis')
ax1.set_xlabel('RA [deg]')
ax1.set_ylabel('DEC [deg]')
plt.colorbar(scatter_plot, ax=ax1, label='PSF size [pixels]')

ax2.set_title(f'PSF size correlation. Band: {my_band}')
ax2.plot(np.degrees(bins), xi*1e4, 'o-', color='darkblue')
ax2.axhline(0, linestyle='--', color='lightgrey')
ax2.set_xscale('log')
ax2.set_ylabel(r'$\xi \times 10^{4}$', labelpad=10)
ax2.set_xlabel(r'$\theta$ [degree]')

plt.show()

Figure 4: Left: spatial distribution of PSF size (pixels) for a half a dgree region aroung Abell 360's BCG, using information from object tables. Right: PSF size correlation function in the same field as the panel on the left.

5. Tangential PSF ellipticity residuals around cluster center in Abel 360.¶

An important systematic check in PSF modeling for weak lensing is to examine whether there are any significant tangential ellipticity residuals around the cluster center, using both the used stars (those included in the PSF model) and the reserved stars (those withheld for validation). Residual tangential patterns in the stellar ellipticities—after subtracting the PSF model—can mimic or contaminate the weak lensing shear signal that is expected to be tangentially aligned around massive structures like galaxy clusters. Detecting and correcting such systematics is essential to ensure that the measured shear truly reflects the underlying gravitational lensing signal rather than imperfections in the PSF modeling.

Calculate the angular separation from the BCG.

In [38]:

r_used = SkyCoord(psf_used_table['coord_ra'],
                  psf_used_table['coord_dec'], unit='deg').separation(
    SkyCoord(ra_bcg, dec_bcg, unit='deg')).arcmin
r_max = np.max(r_used)
bins = np.linspace(0, r_max, 8)

Calculate the tangential and cross components of the residual for used stars.

In [39]:

mid_u, e_t_u, err_u, count_u = compute_et(
    psf_used_table, e1_star_used, e2_star_used, e1_psf_used, e2_psf_used,
    ra_bcg, dec_bcg, bins
)

Calculate the tangential and cross components of the residual for reserved stars.

In [40]:

mid_r, e_t_r, err_r, count_r = compute_et(
    psf_reserved_table, e1_star_reserved, e2_star_reserved,
    e1_psf_reserved, e2_psf_reserved,
    ra_bcg, dec_bcg, bins
)

Plot the profiles.

In [41]:

plt.figure(figsize=(6, 4))
plt.errorbar(mid_u, e_t_u, err_u, fmt='o', label='Used stars')
plt.errorbar(mid_r, e_t_r, err_r, fmt='s', label='Reserved stars')

plt.axhline(0, color='k', linestyle=':')
plt.xlabel('Separation [arcmin]', fontsize=14)
plt.ylabel(r'Residuals: $\langle e_t \rangle$', fontsize=14)
plt.ylim([-0.006, 0.006])
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.legend(fontsize=12)
plt.tight_layout()
plt.show()

Figure 5: Tangential ellipticity PSF residuals around Abell 360's center.

In [ ]: