Skip to content

Photometric Surveys Summary Statistics (Using Tracers)

AngularTwoPoint

Compute angular power spectra \(C_\ell\), pseudo-\(C_\ell\) and COSEBIs from two tracers.

Location: cloelib/summary_statistics/angular_two_point.py

What it does:

  • Takes two Tracer objects
  • Integrates over redshift using the Limber approximation with advanced vectorisation across tomographic bins

Example:

from cloelib.cosmology.camb_cosmology import CAMBBackground, CAMBLinearPerturbations
from cloelib.observables.photo import ShearTracer, PositionsTracer
from cloelib.summary_statistics.angular_two_point import AngularTwoPoint
import numpy as np

# Set up cosmology
bg = CAMBBackground(H0=67.5, Omega_b0=0.0492, ...)

# Create tracers
z = np.linspace(0.01, 3.0, 100)
pert = CAMBLinearPerturbations(background=bg, redshifts=z)
dndz = np.exp(-((z - 0.7) / 0.3)**2)
dndz = dndz / np.trapz(dndz, z)

tracer1 = ShearTracer(perturbations=pert, dndz=dndz[np.newaxis, :], z=z, nuisance_params={...})
tracer2 = ShearTracer(perturbations=pert, dndz=dndz[np.newaxis, :], z=z, nuisance_params={...})

# Compute angular power spectrum
two_point = AngularTwoPoint(tracer1=tracer1, tracer2=tracer2)

ells = np.logspace(1, 3, 20)  # ℓ from 10 to 1000
C_ell = two_point.get_Cl(ells=ells, nl=0, ks=pert.k)

print(f"C_ℓ at ℓ=100: {C_ell[0, 0, 10]:.2e}")

Cross-Correlations:

You can correlate different tracers to calculate different statistics:

# Shear-shear (cosmic shear)
shear_tracer = ShearTracer(...)
C_shear_shear = AngularTwoPoint(shear_tracer, shear_tracer).get_Cl(ells, 0, pert.k)

# Position-position (galaxy clustering)
pos_tracer = PositionsTracer(...)
C_gg = AngularTwoPoint(pos_tracer, pos_tracer).get_Cl(ells, 0, pert.k)

# Shear-position (galaxy-galaxy lensing)
C_g_shear = AngularTwoPoint(pos_tracer, shear_tracer).get_Cl(ells, 0, pert.k)

Tomographic Bins:

# Multiple redshift bins
dndz_bins = np.array([
    np.exp(-((z - 0.5) / 0.2)**2),
    np.exp(-((z - 1.0) / 0.3)**2),
    np.exp(-((z - 1.5) / 0.4)**2),
])

# Normalise! cloelib will always expect the n(z) normalised
dndz_bins = dndz_bins / np.trapz(dndz_bins, z, axis=1)[:, np.newaxis]

tracer = ShearTracer(perturbations=pert, dndz=dndz_bins, z=z, ...)

# Auto and cross-correlations
two_point = AngularTwoPoint(tracer, tracer)
C_ell = two_point.get_Cl(ells, nl=0, ks=pert.k)

COSEBIs: Complete Orthogonal Sets of E/B-Integrals are specialized statistics for cosmic shear. They are available both as methods of AngularTwoPoint (which injects a software provenance tag automatically) and as standalone module-level functions that can be used without instantiating the class, useful when \(C_\ell\) or \(\xi(\pm)\) are already computed, even from real or simulated measurements.

Standalone functions (no tracers required):

from cloelib.summary_statistics.angular_two_point import (
    get_cosebis_from_cl,
    get_cosebis_from_2pcf,
)

# From angular power spectra
cosebis = get_cosebis_from_cl(cells, ells, w_ell, ns)

# From two-point correlation functions
cosebis = get_cosebis_from_2pcf(twopcf, theta, T_plus, T_minus, ns)

Class method (get_cosebis — computes \(C_\ell\) internally, software provenance tag included automatically):

two_point = AngularTwoPoint(tracer1, tracer2)
cosebis = two_point.get_cosebis(ells, nl, ks, w_ell, ns)

Note

Only get_cosebis is available as a class method. get_cosebis_from_2pcf is standalone only — use it directly when \(\xi_\pm(\theta)\) are already available.

Both interfaces require optional dependencies (pylevin, mpmath). The COSEBIs kernels (w_ell, T_plus, T_minus) must be precomputed using helpers from cloelib.auxiliary.cosebi_helpers.

Output format and units: The output matches the Euclid LE3 COSEBI products read by euclidlib.le3.twopcf_wl.cosebis. Each COSEBI has an array of shape (2, 2, n_modes), with EE in array[0, 0] and BB in array[1, 1]. The EB/BE entries are zero in the theory output. thmin/thmax are in arcmin.

The inputs are in radians: theta in get_cosebis_from_2pcf, and the thetagrid passed to get_W_ell, whose metadata["THMIN"]/["THMAX"] are therefore in radians too. For example, for kernels matching a 1–400 arcmin product:

import numpy as np
from cloelib.auxiliary.cosebi_helpers import get_W_ell

thetagrid = np.radians(np.geomspace(1.0, 400.0, 500) / 60)  # arcmin -> rad
w_ell = get_W_ell(thetagrid, Nmax=20, ells=ells, N_thread=4)

Integration grids: The integrals over \(\ell\) (and over \(\theta\) for get_cosebis_from_2pcf) use the actual grid spacing, so ells and theta can be linearly or logarithmically spaced. Simpson's rule is used on linear and log grids, and the trapezoidal rule on any other increasing grid. The ells passed to get_cosebis_from_cl/get_cosebis must be the grid on which w_ell was evaluated, and it must be fine enough to resolve the oscillations of the kernels, whose period in \(\ell\) is roughly \(2\pi/\theta_{\max}\). Check convergence by comparing against a denser grid.

AngularCorrelationFunction

Protocol to set up how to compute real-space angular correlation functions \(\xi(\theta)\).

Location: cloelib/summary_statistics/angular_correlation_function.py

AngularCorrelationFunctionWigner

Alternative implementation using Wigner 3j symbols for computing real-space angular correlation functions \(\xi(\theta)\) from harmonic angular power spectra.

Location: cloelib/summary_statistics/angular_correlation_function_wigner.py

Next Steps

The pipeline is now complete. From here:

Or return to any component: