2.2 — EDFS-eRASS cluster 1 (mass fit + aperture-mass map)

Authors: Prakruth Adari, Xiangchong Li LSST Science Pipelines: Weekly 2025_49 Kernel: Python-dp1 (uses clmm, pyccl, emcee).

Purpose

Take the EDFS source catalog produced by 2_1_edfs_photoz_estimate.ipynb, build a clmm.GalaxyCluster object around the EDFS-eRASS cluster at z = 0.6922, fit M_500c (NFW, fixed concentration) by both curve_fit and MCMC, and produce a Schirmer aperture-mass S/N map.

Steps

1. Load EDFS AnaCal table + cached photo-z point estimates and per-source PDFs (data_edfs_redshift{,_pdfs}.{hdf5,fits}); add zbest = zmode.
2. Define the cluster (ra, dec, z) = (59.487317, -49.000349, 0.6922) (EDFS_eRASS; commented alternative for EDFS_DES).
3. Compute per-source angular separation and projected distance (Da(z_cl) * angle).
4. Mask with (u<27.5, g<27.5, r<24.5, i<27.5, z<27.5, y<27.5), FPFS trace

   0.1, and zbest > 0.78 → background source sample.

5. Combine the 4-band FPFS moments → (e1, e2) and selection response R; rotate to tangential/cross relative to the cluster centre.
6. Bin shear in log-radius (0.2 → 9 Mpc), bootstrap CIs, plot tangential and cross profiles, report √χ² and σ.
7. Build the source-sample n(z) by response-weighted stacking of the per-source PDFs; compute β_s and β_s² from that n(z).
8. Wrap the sample in a clmm.GalaxyCluster, fit M_500c with fixed concentration (NFW, c=2.26) via curve_fit, then refine with a small emcee chain (32 walkers × 100 steps).
9. Histogram the posterior in log10(M_500c) and mark eRASS allowed range (low/high logM).
10. Build a flat-tangent WCS centred on the BCG and evaluate the Schirmer aperture-mass filter on a 161×161 grid; plot E/B-mode S/N histograms and 2D maps.

Inputs

• rail_data/data_edfs_with_mag.hdf5
• rail_data/data_edfs_redshift.hdf5 (cached point estimates from 2.1)
• rail_data/data_edfs_redshift_pdfs.fits (cached PDFs from 2.1)
• catalogs/anacal_catalog_ecdfs.fits (used as sanity reference only)

Outputs

Inline figures: (RA, Dec) + color-magnitude scatter, tangential / cross shear profile, source-vs-full n(z), mass posterior histogram, and E/B-mode aperture-mass S/N maps. No files written.

import os

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate
import scipy.stats as stats

import tables_io
from astropy.table import Table
from astropy.io import fits
from matplotlib import cm
import emcee

import clmm
from clmm import GCData
from clmm.support.sampler import fitters
cosmo = clmm.Cosmology(H0=70.0, Omega_dm0=0.3 - 0.045, Omega_b0=0.045, Omega_k0=0.0)

plt.rcParams.update({
    "text.usetex": False,
    "font.family": "sans-serif",
    "font.sans-serif": ["Times"],
    "figure.labelsize": 14,
    "mathtext.fontset": 'stix',
    "font.family": "STIXGeneral",
    "font.size": 14,
    "figure.titlesize": 14
})

%config InlineBackend.figure_format='retina'

arcsec = 1 / 60**2
rng = np.random.default_rng()
omega_m = .31
omega_de= .69
omega_r = 0
H0 = 70 # km/s/Mpc
c = 3e5 # km/s

Hz = lambda z : c/(H0  * np.sqrt((omega_de + omega_m * (1+z)**3 + omega_r * (1+z)**4)))
chi_dl = lambda z, z0=0 : integrate.quad(Hz, z0, z)[0]
Da = lambda z : chi_dl(z)/(1+z)
Daz1z2 = lambda z1, z2 : 1/(1+z2)*(chi_dl(z2) - chi_dl(z1))
# beta_r = lambda zl, zs : integrate.quad(Hz, zl, zs)[0]/integrate.quad(Hz, 0, zs)[0]
beta_s = lambda zl, zs : 0 if zs < zl else Daz1z2(zl, zs)/Da(zs) * Da(1e3)/Daz1z2(zl, 1e3)

pscratch= "/gpfs02/astro/workarea/xli6/"
Dir = os.path.join(
    pscratch,
    "work/DP1/"
)

edfs_anacal = tables_io.read(
    f'{Dir}/rail_data/data_edfs_with_mag.hdf5',
    tType = "astropyTable"
)
ecdfs_anacal = Table.read(f'{Dir}/catalogs/anacal_catalog_ecdfs.fits')
edfs_redshifts = tables_io.read(
    f'{Dir}/rail_data/data_edfs_redshift.hdf5',
    tType = "astropyTable"
)
edfs_pzs = tables_io.read(
    f'{Dir}/rail_data/data_edfs_redshift_pdfs.fits',
    tType = "astropyTable"
)
with fits.open(f'{Dir}/rail_data/data_edfs_redshift_pdfs.fits') as hdul:
    edfs_pzs = hdul[0].data
    edfs_pzs_meta = hdul[0].header

edfs_anacal['zbest'] = edfs_redshifts['zmode']
ecdfs_anacal['index'] = np.arange(len(ecdfs_anacal))
edfs_anacal['index'] = np.arange(len(edfs_anacal))

from utils import anacal_get_tang_cross, calibrate_shapes
from utils import build_flat_wcs_grid, compute_mass_map, schirmer_filter

mag_cuts = [27.5, 27.5, 24.5, 27.5, 27.5, 27.5]

def get_mask(df, z_cut, SG_cut=0.1, mag_cuts=mag_cuts):
    msk = (
        (df["u_mag_gauss2"] < mag_cuts[0]) &
        (df["g_mag_gauss2"] < mag_cuts[1]) &
        (df["r_mag_gauss2"] < mag_cuts[2]) &
        (df["i_mag_gauss2"] < mag_cuts[3]) &
        (df["z_mag_gauss2"] < mag_cuts[4]) &
        (df["y_mag_gauss2"] < mag_cuts[5]) &
        ((df["i_fpfs1_m00"] + df["i_fpfs1_m20"]) / df["i_fpfs1_m00"] > SG_cut) &
        (df["zbest"] > z_cut)
    )
    return df[msk], msk


def add_cluster_info(df, cluster_info):
    source_phi = np.arctan2(
        df['dec'] - cluster_info[1],
        (cluster_info[0] - df['ra']) * np.cos(np.deg2rad(cluster_info[1])),
    )
    ang_dist = np.sqrt(
        ((df['ra'] - cluster_info[0]) * np.cos(np.deg2rad(cluster_info[1]))) ** 2
        + (df['dec'] - cluster_info[1]) ** 2
    )
    sky_distance = cosmo.eval_da(cluster_info[2]) * ang_dist * (np.pi / 180)
    df['source_phi'] = source_phi
    df['ang_dist'] = ang_dist
    df['sky_dist'] = sky_distance
    return source_phi, ang_dist, sky_distance

# RA, DEC, z
cluster_info = (59.487317, -49.000349, 0.6922) # EDFS_ERASS
#cluster_info = (59.482185, -48.998922, 0.6916004419) # EDFS_DES
_ = add_cluster_info(edfs_anacal, cluster_info)

mag_cuts = np.array([27.5, 27.5, 24.5, 27.5, 27.5, 27.5])
masked_df, _ = get_mask(edfs_anacal, 0.78, mag_cuts=mag_cuts)
e1, e2, res = calibrate_shapes(masked_df)
masked_df['response'] = res

trial_shear = e1 + 1.j*e2
cl_shear = trial_shear * -1*np.exp(-2j*masked_df['source_phi'])

bins_mpc = clmm.make_bins(0.2, 9, nbins=7, method='evenlog10width')

bin_mids = 0.5 * (bins_mpc[1:] + bins_mpc[:-1])
shear_cl = anacal_get_tang_cross(
    cl_shear, masked_df['sky_dist'],
    bins_mpc, res, ci_level=.68,
    r_correction_offset=0.01,
)

chi2_t = np.sum(shear_cl[0] ** 2.0 / ((shear_cl[2][:,0] - shear_cl[2][:,1]) / 2.0) ** 2.0)
chi2_x = np.sum(shear_cl[1] ** 2.0 / ((shear_cl[3][:,0] - shear_cl[3][:,1]) / 2.0) ** 2.0)
print(np.sqrt(chi2_t))
print(np.sqrt(chi2_x))

sigma_t = stats.norm.isf(stats.chi2.sf(chi2_t, len(bins_mpc)))
sigma_x = stats.norm.isf(stats.chi2.sf(chi2_x, len(bins_mpc)))
print(sigma_t)
print(sigma_x)

4.550335961931738
2.230952678660221
2.4101898463000206
-0.7063526724970876

fig, ax = plt.subplots(ncols=2, figsize=(10, 5), layout='constrained')

ax[0].scatter(masked_df['ra'], masked_df['dec'],
              s=0.1, c=masked_df['ang_dist'])
ax[0].set_xlabel("RA")
ax[0].set_ylabel("Dec")

cmap = cm.coolwarm
ax[1].plot(bin_mids, shear_cl[0], '.', label=f'tangential - {sigma_t:0.2f}σ', color=cmap(.05))
ax[1].vlines(bin_mids, shear_cl[2][:,0], shear_cl[2][:,1], color=cmap(0.05))
ax[1].plot(1.05*bin_mids, shear_cl[1], '.', label=f'cross - {sigma_x:0.2f}σ', color=cmap(.95))
ax[1].vlines(1.05*bin_mids, shear_cl[3][:,0], shear_cl[3][:,1], color=cmap(0.95))
ax[1].axhline(0, ls='--', color='k', alpha=.5)
ax[1].semilogx()
ax[1].set_xlim([0.2, 8.5])
ax[1].set_ylabel("Reduced shear")
ax[1].set_xlabel("R (Mpc)")
ax[1].legend(frameon=False, loc='upper right')

z = cluster_info[2]
fig.suptitle(f"EDFS-eRASS Cluster at {z=:0.2f}")

Text(0.5, 0.98, 'EDFS-eRASS Cluster at z=0.69')

N(z) Conversion

test_zs = np.arange(0, 5.01, 0.01)
response_stacked_nz = np.zeros(len(test_zs))

nz_ndxs = masked_df['index']
response_stacked_nz = np.dot(masked_df['response'], edfs_pzs[nz_ndxs])
nz_normalization = integrate.trapezoid(response_stacked_nz, test_zs)
normalized_nz = response_stacked_nz/nz_normalization

full_nz = np.sum(edfs_pzs, axis=0)
full_norm = integrate.trapezoid(full_nz, test_zs)
full_norm_nz = full_nz/full_norm

cmap = cm.coolwarm
plt.plot(test_zs, normalized_nz, '-', label='Source Sample', color=cmap(0.95))
plt.plot(test_zs, full_norm_nz, '-', label='Full Sample', color=cmap(0.05))
plt.xlabel("Redshift")
plt.legend()
plt.ylabel("Normalized $N(z)$")
plt.title("EDFS Redshift Distribution")
plt.axvline(cluster_info[2], ls='--', color='k', alpha=0.3)

<matplotlib.lines.Line2D at 0x7fdeb6ba4590>

def get_beta_betasqr_pdf(nz, beta_s, beta_sqr, zs):
    normalization = integrate.trapezoid(nz, zs)
    avg_beta_s = integrate.trapezoid(beta_s * nz, zs) / normalization
    avg_beta_s_sqr = integrate.trapezoid(beta_sqr * nz, zs) / normalization
    return avg_beta_s, avg_beta_s_sqr

test_xs = np.linspace(0, 5, 501)

cluster_info[2]

0.6922

beta_s_vals = np.array([beta_s(cluster_info[2], xx) for xx in test_xs])
beta_sqr_vals = beta_s_vals**2

cluster_beta_s, cluster_beta_s_sqr = get_beta_betasqr_pdf(normalized_nz, beta_s_vals, beta_sqr_vals, test_xs)

cluster_beta_s

np.float64(0.3446779119901067)

cluster_beta_s_sqr

np.float64(0.1582670753673673)

CLMM Mass Fit

def fit_mass(predict_function, fit_data, **kwargs):
    popt, pcov = fitters["curve_fit"](
        predict_function,
        fit_data["radius"],
        fit_data["gt"],
        fit_data["gt_err"],
        **kwargs,
    )
    return {
        "vals": popt,
        "err": np.sqrt(np.diag(pcov)),
        "cov": pcov,
    }


def fixed_concentration_fit(profile, beta, betasqr, cdelta, halo_profile="nfw"):
    fn = lambda radius, logm: clmm.compute_reduced_tangential_shear(
                r_proj=radius,  # Radial component of the profile
                mdelta=10**logm,  # Mass of the cluster [M_sun]
                cdelta=cdelta,  # Concentration of the cluster
                z_cluster=cluster_info[2],  # Redshift of the cluster
                z_src=(beta, betasqr),  # tuple of (bs_mean, bs2_mean)
                z_src_info="beta",
                approx="order2",
                cosmo=cosmo,
                delta_mdef=500,
                massdef="critical",
                halo_profile_model=halo_profile
                )

    _bounds = [10.0, 17.0]
    mass_fit = fit_mass(fn, profile.profile, bounds=_bounds)
    return fn, mass_fit
    
def mcmc_fit_1d(profile, fn, init, nwalkers=32, nsteps=100):
    fit_data = profile.profile
    def log_prob(pars):

        if pars[0] < 10:
            return -np.inf
        if pars[0] > 16:
            return -np.inf

        model = fn(fit_data["radius"], *pars)
        diff = model - fit_data["gt"]
        sig = fit_data["gt_err"]
        lnlike = -0.5 * diff**2 / sig**2

        return lnlike.sum()

    sampler = emcee.EnsembleSampler(nwalkers, 1, log_prob)

    p0 = np.random.normal(init["vals"], init["err"] * 0.3, (nwalkers, 1))
    p0[p0 < 0] *= -1
    sampler.run_mcmc(p0, nsteps, progress=True)
    return {
        "mcmc_chain": sampler.get_chain().flatten(),
        "mcmc_loglike": sampler.get_log_prob().flatten(),
    }

global_R_fit = float(np.mean(res))
print(f"global R (mass fit) = {global_R_fit:.4f}")

galcat = GCData()
galcat['ra'] = masked_df['ra']
galcat['dec'] = masked_df['dec']
galcat['e1'] = e1 / global_R_fit   # AnaCal response correction
galcat['e2'] = e2 / global_R_fit
galcat['id'] = np.arange(len(masked_df))
galcat['z'] = masked_df['zbest']

global R (mass fit) = 0.2458

/gpfs/mnt/gpfs02/astro/workarea/xli6/superonionSDCC/miniconda3/envs/analysis/lib/python3.12/site-packages/clmm/gcdata.py:30: UserWarning: coordinate_system not set, defaulting to 'euclidean'.
  warnings.warn("coordinate_system not set, defaulting to 'euclidean'.")

ra_bcg = cluster_info[0]
dec_bcg = cluster_info[1]

cluster_id = "EDFS z0.7"
# gc_object1 = clmm.GalaxyCluster(cluster_id, ra_bcg, dec_bcg, 0.22, galcat, coordinate_system='euclidean')
gc_object1 = clmm.GalaxyCluster(cluster_id, ra_bcg, dec_bcg, cluster_info[2], galcat)

gc_object1.compute_tangential_and_cross_components(add=True);

# bins_mpc = np.array([0.3, 0.7       , 1.06495522, 1.62018517, 2.46489237, 3.75    , 5, 6 ])
bins_mpc = clmm.make_bins(0.2,9,nbins=7, method='evenlog10width')
bins_mpc = clmm.make_bins(0.3,7,nbins=5, method='evenlog10width')


gc_object1.make_radial_profile(bins=bins_mpc, bin_units='Mpc', add=True, cosmo=cosmo, 
                               overwrite=True, use_weights=False, gal_ids_in_bins=False);

profile = gc_object1.profile

low_logM = np.log10(27.67 * 1e13)
high_logM = np.log10(38.84 * 1e13)

fn, init = fixed_concentration_fit(
    gc_object1, 
    cluster_beta_s, 
    cluster_beta_s_sqr,
    2.26, "nfw",
)

mcmc_results = mcmc_fit_1d(gc_object1, fn, init, nwalkers=32, nsteps=100)
np_mcmc = np.zeros((32*100, 2))
np_mcmc[:, 0] = mcmc_results['mcmc_chain']
np_mcmc[:, 1] = mcmc_results['mcmc_loglike']

  0%|          | 0/100 [00:00<?, ?it/s]

  2%|▏         | 2/100 [00:00<00:09, 10.49it/s]

  4%|▍         | 4/100 [00:00<00:08, 10.77it/s]

  6%|▌         | 6/100 [00:00<00:08, 10.78it/s]

  8%|▊         | 8/100 [00:00<00:08, 10.86it/s]

 10%|█         | 10/100 [00:00<00:08, 10.95it/s]

 12%|█▏        | 12/100 [00:01<00:08, 10.99it/s]

 14%|█▍        | 14/100 [00:01<00:07, 11.05it/s]

 16%|█▌        | 16/100 [00:01<00:07, 11.15it/s]

 18%|█▊        | 18/100 [00:01<00:07, 11.06it/s]

 20%|██        | 20/100 [00:01<00:07, 10.92it/s]

 22%|██▏       | 22/100 [00:02<00:07, 10.92it/s]

 24%|██▍       | 24/100 [00:02<00:06, 10.99it/s]

 26%|██▌       | 26/100 [00:02<00:06, 10.91it/s]

 28%|██▊       | 28/100 [00:02<00:06, 11.00it/s]

 30%|███       | 30/100 [00:02<00:06, 11.07it/s]

 32%|███▏      | 32/100 [00:02<00:06, 11.04it/s]

 34%|███▍      | 34/100 [00:03<00:05, 11.26it/s]

 36%|███▌      | 36/100 [00:03<00:05, 11.43it/s]

 38%|███▊      | 38/100 [00:03<00:05, 11.46it/s]

 40%|████      | 40/100 [00:03<00:05, 11.51it/s]

 42%|████▏     | 42/100 [00:03<00:05, 11.40it/s]

 44%|████▍     | 44/100 [00:03<00:04, 11.32it/s]

 46%|████▌     | 46/100 [00:04<00:04, 11.25it/s]

 48%|████▊     | 48/100 [00:04<00:04, 10.97it/s]

 50%|█████     | 50/100 [00:04<00:04, 10.95it/s]

 52%|█████▏    | 52/100 [00:04<00:04, 10.64it/s]

 54%|█████▍    | 54/100 [00:04<00:04, 10.53it/s]

 56%|█████▌    | 56/100 [00:05<00:04, 10.51it/s]

 58%|█████▊    | 58/100 [00:05<00:04, 10.39it/s]

 60%|██████    | 60/100 [00:05<00:03, 10.32it/s]

 62%|██████▏   | 62/100 [00:05<00:03, 10.43it/s]

 64%|██████▍   | 64/100 [00:05<00:03, 10.48it/s]

 66%|██████▌   | 66/100 [00:06<00:03, 10.46it/s]

 68%|██████▊   | 68/100 [00:06<00:03, 10.43it/s]

 70%|███████   | 70/100 [00:06<00:02, 10.51it/s]

 72%|███████▏  | 72/100 [00:06<00:02, 10.77it/s]

 74%|███████▍  | 74/100 [00:06<00:02, 10.79it/s]

 76%|███████▌  | 76/100 [00:06<00:02, 10.82it/s]

 78%|███████▊  | 78/100 [00:07<00:02, 10.62it/s]

 80%|████████  | 80/100 [00:07<00:01, 10.57it/s]

 82%|████████▏ | 82/100 [00:07<00:01, 10.58it/s]

 84%|████████▍ | 84/100 [00:07<00:01, 10.63it/s]

 86%|████████▌ | 86/100 [00:07<00:01, 10.70it/s]

 88%|████████▊ | 88/100 [00:08<00:01, 10.62it/s]

 90%|█████████ | 90/100 [00:08<00:00, 10.57it/s]

 92%|█████████▏| 92/100 [00:08<00:00, 10.56it/s]

 94%|█████████▍| 94/100 [00:08<00:00, 10.44it/s]

 96%|█████████▌| 96/100 [00:08<00:00, 10.43it/s]

 98%|█████████▊| 98/100 [00:09<00:00, 10.42it/s]

100%|██████████| 100/100 [00:09<00:00, 10.38it/s]

100%|██████████| 100/100 [00:09<00:00, 10.77it/s]

fig = plt.figure(figsize=(7,3))

cmap = cm.coolwarm


mass_bins = np.arange(13.1, 15.5, 0.05)

chain = np_mcmc[:,0]
ll = np_mcmc[:,1]
fin_filt = np.isfinite(ll)
chain = chain[fin_filt]
ll = ll[fin_filt]

plt.hist(chain, bins=mass_bins, weights=np.exp(ll), density=True, linewidth=2,
         histtype='step', color=cmap(0.05), linestyle='-')
plt.hist(chain, bins=mass_bins, weights=np.exp(ll), density=True, linewidth=2,
         histtype='stepfilled', color=cmap(0.05), linestyle='-', alpha=0.2)

plt.xlim(13, 16)
# plt.legend(frameon=False)
plt.xlabel(r"$\log{M_{500c}/M_\odot}$")
# plt.axvline(logM, ls='--', color='k', alpha=0.5)
plt.axvline(low_logM, ls='--', color='k', alpha=0.5)
plt.axvline(high_logM, ls='--', color='k', alpha=0.5)
plt.ylabel("Normalized Counts")

Text(0, 0.5, 'Normalized Counts')

Schirmer aperture-mass S/N map

Build a flat-tangent WCS centred on the BCG, project sources into degree-scale (x, y), then evaluate the Schirmer-filter aperture mass on a 161×161 grid. The S/N is taken w.r.t. the variance map.

N, rr = 161, 0.16
x, y, x_grid, y_grid = build_flat_wcs_grid(
    ra_bcg, dec_bcg, masked_df['ra'], masked_df['dec'],
    n=N, half_extent_deg=rr,
)
global_R = float(np.mean(res))
print(f"global R = {global_R:.4f}")
g1 = e1 / global_R
g2 = e2 / global_R
weights = np.ones(len(x))

e_ap, b_ap, v_ap = compute_mass_map(
    x_grid, y_grid, x, y, g1, g2, weights,
    schirmer_filter,
    filter_kwargs={'aperture_size': 0.15, 'a': 6, 'b': 150},
)

global R = 0.2458

test_xs = np.linspace(-5, 5, num=1001)
normal_dist = stats.norm(0, 1).pdf(test_xs)
sn_v = np.sqrt(np.mean(v_ap))
plt.hist((b_ap / sn_v).flatten(), bins=51, label='B Mode',
         histtype='step', density=True)
plt.hist((e_ap / sn_v).flatten(), bins=51, label='E Mode',
         histtype='step', density=True)
plt.plot(test_xs, normal_dist, '--')
plt.legend()

<matplotlib.legend.Legend at 0x7fdee1827c80>

fig, axs = plt.subplots(ncols=2, figsize=(10, 6))
extent = (rr, -rr, -rr, rr)
sn_v = np.sqrt(np.mean(v_ap))

axs[0].set_title('E-Mode SN')
MapE = axs[0].imshow(e_ap / sn_v, vmin=-2, vmax=3.5, origin='lower', extent=extent)
y_p, x_p = np.where(e_ap / sn_v == np.max(e_ap / sn_v))
print((e_ap / sn_v)[y_p, x_p])
x_vals = np.linspace(extent[0], extent[1], N)
y_vals = np.linspace(extent[2], extent[3], N)
axs[0].scatter(x_vals[x_p], y_vals[y_p], s=5.0, color='red')

axs[1].set_title('B-Mode SN')
axs[1].imshow(b_ap / sn_v, vmin=-2, vmax=3.5, origin='lower', extent=extent)

fig.colorbar(MapE, ax=axs, fraction=0.025)
fig.supxlabel('RA',  y=0.14)
fig.supylabel('DEC', x=0.05)

[4.47059291]

Text(0.05, 0.5, 'DEC')