import io
import urllib.request
from datetime import datetime
from math import ceil, log10
import matplotlib.image as mpimg
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import gridspec
from matplotlib.backends.backend_pdf import PdfPages
from scipy.interpolate import interp1d
import lephare as lp
# Compatibility shim
trapezoid = np.trapezoid if hasattr(np, "trapezoid") else np.trapz
__all__ = ["PlotUtils"]
[docs]
class PlotUtils:
"""
Utility class for generating diagnostic plots and statistical summaries
from LePHARE output tables.
This class provides a suite of methods for analyzing photometric redshift
and physical parameters catalogs generated by LePHARE, including comparisons between photometric
and spectroscopic redshifts, redshift error distributions, chi-squared
diagnostics, magnitude and color relations, and other common validation plots.
Examples
--------
Typical usage example::
from astropy.table import Table
import lp
t = Table.read("outputphotoz.fits")
utils = lp.PlotUtils(
t,
sel_filt=3,
pos_filt=[0, 1, 2, 4, 5, 5],
range_z=[0, 0.5, 1, 1.5, 3],
range_mag=[19, 20.5, 21.5, 22.5, 25]
)
utils.zml_zs()
utils.zp_zs()
utils.zml_zp()
utils.dist_z()
utils.dist_chi2()
utils.dist_filt()
utils.dist_model()
utils.dist_ebv()
utils.secondpeak()
utils.bzk()
utils.cumulative68()
utils.check_error()
utils.errormag()
utils.errorz()
utils.pit_qq()
utils.absmag_z()
utils.rf_color()
utils.william()
utils.dist_mass()
utils.dist_sfr()
utils.dist_ssfr()
utils.mass_med_best()
utils.sfr_med_best()
utils.mass_sfr()
utils.mass_z()
utils.sfr_z()
utils.lnuv_sfr()
utils.absmagu_sfr()
utils.absmagk_mass()
utils.masstolight_z()
utils.save_all_plots_pdf()
Parameters
----------
t : astropy.table.Table
Table containing the output from a LePHARE run. Must include columns such as
``Z_BEST``, ``Z_MED``, ``MAG_OBS()``, ``MAG_ABS()``, ``CHI_BEST``, etc.
sel_filt : int, optional
Index (starting at 0) of the filter used to select objects by observed magnitude.
Default is 0.
pos_filt : list of int, optional
Indices (starting at 0) of the filters corresponding to u, g, r, z, J, and Ks bands.
Used to generate color-based plots. Default is ``[0, 0, 0, 0, 0, 0]``.
range_z : list of float, optional
Redshift bin edges to use in diagnostics. If not provided, quartiles of the
spectroscopic redshift distribution (0, 0.25, 0.5, 0.75, 1 quantiles) are used.
range_mag : list of float, optional
Observed magnitude bin edges to use in diagnostics. If not provided, quartiles of
the selected-band magnitude distribution (0, 0.25, 0.5, 0.75, 1 quantiles) are used.
Attributes
----------
Id : ndarray
Object identifiers.
zp, zml, zs : ndarray
Best-fit, median, and spectroscopic redshifts, respectively.
chi, chi2 : ndarray
Chi-squared values for best-fit and secondary models.
mag, magu, magb, magr, magz, magj, magk : ndarray
Observed magnitudes in various filters.
mabsu, mabsb, mabsr, mabsz, mabsj, mabsk : ndarray
Absolute magnitudes in various filters.
range_z, range_mag : ndarray
Bin edges used for redshift and magnitude-based diagnostics.
cond, condgal, condstar, condspec : ndarray of bool
Boolean selection masks for general, galaxy, star, and spectroscopic subsets.
z_min, z_max, mag_min, mag_max : float
Minimum and maximum values for redshift and magnitude selection.
nbRowM, nbColM, nbRowZ, nbColZ : int
Number of subplot rows and columns used in magnitude and redshift binnings.
Notes
-----
- The filter numbering starts at 0.
- The table must contain all expected columns from a LePHARE output table.
- The class methods produce a variety of plots to assess photometric redshift
performance and SED fitting reliability.
"""
def __init__(self, t, sel_filt=0, pos_filt=None, range_z=None, range_mag=None, pdf_col=None):
if pos_filt is None:
pos_filt = [0, 0, 0, 0, 0, 0]
if range_z is None:
range_z = [-1]
if range_mag is None:
range_mag = [-1]
if pdf_col is None:
pdf_col = "PDF_BAY_ZG()"
# Get the number of columns in MAG_OBS
mag_obs = t["MAG_OBS()"]
n_cols = mag_obs.shape[1] # Number of columns
# Check if sel_filt is within valid range
if sel_filt >= n_cols or sel_filt < 0:
print("Warning: sel_filt out of bounds. Setting to 0.")
self.sel_filt = 0 # Valeur par défaut
# Check if each element in pos_filt is within valid range
for i, filt in enumerate(pos_filt):
if filt >= n_cols or filt < 0:
print(f"Warning: pos_filt[{i}] out of bounds. Setting to 0.")
pos_filt[i] = 0 # Valeur par défaut
# Number of the filter start at 0
[docs]
self.sel_filt = sel_filt # filter for the selection in mag
[docs]
self.uFilt = pos_filt[0]
[docs]
self.bFilt = pos_filt[1]
[docs]
self.rFilt = pos_filt[2]
[docs]
self.zFilt = pos_filt[3]
[docs]
self.jFilt = pos_filt[4]
[docs]
self.kFilt = pos_filt[5]
# Check the position of the filter is compatible with the maximum number of filter
# Read the input file
[docs]
self.zl68 = t["Z_BEST68_LOW"]
[docs]
self.zu68 = t["Z_BEST68_HIGH"]
[docs]
self.zmll68 = t["Z_MED68_LOW"]
[docs]
self.zmlu68 = t["Z_MED68_HIGH"]
[docs]
self.chi = t["CHI_BEST"]
[docs]
self.mod = t["MOD_BEST"]
[docs]
self.law = t["EXTLAW_BEST"]
[docs]
self.ebv = t["EBV_BEST"]
[docs]
self.chi2 = t["CHI_SEC"]
[docs]
self.mod2 = t["MOD_SEC"]
[docs]
self.ebv2 = t["EBV_SEC"]
[docs]
self.chiq = t["CHI_QSO"]
[docs]
self.modq = t["MOD_QSO"]
[docs]
self.mods = t["MOD_STAR"]
[docs]
self.chis = t["CHI_STAR"]
[docs]
self.mag = t["MAG_OBS()"][:, self.sel_filt]
[docs]
self.magu = t["MAG_OBS()"][:, self.uFilt]
[docs]
self.magb = t["MAG_OBS()"][:, self.bFilt]
[docs]
self.magr = t["MAG_OBS()"][:, self.rFilt]
[docs]
self.magz = t["MAG_OBS()"][:, self.zFilt]
[docs]
self.magj = t["MAG_OBS()"][:, self.jFilt]
[docs]
self.magk = t["MAG_OBS()"][:, self.kFilt]
[docs]
self.mabsu = t["MAG_ABS()"][:, self.uFilt]
[docs]
self.mabsb = t["MAG_ABS()"][:, self.bFilt]
[docs]
self.mabsr = t["MAG_ABS()"][:, self.rFilt]
[docs]
self.mabsz = t["MAG_ABS()"][:, self.zFilt]
[docs]
self.mabsj = t["MAG_ABS()"][:, self.jFilt]
[docs]
self.mabsk = t["MAG_ABS()"][:, self.kFilt]
[docs]
self.scale = t["SCALE_BEST"]
[docs]
self.nbFilt = t["NBAND_USED"]
[docs]
self.context = t["CONTEXT"]
[docs]
self.ageb = t["AGE_BEST"]
[docs]
self.agel = t["AGE_INF"]
[docs]
self.agem = t["AGE_MED"]
[docs]
self.ages = t["AGE_SUP"]
[docs]
self.ldustb = t["LDUST_BEST"]
[docs]
self.ldustl = t["LDUST_INF"]
[docs]
self.ldustm = t["LDUST_MED"]
[docs]
self.ldusts = t["LDUST_SUP"]
[docs]
self.ltirb = t["LUM_TIR_BEST"]
[docs]
self.ltirl = t["LUM_TIR_INF"]
[docs]
self.ltirm = t["LUM_TIR_MED"]
[docs]
self.ltirs = t["LUM_TIR_SUP"]
[docs]
self.massb = t["MASS_BEST"]
[docs]
self.massl = t["MASS_INF"]
[docs]
self.massm = t["MASS_MED"]
[docs]
self.masss = t["MASS_SUP"]
[docs]
self.sfrb = t["SFR_BEST"]
[docs]
self.sfrl = t["SFR_INF"]
[docs]
self.sfrm = t["SFR_MED"]
[docs]
self.sfrs = t["SFR_SUP"]
[docs]
self.ssfrb = t["SSFR_BEST"]
[docs]
self.ssfrl = t["SSFR_INF"]
[docs]
self.ssfrm = t["SSFR_MED"]
[docs]
self.ssfrs = t["SSFR_SUP"]
[docs]
self.lnuv = t["LUM_NUV_BEST"]
[docs]
self.lr = t["LUM_R_BEST"]
[docs]
self.lk = t["LUM_K_BEST"]
[docs]
self.pdfs = np.array(t[pdf_col])
# Define the panels with the binning in redshift an magnitude
if len(range_z) <= 1:
self.range_z = np.quantile(self.zp[(self.zp > -1) & (self.zp < 9)], [0, 0.25, 0.5, 0.75, 1])
else:
self.range_z = range_z
if len(range_mag) <= 1:
self.range_mag = np.quantile(self.mag[(self.mag > 10) & (self.mag < 40)], [0, 0.25, 0.5, 0.75, 1])
else:
self.range_mag = range_mag
# Define the maximums in redshift and magnitude
self.z_min, self.z_max = np.amin(self.range_z), np.amax(self.range_z)
self.mag_min, self.mag_max = np.amin(self.range_mag), np.amax(self.range_mag)
# Define the number panels
# case with several bin of mag
if len(self.range_mag) > 2:
# Number of row with 2 columns
self.nbRowM = int(ceil(float(len(self.range_mag) - 1) / 2.0))
self.nbColM = 2
else:
self.nbColM = 1
self.nbRowM = 1
# case with several bin of z
if len(self.range_z) > 2:
# Number of row with 2 columns
self.nbRowZ = int(ceil(float(len(self.range_z) - 1) / 2.0))
self.nbColZ = 2
else:
self.nbColZ = 1
self.nbRowZ = 1
# Condition for selection
# general condition to select the galaxies in the expected z/mag range
[docs]
self.cond = (
(self.zp > self.z_min)
& (self.zp < self.z_max)
& (self.zml > self.z_min)
& (self.zml < self.z_max)
& (self.mag > self.mag_min)
& (self.mag < self.mag_max)
)
# condition to select stars
[docs]
self.condstar = self.chis < self.chi
# condition to select galaxies
[docs]
self.condgal = ~self.condstar
# condition to select spectroscopic redshifts
[docs]
self.condspec = (self.zs > 0) & (self.zs < 9)
return
[docs]
def save_photoz_plots_pdf(self, filename="all_plots.pdf", **kwargs):
"""
Generate all photoz plotting methods in the class and save them in a single PDF.
Parameters
----------
filename : str, optional
The output PDF filename (default "all_photoz_plots.pdf").
**kwargs : dict, optional
Keyword arguments to pass to the plotting methods. Common keys include:
- bins : int
Number of bins for histograms.
- show_pit : bool
Whether to include PIT histogram in applicable plots.
- show_qq : bool
Whether to include QQ plots in applicable methods.
- savefig : bool
Whether to save individual plots.
Any method that does not accept a particular key will ignore it.
"""
# Get all callable public methods
plot_methods = [
self.title_page,
self.zml_zs,
self.zp_zs,
self.zml_zp,
self.dist_z,
self.dist_chi2,
self.dist_filt,
self.dist_model,
self.dist_ebv,
self.secondpeak,
self.bzk,
self.cumulative68,
self.check_error,
self.errormag,
self.errorz,
self.pit_qq,
]
with PdfPages(filename) as pdf:
for method in plot_methods:
name = method.__name__
print(f"Running plot method: {name}()")
plt.close("all") # clear previous figures
try:
# Try calling with kwargs
method(**kwargs)
except TypeError:
# Fallback if method doesn’t accept kwargs
method()
figs = [plt.figure(i) for i in plt.get_fignums()]
if not figs:
print(f"No figure created by {name}(), skipping.")
continue
for fig in figs:
if not fig.axes:
plt.close(fig)
continue
pdf.savefig(fig, bbox_inches="tight", dpi=300)
plt.close("all")
print(f"All plots saved to {filename}")
[docs]
def save_phys_plots_pdf(self, filename="all_phys_plots.pdf", **kwargs):
"""
Generate all physical parameter plots in the class and save them in a single PDF.
Parameters
----------
filename : str, optional
The output PDF filename (default "all_phys_plots.pdf").
**kwargs : dict, optional
Keyword arguments to pass to the plotting methods. Common keys include:
- bins : int
Number of bins for histograms.
- show_pit : bool
Whether to include PIT histogram in applicable plots.
- show_qq : bool
Whether to include QQ plots in applicable methods.
- savefig : bool
Whether to save individual plots.
Any method that does not accept a particular key will ignore it.
"""
# Get all callable public methods
plot_methods = [
self.title_page,
self.dist_z,
self.dist_chi2,
self.dist_filt,
self.dist_model,
self.dist_ebv,
self.absmag_z,
self.rf_color,
self.william,
self.dist_mass,
self.dist_sfr,
self.dist_ssfr,
self.mass_med_best,
self.sfr_med_best,
self.mass_sfr,
self.mass_z,
self.sfr_z,
self.lnuv_sfr,
self.absmagu_sfr,
self.absmagk_mass,
self.masstolight_z,
]
with PdfPages(filename) as pdf:
for method in plot_methods:
name = method.__name__
print(f"Running plot method: {name}()")
plt.close("all") # clear previous figures
try:
# Try calling with kwargs
method(**kwargs)
except TypeError:
# Fallback if method doesn’t accept kwargs
method()
figs = [plt.figure(i) for i in plt.get_fignums()]
if not figs:
print(f"No figure created by {name}(), skipping.")
continue
for fig in figs:
if not fig.axes:
plt.close(fig)
continue
pdf.savefig(fig, bbox_inches="tight")
plt.close("all")
print(f"All plots saved to {filename}")
[docs]
def title_page(self):
"""Create a title page for diagnostic plots.
This method generates a simple title page summarizing the LePHARE
diagnostic plots. It includes the LePHARE version and the current date.
The information is displayed centered on the page.
The axes are retained (not hidden) so that further annotations or
plot elements can be added later if needed.
"""
today = datetime.now().strftime("%Y-%m-%d")
# Message text
message = f"LePHARE Diagnostics\n" f"LePHARE version: {lp.__version__}\n" f"Date: {today}"
fig, ax = plt.subplots(figsize=(8.5, 11), dpi=300) # standard A4-ish size
# Try to download and display the logo
logo_url = "https://raw.githubusercontent.com/lephare-photoz/lephare-logo/main/Logo/On%20White%20Background/Colour/Digital/LePhareLogo_RGB.png"
try:
with urllib.request.urlopen(logo_url) as url:
image_data = io.BytesIO(url.read())
logo = mpimg.imread(image_data)
# Original dimensions
orig_width, orig_height = 3463, 3000
aspect = orig_height / orig_width # height / width
# Desired width fraction of figure
width_frac = 0.3 # logo spans 50% of figure width
x_center = 0.5
y_top = 0.95 # top margin
# Compute extent while preserving aspect ratio
x_half = width_frac / 2
y_height = width_frac * aspect
ax.imshow(
logo, extent=[x_center - x_half, x_center + x_half, y_top - y_height, y_top], aspect="auto"
)
except Exception as e:
print(f"Could not load LePHARE logo: {e}")
ax.axis("off") # Remove all axes, ticks, and labels
ax.set_frame_on(False) # Remove the border/frame
ax.text(0.5, 0.5, message, ha="center", va="center", fontsize=14, transform=ax.transAxes)
# Optionally, adjust limits for aesthetics
ax.set_xlim(0, 1)
ax.set_ylim(0, 1)
[docs]
def zml_zs(self):
"""
Plot photometric redshifts (median of the PDF) versus spectroscopic redshifts.
This diagnostic compares the median of PDF(z) (`Z_MED`)
derived by LePHARE with the spectroscopic redshifts (`ZSPEC`) for all
objects meeting the selection criteria. It is commonly used to assess
the accuracy, scatter, and presence of catastrophic outliers in
photometric redshift estimates.
The plot is divided into subpanels based on magnitude bins defined in
`range_mag`. Each panel displays `Z_MED` versus `ZSPEC`, with reference
lines indicating the 1:1 relation and typical deviation thresholds.
Notes
-----
- Only objects with valid spectroscopic redshifts (`0 < ZSPEC < 9`) and
within the defined redshift and magnitude ranges are included.
- The number and layout of subpanels depend on the magnitude binning.
- This plot is typically one of the first diagnostics used to check the
performance of photometric redshift estimation.
See Also
--------
zml_zp : Compare median and best-fit photometric redshifts.
zp_zs : Compare best-fit photometric redshifts versus spectroscopic redshifts.
dist_z : Show distribution of redshift differences.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.zml_zs()
# Produces a Z_MED vs ZSPEC scatter plot by magnitude bin.
"""
plt.clf()
f, axarr = plt.subplots(
self.nbRowM, self.nbColM, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
# No space between the figures
plt.subplots_adjust(hspace=0, wspace=0)
# label of the figure
f.text(0.5, 0.04, "$z_{spec}$", ha="center", fontsize=15)
f.text(0.04, 0.5, r"$z_{phot}\; median\; PDF(z)$", va="center", rotation="vertical", fontsize=15)
# Loop over the magnitude bins
for rm in range(len(self.range_mag) - 1):
# Define the subplots and pass the panel
ax = axarr[int(rm // 2.0), int(rm % 2)]
# mag limits
magmin = self.range_mag[rm]
magmax = self.range_mag[rm + 1]
ax.axis([self.z_min, self.z_max, self.z_min, self.z_max])
conda = self.cond & self.condspec & (self.mag > magmin) & (self.mag < magmax)
ax.scatter(
self.zs[conda], self.zml[conda], s=1, color="b", alpha=0.5, marker="s", rasterized=True
)
# Check that we have some sources before performing statistics
ngal = len(self.zml[conda])
if ngal > 0:
# statistics
arg_bias = self.zml[conda] - self.zs[conda]
arg_std = arg_bias / (1.0 + self.zs[conda])
nmad = 1.4821 * np.median(abs(arg_std))
cond_outl = abs(arg_std) > 0.15
outl_rate = len(arg_std[cond_outl]) / float(ngal)
ax.annotate(
r"$"
+ str(round(magmin, 2))
+ " < mag < "
+ str(round(magmax, 2))
+ "$ \n"
+ "$N_{gal} = "
+ str(ngal)
+ "$ \n"
+ r"$\eta ="
+ str(round(100 * outl_rate, 2))
+ " \\%$\n"
+ r"$ \sigma_{\Delta z /(1+z)} = "
+ str(round(nmad, 5))
+ "$",
xy=(0.05 * self.z_max, 0.65 * self.z_max),
color="black",
fontsize=12,
)
# Trace the limits 0.15(1+z)
x_zs = np.array([0, 6])
ax.plot(x_zs, x_zs * 1.15 + 0.15, "c--")
ax.plot(x_zs, x_zs, "r-")
ax.plot(x_zs, x_zs * 0.85 - 0.15, "c--")
return
[docs]
def zp_zs(self):
"""
Plot photometric redshifts (minimum chi-squared value) versus spectroscopic redshifts.
This diagnostic compares the best-fit photometric redshifts (`Z_BEST`)
obtained from the minimum χ² solution in LePHARE with the spectroscopic
redshifts (`ZSPEC`). It is useful for assessing the quality and bias of
the primary redshift solution and identifying catastrophic outliers.
The figure is divided into subpanels according to the magnitude bins
defined in `range_mag`. Each panel shows `Z_BEST` versus `ZSPEC` along
with reference lines for the one-to-one relation and typical redshift
deviation thresholds (e.g., Δz / (1 + z) = ±0.05).
Notes
-----
- Only sources with valid spectroscopic redshifts (`0 < ZSPEC < 9`) and
within the defined redshift and magnitude ranges are included.
- The panel layout is determined by the number of magnitude bins.
- This plot complements :func:`zml_zs` by using the minimum χ² redshift
rather than the PDF median.
See Also
--------
zml_zs : Compare median photometric redshifts versus spectroscopic redshifts.
zml_zp : Compare median and best-fit photometric redshifts.
dist_z : Show distribution of redshift differences.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.zp_zs()
# Produces a Z_BEST vs ZSPEC scatter plot by magnitude bin.
"""
plt.clf()
f, axarr = plt.subplots(
self.nbRowM, self.nbColM, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
# No space between the figures
plt.subplots_adjust(hspace=0, wspace=0)
# label of the figure
f.text(0.5, 0.04, "$z_{spec}$", ha="center", fontsize=15)
f.text(0.04, 0.5, r"$z_{phot}\; median\; PDF(z)$", va="center", rotation="vertical", fontsize=15)
# Loop over the magnitude bins
for rm in range(len(self.range_mag) - 1):
# Define the subplots and pass the panel
ax = axarr[int(rm // 2.0), int(rm % 2)]
# mag limits
magmin = self.range_mag[rm]
magmax = self.range_mag[rm + 1]
ax.axis([self.z_min, self.z_max, self.z_min, self.z_max])
conda = self.cond & self.condspec & (self.mag > magmin) & (self.mag < magmax)
ax.scatter(self.zs[conda], self.zp[conda], s=1, color="b", alpha=0.5, marker="s", rasterized=True)
# Check that we have some sources before performing statistics
ngal = len(self.zp[conda])
if ngal > 0:
# statistics
arg_bias = self.zp[conda] - self.zs[conda]
arg_std = arg_bias / (1.0 + self.zs[conda])
nmad = 1.4821 * np.median(abs(arg_std))
cond_outl = abs(arg_std) > 0.15
outl_rate = len(arg_std[cond_outl]) / float(ngal)
ax.annotate(
r"$"
+ str(round(magmin, 2))
+ " < mag < "
+ str(round(magmax, 2))
+ "$ \n"
+ "$N_{gal} = "
+ str(ngal)
+ "$ \n"
+ r"$\eta ="
+ str(round(100 * outl_rate, 2))
+ " \\%$\n"
+ r"$ \sigma_{\Delta z /(1+z)} = "
+ str(round(nmad, 5))
+ "$",
xy=(0.05 * self.z_max, 0.65 * self.z_max),
color="black",
fontsize=12,
)
# Trace the limits 0.15(1+z)
x_zs = np.array([0, 6])
ax.plot(x_zs, x_zs * 1.15 + 0.15, "c--")
ax.plot(x_zs, x_zs, "r-")
ax.plot(x_zs, x_zs * 0.85 - 0.15, "c--")
return
[docs]
def zml_zp(self):
"""
Plot photometric redshifts (median of the PDF) versus photometric redshifts (minimum χ²).
This diagnostic compares the two photometric redshift estimates produced by
LePHARE: the median of the posterior probability distribution (`Z_MED`)
and the best-fit value from the minimum χ² solution (`Z_BEST`).
The comparison helps identify systematic biases between the two estimators
and assess the stability of redshift determinations.
The figure is divided into subpanels based on magnitude bins defined in
`range_mag`. Each subpanel shows `Z_MED` versus `Z_BEST` with a 1:1
reference line and scatter indicating consistency between the two measures.
Notes
-----
- Only sources within the defined redshift and magnitude ranges are included.
- The plot can reveal populations where the PDF and best-fit solutions diverge,
such as multi-peaked PDFs or poor fits.
- Strong deviations from the 1:1 line indicate inconsistency between the
χ²-based and PDF-based estimates.
See Also
--------
zml_zs : Compare median photometric redshifts versus spectroscopic redshifts.
zp_zs : Compare best-fit photometric redshifts versus spectroscopic redshifts.
dist_z : Show distribution of redshift distribution.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.zml_zp()
# Produces a Z_MED vs Z_BEST comparison plot by magnitude bin.
"""
plt.clf()
f, axarr = plt.subplots(
self.nbRowM, self.nbColM, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
# No space between the figures
plt.subplots_adjust(hspace=0, wspace=0)
# label of the figure
f.text(0.5, 0.04, r"$z_{phot}\; minimum\; \chi^2$", ha="center", fontsize=15)
f.text(0.04, 0.5, r"$z_{phot}\; median\; PDF(z)$", va="center", rotation="vertical", fontsize=15)
# Loop over the magnitude bins
for rm in range(len(self.range_mag) - 1):
# Define the subplots and pass the panel
ax = axarr[int(rm // 2.0), int(rm % 2)]
# mag limits
magmin = self.range_mag[rm]
magmax = self.range_mag[rm + 1]
# set the axis
ax.axis([self.z_min, self.z_max, self.z_min, self.z_max])
# new condition with the magnitude range
conda = self.cond & self.condgal & (self.mag > magmin) & (self.mag < magmax)
# Plot photo-z versus spec-z
ax.scatter(self.zp[conda], self.zml[conda], s=1, color="b", alpha=0.5, marker="s")
# Trace the limits 0.15(1+z)
x_zs = np.array([0, 6])
ax.plot(x_zs, x_zs * 1.15 + 0.15, "c--")
ax.plot(x_zs, x_zs, "r-")
ax.plot(x_zs, x_zs * 0.85 - 0.15, "c--")
# labels
ax.annotate(
"$" + str(round(magmin, 2)) + " < mag < " + str(round(magmax, 2)) + "$",
xy=(0.1 * self.z_max, 0.8 * self.z_max),
color="black",
fontsize=15,
)
return
[docs]
def dist_z(self, nstep=20):
"""
Plot the redshift distribution of the spectroscopic and photometric samples.
This diagnostic shows histograms of the spectroscopic redshifts (`ZSPEC`),
the photometric redshifts from the minimum χ² solution (`Z_BEST`), and the
median of the PDF (`Z_MED`). It provides an overview of the redshift coverage
and potential biases between the photometric and spectroscopic distributions.
Parameters
----------
nstep : int, optional
Number of bins to use in the redshift histograms. Default is 20.
Notes
-----
- Only objects satisfying the redshift and magnitude selection criteria are used.
- Overlaid histograms allow direct comparison of the redshift distributions.
- Discrepancies between spectroscopic and photometric distributions may indicate
systematic shifts or incomplete coverage in certain redshift ranges.
See Also
--------
zml_zs : Compare median photometric redshifts versus spectroscopic redshifts.
zp_zs : Compare best-fit photometric redshifts versus spectroscopic redshifts.
zml_zp : Compare median and best-fit photometric redshifts.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.dist_z(nstep=30)
# Produces histograms of ZSPEC, Z_BEST, and Z_MED distributions.
"""
## Create a figure with an array with nbRowM*nbColM subpanels
plt.clf()
f, axarr = plt.subplots(
self.nbRowM, self.nbColM, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
# No space between the figures
plt.subplots_adjust(hspace=0, wspace=0)
# label of the figure
f.text(0.5, 0.04, "$Redshift$", ha="center", fontsize=15)
f.text(0.04, 0.5, "$N_{normalized}$", va="center", rotation="vertical", fontsize=15)
# Loop over the magnitude bins
for rm in range(len(self.range_mag) - 1):
# No legend
leg = 0
if rm == 1:
leg = 1 # legend in the first panel
# Define the subplots and pass the panel
ax = axarr[int(rm // 2.0), int(rm % 2)]
# mag limits
magmin = self.range_mag[rm]
magmax = self.range_mag[rm + 1]
# set the axis
ax.axis([0, self.z_max, 0, 2.9])
# new condition with the magnitude range
conda = self.cond & self.condgal & (self.mag > magmin) & (self.mag < magmax) & (self.zml > 0)
# Check that some object are present
if len(self.zp[conda]) > 0:
# Histogram with the photometric redshifts median PDF
ax.hist(
self.zml[conda],
bins=nstep,
histtype="stepfilled",
density=1,
color="b",
label=r"$z_{phot}\; median\; PDF(z)$",
)
# Histogram with the photometric redshifts minimum chi2
ax.hist(
self.zp[conda],
bins=nstep,
histtype="stepfilled",
density=1,
color="r",
alpha=0.5,
label=r"$z_{phot}\; minimum\; \chi^2$",
)
# new condition with the magnitude range and spectro-z
conda = self.cond & self.condspec & (self.mag > magmin) & (self.mag < magmax)
# Check that some object are present
if len(self.zp[conda]) > 0:
# Histogram with the photometric redshifts zbest
ax.hist(
self.zs[conda],
bins=nstep,
histtype="stepfilled",
density=1,
color="g",
alpha=0.2,
label="$z_{spectro}$",
)
# labels
ax.annotate(
"$" + str(round(magmin, 2)) + " < mag < " + str(round(magmax, 2)) + "$",
xy=(0.5 * self.z_max, 1.5),
color="black",
fontsize=14,
)
# print the legend
if leg == 1:
ax.legend(prop={"size": 8})
return
[docs]
def dist_chi2(self, nstep=20):
"""
Plot the reduced χ² distribution as a function of observed magnitude.
This diagnostic shows the distribution of the minimum χ² values
(`CHI_BEST`) for different magnitude bins, helping to assess the quality
of SED fits across brightness ranges. Poorer fits at faint magnitudes
may indicate increased photometric noise or template limitations.
We estimate the reduced χ² by dividing the minimum χ² by the number of
filters used minus 1
Parameters
----------
nstep : int, optional
Number of bins to use in the χ² histograms. Default is 20.
Notes
-----
- The figure is typically divided into subpanels according to `range_mag`.
- Only objects meeting the redshift and magnitude selection criteria are included.
- A broad or asymmetric χ² distribution at a given magnitude may suggest
underestimated photometric errors or mismatched templates.
See Also
--------
dist_z : Show redshift distributions.
dist_model : Show distribution of best-fit SED model identifiers.
dist_ebv : Show distribution of dust extinction values.
check_error : Verify internal photometric error estimates.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.dist_chi2(nstep=30)
# Produces χ² histograms by magnitude bin.
"""
## Create a figure with an array with nbRowM*nbColM subpanels
plt.clf()
f, axarr = plt.subplots(
self.nbRowM, self.nbColM, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
# No space between the figures
plt.subplots_adjust(hspace=0, wspace=0)
# label of the figure
f.text(0.5, 0.04, r"log($\chi^2$)", ha="center", fontsize=15)
f.text(0.04, 0.5, "$N_{normalized}$", va="center", rotation="vertical", fontsize=15)
# Loop over the magnitude bins
for rm in range(len(self.range_mag) - 1):
# No legend
leg = 0
if rm == 2:
leg = 1 # legend in the first panel
# Define the subplots and pass the panel
ax = axarr[int(rm // 2.0), int(rm % 2)]
# mag limits
magmin = self.range_mag[rm]
magmax = self.range_mag[rm + 1]
# set the axis
ax.axis([-2, 2.9, 0, 1.49])
# new condition with the magnitude range for galaxies, chi gal< chi star
conda = (
self.cond
& self.condgal
& (self.mag > magmin)
& (self.mag < magmax)
& (self.nbFilt > 2)
& (self.chi > 0)
)
if len(self.zp[conda]) > 1:
chireduit = self.chi[conda] / (self.nbFilt[conda] - 1)
# Histogram of log (reduced chi2 gal)
ax.hist(
np.log10(chireduit),
bins=nstep,
histtype="stepfilled",
density=1,
color="b",
label=r"$\chi^2_{gal} \; if \; \chi^2_{gal}<\chi^2_{stars}$",
)
# new condition with the magnitude range for stars, chi star < chi gal
conda = (
self.cond
& self.condstar
& (self.mag > magmin)
& (self.mag < magmax)
& (self.nbFilt > 2)
& (self.chi > 0)
)
if len(self.zp[conda]) > 1:
chireduit = self.chis[conda] / (self.nbFilt[conda] - 1)
# Histogram of log (reduced chi2 gal)
ax.hist(
np.log10(chireduit),
bins=nstep,
histtype="stepfilled",
density=1,
color="r",
alpha=0.5,
label=r"$\chi^2_{star} \; if \; \chi^2_{gal}>\chi^2_{stars}$",
)
# print the legend
if leg == 1:
ax.legend(prop={"size": 8})
# labels
ax.annotate(
"$" + str(round(magmin, 2)) + " < mag < " + str(round(magmax, 2)) + "$",
xy=(-1.8, 1.2),
color="black",
fontsize=15,
)
return
[docs]
def dist_filt(self, nstep=10):
"""
Plot the distribution of the number of filters used in each fit.
This diagnostic shows how many photometric bands were included by LePHARE
when fitting each object. It helps assess data completeness and detect
cases where limited filter coverage may affect the reliability of photometric
redshifts or SED parameters.
Parameters
----------
nstep : int, optional
Number of bins to use in the histogram of `NBAND_USED`. Default is 10.
Notes
-----
- Uses the `NBAND_USED` column from the LePHARE output table.
- Only sources satisfying the redshift and magnitude selection criteria are included.
- A concentration of objects with a low number of filters may indicate missing
photometry or incomplete coverage in certain bands.
See Also
--------
dist_chi2 : Show χ² distribution by magnitude bin.
dist_model : Show distribution of best-fit model identifiers.
dist_ebv : Show distribution of dust extinction values.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.dist_filt(nstep=15)
# Produces a histogram of the NBAND_USED distribution.
"""
## Create a figure with an array with nbRowM*nbColM subpanels
plt.clf()
f, axarr = plt.subplots(
self.nbRowZ, self.nbColZ, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
# No space between the figures
plt.subplots_adjust(hspace=0, wspace=0)
# label of the figure
f.text(0.5, 0.04, "Number of filters", ha="center", fontsize=15)
f.text(0.04, 0.5, "N", va="center", rotation="vertical", fontsize=15)
# Loop over the redshift bins
for rm in range(len(self.range_z) - 1):
# Define the subplots and pass the panel
ax = axarr[int(rm // 2.0), int(rm % 2)]
# mag limits
zmin = self.range_z[rm]
zmax = self.range_z[rm + 1]
# set the axis
ax.axis([0, 50, 0, 3])
# new condition with the magnitude range
conda = self.cond & (self.zp > zmin) & (self.zp < zmax)
# Histogram with the number of filters
ax.hist(self.nbFilt[conda], bins=nstep, histtype="stepfilled", density=1, color="r")
# labels
ax.annotate(
"$" + str(round(zmin, 2)) + " < z < " + str(round(zmax, 2)) + "$",
xy=(0.1, 1),
color="black",
fontsize=15,
)
return
[docs]
def dist_model(self, nstep=20):
"""
Plot the distribution of best-fit SED model identifiers.
This diagnostic displays the frequency of SED templates (`MOD_BEST`)
selected by LePHARE as the best-fitting models across the sample. It helps
visualize which galaxy, AGN, or stellar templates dominate the fits and
can reveal biases toward certain template types.
Parameters
----------
nstep : int, optional
Number of bins to use in the histogram of model identifiers. Default is 20.
Notes
-----
- Uses the `MOD_BEST` column from the LePHARE output table.
- Only sources satisfying the redshift and magnitude selection criteria are included.
- The model identifiers correspond to entries in the LePHARE template library.
- A dominance of certain templates may indicate population trends or a lack of
diversity in the input SED library.
See Also
--------
dist_ebv : Show distribution of best-fit dust extinction values.
dist_chi2 : Show χ² distribution by magnitude bin.
dist_filt : Show number of filters used in each fit.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.dist_model(nstep=25)
# Produces a histogram of MOD_BEST template identifiers.
"""
## Create a figure with an array with nbRowM*nbColM subpanels
plt.clf()
f, axarr = plt.subplots(
self.nbRowZ, self.nbColZ, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
# No space between the figures
plt.subplots_adjust(hspace=0, wspace=0)
# label of the figure
f.text(0.5, 0.04, "best fit template", ha="center", fontsize=15)
f.text(0.04, 0.5, "N", va="center", rotation="vertical", fontsize=15)
# Loop over the redshift bins
for rm in range(len(self.range_z) - 1):
# Define the subplots and pass the panel
zmin = self.range_z[rm]
zmax = self.range_z[rm + 1]
ax = axarr[int(rm // 2.0), int(rm % 2)]
# set the axis
ax.axis([0, 70, 0, 0.3])
# new condition with the magnitude range
conda = self.cond & self.condgal & (self.zp > zmin) & (self.zp < zmax)
# check if some objects exist
if len(self.mod[conda]) > 1:
# Histogram of the models
ax.hist(self.mod[conda], bins=nstep, histtype="stepfilled", density=1, color="b")
# labels
ax.annotate(
"$" + str(round(zmin, 2)) + " < z < " + str(round(zmax, 2)) + "$",
xy=(10, 0.25),
color="black",
fontsize=15,
)
return
[docs]
def dist_ebv(self, nstep=10):
"""
Plot the distribution of best-fit dust extinction values (E(B-V)).
This diagnostic shows the range and frequency of dust extinction
(`EBV_BEST`) values derived by LePHARE for the sample. It helps assess
typical dust content and the effects of extinction on photometric redshift
or SED fitting results.
Parameters
----------
nstep : int, optional
Number of bins to use in the histogram of E(B-V) values. Default is 10.
Notes
-----
- Uses the `EBV_BEST` column from the LePHARE output table.
- Only sources meeting the redshift and magnitude selection criteria are included.
- Peaks or concentrations in E(B-V) may indicate preferred extinction values
in the template set or the typical dust content of the population.
See Also
--------
dist_model : Show distribution of best-fit SED model identifiers.
dist_chi2 : Show χ² distribution by magnitude bin.
dist_filt : Show number of filters used in each fit.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.dist_ebv(nstep=15)
# Produces a histogram of EBV_BEST values.
"""
## Create a figure with an array with nbRowZ*nbColZ subpanels
plt.clf()
f, axarr = plt.subplots(
self.nbRowZ, self.nbColZ, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
# No space between the figures
plt.subplots_adjust(hspace=0, wspace=0)
# label of the figure
f.text(0.5, 0.04, "E(B-V)", ha="center", fontsize=15)
f.text(0.04, 0.5, "N", va="center", rotation="vertical", fontsize=15)
# Loop over the redshift bins
for rm in range(len(self.range_z) - 1):
# Define the subplots and pass the panel
zmin = self.range_z[rm]
zmax = self.range_z[rm + 1]
ax = axarr[int(rm // 2.0), int(rm % 2)]
# set the axis
ax.axis([0, 0.6, 0, 20])
# new condition with the magnitude range
conda = self.cond & self.condgal & (self.zp > zmin) & (self.zp < zmax)
# check if some objects exist
if len(self.ebv[conda]) > 1:
# Histogram with the E(B-V)
ax.hist(self.ebv[conda], bins=nstep, histtype="stepfilled", density=1, color="b")
# labels
ax.annotate(
"$" + str(round(zmin, 2)) + " < z < " + str(round(zmax, 2)) + "$",
xy=(0.02, 15),
color="black",
fontsize=15,
)
return
[docs]
def secondpeak(self, nstep=20):
"""
Plot the redshift distribution of the secondary photometric redshift peak.
This diagnostic examines objects where LePHARE identifies a secondary
photometric redshift solution (`Z_SEC`) and compares it with the
primary solutions (`Z_BEST` and `Z_MED`). It helps identify
multi-modal PDFs and potential catastrophic redshift outliers.
Parameters
----------
nstep : int, optional
Number of bins to use in the histogram of secondary peak redshifts.
Default is 20.
Notes
-----
- Only sources with a valid secondary peak (`Z_SEC`) are included.
- The figure may be divided into subpanels based on magnitude bins.
- Comparison of primary and secondary peaks can reveal degeneracies in SED fitting
and ambiguous photometric redshift solutions.
- Useful for assessing the reliability of photometric redshift catalogs
and flagging objects with uncertain redshifts.
See Also
--------
zml_zs : Compare median photometric redshift versus spectroscopic redshift.
zp_zs : Compare best-fit photometric redshift versus spectroscopic redshift.
zml_zp : Compare median and best-fit photometric redshifts.
dist_z : Show overall redshift distributions.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.secondpeak(nstep=25)
# Produces a histogram of secondary peak redshifts compared with primary solutions.
"""
## Create a figure with an array with nbRowM*nbColM subpanels
plt.clf()
f, axarr = plt.subplots(
self.nbRowM, self.nbColM, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
# No space between the figures
plt.subplots_adjust(hspace=0, wspace=0)
# label of the figure
f.text(0.5, 0.04, "$Redshift$", ha="center", fontsize=15)
f.text(0.04, 0.5, "$N_{normalized}$", va="center", rotation="vertical", fontsize=15)
# Loop over the magnitude bins
for rm in range(len(self.range_mag) - 1):
# No legend
leg = 0
if rm == 4:
leg = 1 # legend in the first panel
# Define the subplots and pass the panel
ax = axarr[int(rm // 2.0), int(rm % 2)]
# mag limits
magmin = self.range_mag[rm]
magmax = self.range_mag[rm + 1]
# set the axis
ax.axis([0, self.z_max, 0, 5])
# new condition with the magnitude range
conda = self.cond & self.condgal & (self.mag > magmin) & (self.mag < magmax) & (self.zp2 > 0)
# If some sources
if len(self.zp2[conda]) > 1:
# Histogram with the photometric redshifts median PDF
ax.hist(
self.zp2[conda],
bins=nstep,
histtype="stepfilled",
density=1,
color="b",
label=r"$z_{phot}\; second\; peak$",
)
# Histogram with the photometric redshifts minimum chi2
ax.hist(
self.zp[conda],
bins=nstep,
histtype="stepfilled",
density=1,
color="r",
alpha=0.5,
label=r"$z_{phot}\; minimum \; \chi^2$",
)
# labels
ax.annotate(
"$" + str(round(magmin, 2)) + " < mag < " + str(round(magmax, 2)) + "$",
xy=(0.2, 3.5),
color="black",
fontsize=13,
)
# print the legend
if leg == 1:
ax.legend(prop={"size": 8})
return
[docs]
def bzk(self):
"""
Plot BzK diagram for galaxy/star classification per redshift bin.
This diagnostic produces BzK color-color diagrams to separate star-forming
galaxies, passive galaxies, and stars based on the observed B, z, and K
magnitudes. The plots are generated for different redshift bins defined
in `range_z`, allowing assessment of population segregation at different epochs.
Notes
-----
- Uses the observed magnitudes in B (`magb`), z (`magz`), and K (`magk`) filters.
- Only sources meeting the redshift and magnitude selection criteria are included.
- Overlaid regions indicate typical BzK selection areas for star-forming galaxies,
passive galaxies, and stars.
- Helps identify misclassified objects or potential contamination by stars
in galaxy samples.
See Also
--------
dist_filt : Show number of filters used in each fit.
dist_model : Show distribution of best-fit SED model identifiers.
secondpeak : Examine secondary photometric redshift peaks.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.bzk()
# Produces BzK color-color diagrams per redshift bin.
"""
plt.clf()
f, axarr = plt.subplots(2, 2, sharex=True, sharey=True, figsize=(12, 8), squeeze=False)
plt.subplots_adjust(hspace=0, wspace=0)
f.text(0.5, 0.04, "$B-z$", ha="center", fontsize=15)
f.text(0.04, 0.5, "$z-K$", va="center", rotation="vertical", fontsize=15)
bzk_conf = [[0.001, 1.4, 0, 0], [1.4, 2.7, 0, 1], [2.7, 6, 1, 0], [-0.001, 6, 1, 1]]
for rm in range(4):
zmin = bzk_conf[rm][0]
zmax = bzk_conf[rm][1]
ax = axarr[bzk_conf[rm][2], bzk_conf[rm][3]]
ax.axis([-0.5, 5.5, -1.0, 4.5])
if rm < 3:
conda = self.cond & self.condgal & (self.zp > zmin) & (self.zp < zmax)
else:
conda = self.cond & self.condstar
if (self.bFilt >= 0) and (self.zFilt >= 0) and (self.kFilt >= 0):
colx = self.magb - self.magz
coly = self.magz - self.magk
ax.scatter(colx[conda], coly[conda], s=1, color="b", alpha=0.2, marker="s", rasterized=True)
ax.plot([-1, 6], [-0.5, 6])
if rm < 3:
ax.annotate(f"${zmin:.2f} < z < {zmax:.2f}$", xy=(0.1, 3), color="black", fontsize=15)
else:
ax.annotate(r"$\chi_s < \chi_g$", xy=(0.1, 3), color="black", fontsize=15)
return
[docs]
def absmag_z(self):
"""
Plot B-band absolute magnitude versus redshift for the sample.
This diagnostic shows the relation between absolute magnitudes
and redshift, providing insight into luminosity evolution, selection effects, and completeness
limits of the survey.
Notes
-----
- Uses the B-band absolute magnitudes
- Only objects within the redshift and magnitude selection ranges are included.
- The figure can be divided into subpanels based on magnitude or redshift bins.
- Helps identify trends such as brightening or dimming with redshift
and the distribution of galaxies in the luminosity-redshift plane.
See Also
--------
rf_color : Rest-frame color versus redshift.
bzk : BzK color-color diagram for star/galaxy classification.
dist_z : Redshift distributions.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.absmag_z()
# Produces plots of absolute magnitude versus redshift.
"""
plt.clf()
f, axarr = plt.subplots(
self.nbRowM, self.nbColM, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
plt.subplots_adjust(hspace=0, wspace=0)
f.text(0.5, 0.04, "$Redshift$", ha="center", fontsize=15)
f.text(0.04, 0.5, "$M_B$", va="center", rotation="vertical", fontsize=15)
for rm in range(len(self.range_mag) - 1):
ax = axarr[int(rm // 2.0), int(rm % 2)]
magmin = self.range_mag[rm]
magmax = self.range_mag[rm + 1]
ax.axis([self.z_min, self.z_max, -16, -25.9])
conda = self.cond & self.condgal & (self.mag > magmin) & (self.mag < magmax)
ax.scatter(self.zp[conda], self.mabsb[conda], s=1, color="b", alpha=0.2, marker="s")
ax.annotate(f"${magmin:.2f} < mag < {magmax:.2f}$", xy=(0.1, -25), color="black", fontsize=15)
return
[docs]
def rf_color(self):
"""
Plot rest-frame colors U-R versus absolute magnitude R.
This diagnostic shows the relation between rest-frame colors
and absolute magnitudes in different bands, allowing analysis of galaxy
populations, color bimodality, and stellar population trends.
Notes
-----
- Uses absolute magnitudes (`MAG_ABS()`) and rest-frame colors derived from
the best-fit templates.
- Only objects within the redshift and magnitude selection ranges are included.
- Can reveal red sequence and blue cloud populations in the color–magnitude diagram.
- Subpanels can be used to separate objects by magnitude or redshift bins.
See Also
--------
absmag_z : Absolute magnitude versus redshift.
bzk : BzK color-color diagram for star/galaxy classification.
dist_model : Distribution of best-fit SED models.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.rf_color()
# Produces rest-frame color versus absolute magnitude plots.
"""
plt.clf()
f, axarr = plt.subplots(
self.nbRowM, self.nbColM, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
plt.subplots_adjust(hspace=0, wspace=0)
f.text(0.5, 0.04, "$M_R$", ha="center", fontsize=15)
f.text(0.04, 0.5, "$M_U-M_R$", va="center", rotation="vertical", fontsize=15)
for rm in range(len(self.range_z) - 1):
zmin = self.range_z[rm]
zmax = self.range_z[rm + 1]
ax = axarr[int(rm // 2.0), int(rm % 2)]
ax.axis([-24.9, -16, -0.7, 2.4])
conda = self.cond & self.condgal & (self.zp > zmin) & (self.zp < zmax)
if self.uFilt >= 0 and self.rFilt >= 0:
magx = self.mabsr
coly = self.mabsu - self.mabsr
ax.scatter(magx[conda], coly[conda], s=1, color="b", alpha=0.2, marker="s", rasterized=True)
ax.annotate(f"${zmin:.2f} < z < {zmax:.2f}$", xy=(-20, 1), color="black", fontsize=15)
return
[docs]
def william(self):
"""
Create a rest-frame U-R, R-J color-color plot.
This diagnostic reproduces the rest-frame color-color diagrams of
''William'', which are commonly used to separate galaxy populations
(e.g., quiescent vs star-forming) and study stellar population properties.
Notes
-----
- Uses rest-frame magnitudes or colors derived from the best-fit templates.
- Only objects within the redshift and magnitude selection ranges are included.
- Subpanels can be applied based on redshift or magnitude bins to study
evolution with cosmic time.
- Helps visualize the bimodality of galaxy populations in rest-frame color space.
See Also
--------
rf_color : Rest-frame color versus absolute magnitude.
absmag_z : Absolute magnitude versus redshift.
bzk : BzK color-color diagram for star/galaxy classification.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.william()
# Produces a William style rest-frame color-color plot.
"""
plt.clf()
f, axarr = plt.subplots(
self.nbRowM, self.nbColM, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
plt.subplots_adjust(hspace=0, wspace=0)
f.text(0.5, 0.04, "$M_R-M_J$", ha="center", fontsize=15)
f.text(0.04, 0.5, "$M_U-M_R$", va="center", rotation="vertical", fontsize=15)
for rm in range(len(self.range_z) - 1):
zmin = self.range_z[rm]
zmax = self.range_z[rm + 1]
ax = axarr[int(rm // 2.0), int(rm % 2)]
ax.axis([-2.1, 1.9, 0, 2.5])
conda = self.cond & self.condgal & (self.zp > zmin) & (self.zp < zmax)
if self.uFilt >= 0 and self.rFilt >= 0 and self.jFilt >= 0:
colx = self.mabsr - self.mabsj
coly = self.mabsu - self.mabsr
ax.scatter(colx[conda], coly[conda], s=1, color="b", alpha=0.2, marker="s", rasterized=True)
ax.annotate(f"${zmin:.2f} < z < {zmax:.2f}$", xy=(0.1, 2), color="black", fontsize=15)
return
[docs]
def cumulative68(self, nstep=20):
"""
Check that 68% of spectroscopic redshifts fall within the 68% photometric error interval.
This diagnostic evaluates the accuracy of photometric redshift uncertainties
by computing the fraction of objects for which the spectroscopic redshift
(`ZSPEC`) lies within the 68% confidence interval of the photometric redshift
PDF (`Z_MED68_LOW` to `Z_MED68_HIGH` or `Z_BEST68_LOW` to `Z_BEST68_HIGH`).
It helps validate the reliability of the estimated errors.
Parameters
----------
nstep : int, optional
Number of bins to use when plotting cumulative distributions. Default is 20.
Notes
-----
- Only objects with valid spectroscopic redshifts are considered.
- A well-calibrated photometric redshift error distribution should contain
~68% of spectroscopic redshifts within the 68% interval.
- Deviations from 68% indicate under- or overestimated photometric uncertainties.
See Also
--------
errormag : Compare magnitude errors versus χ².
errorz : Compare photometric redshift errors with scatter.
zml_zs : Compare median photometric redshift versus spectroscopic redshift.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.cumulative68(nstep=30)
# Evaluates the fraction of spectroscopic redshifts within the 68% photometric error interval.
"""
plt.clf()
f, axarr = plt.subplots(
self.nbRowM, self.nbColM, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
plt.subplots_adjust(hspace=0, wspace=0)
f.text(0.5, 0.04, "$abs((z_{spec}-z_{phot})/z_{phot} uncertainties)$", ha="center", fontsize=15)
f.text(0.04, 0.5, "$cumulative(N)$", va="center", rotation="vertical", fontsize=15)
for rm in range(len(self.range_mag) - 1):
ax = axarr[int(rm // 2.0), int(rm % 2)]
magmin = self.range_mag[rm]
magmax = self.range_mag[rm + 1]
ax.axis([0, 9.9, 0, 0.99])
diffzml = abs(self.zml - self.zs) / np.maximum(
abs(self.zmlu68 - self.zml), abs(self.zmll68 - self.zml)
)
diffzp = abs(self.zp - self.zs) / np.maximum(abs(self.zu68 - self.zp), abs(self.zl68 - self.zp))
condazml = (
self.cond
& self.condspec
& (self.mag > magmin)
& (self.mag < magmax)
& (diffzml < 10)
& (self.zml > 0)
)
condazp = (
self.cond
& self.condspec
& (self.mag > magmin)
& (self.mag < magmax)
& (diffzp < 10)
& (self.zp > 0)
)
if len(self.zml[condazml]) > 1 and len(self.zp[condazp]) > 1:
ax.hist(
diffzp[condazp],
bins=nstep,
histtype="step",
density=1,
color="b",
cumulative=True,
label=r"$z_{phot}\; minimum \; \chi^2$",
)
ax.hist(
diffzml[condazml],
bins=nstep,
histtype="step",
density=1,
color="r",
cumulative=True,
label=r"$z_{phot}\; median \; PDF$",
)
ax.axvline(x=1, color="r", linestyle="--")
ax.axhline(y=0.68, color="r", linestyle="--")
ax.legend(prop={"size": 8})
ax.annotate(f"${magmin:.2f} < mag < {magmax:.2f}$", xy=(0.1, 0.80), color="black", fontsize=15)
return
[docs]
def check_error(self):
"""
Check photometric redshift errors as a function of magnitude.
This diagnostic evaluates how the estimated photometric redshift
uncertainties change as a function of observed magnitude.
Notes
-----
- Uses the columns `Z_BEST68_LOW`, `Z_BEST68_HIGH`, `Z_MED68_LOW`, `Z_MED68_HIGH`,
and `ZSPEC`.
- The plot shows error distributions as a function of the selected magnitude band.
- Can reveal trends where photometric errors are underestimated at faint magnitudes
or overestimated at bright magnitudes.
See Also
--------
cumulative68 : Check fraction of spectroscopic redshifts within the 68% photometric error.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.check_error()
# Produces a plot of photometric redshift errors versus observed magnitude.
"""
# Clear the figure
plt.clf()
# Set the axis
plt.axis([15, 27, -1, 1])
# 1 sigma error on the photo-z
diffu = self.zu68 - self.zp
diffl = self.zl68 - self.zp
conda = self.cond & self.condgal
# Error versus mag
plt.scatter(self.mag[conda], diffu[conda], s=1, color="b", alpha=0.5, marker="s", rasterized=True)
plt.scatter(self.mag[conda], diffl[conda], s=1, color="r", alpha=0.5, marker="s", rasterized=True)
# 68%
plt.plot([15, 27], [0, 0], color="r") # Adapte les limites x à ton graphique
plt.xlabel("magnitude", fontsize=18, labelpad=13)
plt.ylabel(r"$68\% \; photo-z \; error$", fontsize=18, labelpad=13)
return
[docs]
def errormag(self):
"""
Plot photometric redshift errors as a function of observed magnitude.
This diagnostic shows how the estimated photometric redshift errors
(e.g., derived from the PDF width or χ² fitting) vary with object brightness.
It is useful for assessing the reliability of redshift estimates across
the magnitude range.
Notes
-----
- Uses the selected magnitude (`mag`) and corresponding photometric redshift
errors from the LePHARE output.
- Only objects within the redshift and magnitude selection criteria are included.
- Typically, errors increase for fainter objects due to lower signal-to-noise.
See Also
--------
check_error : Evaluate 68% error coverage as a function of magnitude.
errorz : Compare photometric redshift errors as a function of redshift.
cumulative68 : Check fraction of spectroscopic redshifts within the 68% photometric error.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.errormag()
# Produces a plot of photometric redshift errors versus magnitude.
"""
plt.clf()
f, axarr = plt.subplots(
self.nbRowM, self.nbColM, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
plt.subplots_adjust(hspace=0, wspace=0)
f.text(0.5, 0.04, "$mag$", ha="center", fontsize=15)
f.text(0.04, 0.5, r"$z_{phot} \; uncertainties$", va="center", rotation="vertical", fontsize=15)
for rm in range(len(self.range_z) - 1):
zmin = self.range_z[rm]
zmax = self.range_z[rm + 1]
ax = axarr[(rm // 2) % self.nbRowM, rm % 2]
mstep = 0.2
mstart = 16
mnbstep = 100
medpvec = []
mednvec = []
magvec = []
ax.axis([15, 27, -0.99, 1])
for im in range(mnbstep):
diffp = self.zmlu68 - self.zml
diffn = self.zmll68 - self.zml
mlimi = mstart + mstep * im
mlims = mstart + mstep * (im + 1)
conda = (
self.cond
& self.condgal
& (self.mag > mlimi)
& (self.mag <= mlims)
& (self.zml > zmin)
& (self.zml < zmax)
)
if len(self.zml[conda]) > 3:
medp = np.median(diffp[conda])
medn = np.median(diffn[conda])
medpvec.append(medp)
mednvec.append(medn)
magvec.append(mstart + mstep * (im + 0.5))
ax.plot(magvec, medpvec, color="b", linestyle="-")
ax.plot(magvec, mednvec, color="r", linestyle="-")
ax.annotate(f"${zmin:.2f} < z < {zmax:.2f}$", xy=(mstart, 0.60), color="black", fontsize=15)
return
[docs]
def errorz(self):
"""
Plot photometric redshift errors as a function of redshift.
This diagnostic shows how the estimated photometric redshift uncertainties
(e.g., derived from the PDF width or χ² fitting) vary with redshift.
It is useful for identifying redshift ranges where photometric redshifts
may be less reliable or exhibit larger scatter.
Notes
-----
- Uses the photometric redshifts (`Z_BEST` or `Z_MED`) and their
associated uncertainties from the LePHARE output.
- Only objects within the redshift and magnitude selection criteria are included.
- Errors often increase at high redshift due to fainter magnitudes and
fewer spectral features in the observed bands.
See Also
--------
errormag : Photometric redshift errors versus observed magnitude.
check_error : Evaluate 68% error coverage as a function of magnitude.
cumulative68 : Check fraction of spectroscopic redshifts within the 68% photometric error.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.errorz()
# Produces a plot of photometric redshift errors versus redshift.
"""
plt.clf()
f, axarr = plt.subplots(
self.nbRowM, self.nbColM, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
plt.subplots_adjust(hspace=0, wspace=0)
f.text(0.5, 0.04, "$z_{phot}$", ha="center", fontsize=15)
f.text(0.04, 0.5, r"$z_{phot} \; uncertainties$", va="center", rotation="vertical", fontsize=15)
for rm in range(len(self.range_mag) - 1):
ax = axarr[(rm // 2) % self.nbRowM, rm % 2]
magmin = self.range_mag[rm]
magmax = self.range_mag[rm + 1]
zstep = 0.1
zstart = 0
znbstep = 60
medpvec = []
mednvec = []
zvec = []
ax.axis([0, 5.9, -0.99, 1])
for iz in range(znbstep):
diffp = self.zmlu68 - self.zml
diffn = self.zmll68 - self.zml
zlimi = zstart + zstep * iz
zlims = zstart + zstep * (iz + 1)
conda = (
self.cond
& self.condgal
& (self.zml > zlimi)
& (self.zml <= zlims)
& (self.mag > magmin)
& (self.mag < magmax)
)
if len(self.zml[conda]) > 3:
medp = np.median(diffp[conda])
medn = np.median(diffn[conda])
medpvec.append(medp)
mednvec.append(medn)
zvec.append(zstart + zstep * (iz + 0.5))
ax.plot(zvec, medpvec, color="b", linestyle="-")
ax.plot(zvec, mednvec, color="r", linestyle="-")
ax.annotate(
f"${magmin:.2f} < mag < {magmax:.2f}$",
xy=(zstart + 0.2, 0.60),
color="black",
fontsize=15,
)
return
[docs]
def pit_qq(
self,
pdfs=None,
zgrid=None,
ztrue=None,
bins=None,
title=None,
show_pit=True,
show_qq=True,
savefig=False,
) -> str:
"""
Generate a PIT and QQ quantile-quantile plot for photometric redshift PDFs.
This function creates a visualization to assess the calibration of
predicted PDFs against true redshifts. It can display a Probability
Integral Transform (PIT) histogram and/or a QQ plot, optionally saving
the figure.
Parameters
----------
pdfs : np.ndarray, shape (m, n), optional
Array of PDFs for m objects over n z-grid points. If None, use the
PDFs from the class Metrics object.
zgrid : np.ndarray, shape (n,), optional
The z-axis corresponding to the PDF bins. Required if `pdfs` is provided.
ztrue : np.ndarray, shape (m,), optional
True redshift values for the m objects. Required if `pdfs` is provided.
bins : int, optional
Number of bins for the PIT histogram. If None, use the default
number of quantiles defined in the sample (`sample.n_quant`).
title : str, optional
Title for the plot. If None, a formatted sample name will be used.
show_pit : bool, default=True
Whether to include the PIT histogram in the figure.
show_qq : bool, default=True
Whether to include the QQ plot in the figure.
savefig : bool, default=False
Whether to save the plot to a PNG file. If True, the filename will
be automatically generated.
Returns
-------
str
Path to the saved figure if `savefig=True`; otherwise an empty string.
Notes
-----
- PIT (Probability Integral Transform) values should be uniformly
distributed if the PDFs are well-calibrated.
- QQ plot compares the quantiles of the PIT distribution against a
uniform distribution to visually assess calibration.
"""
if pdfs is None:
pdfs = self.pdfs
if zgrid is None:
zgrid = np.linspace(0, 6, pdfs.shape[1])
if ztrue is None:
ztrue = self.zs
if bins is None:
bins = 100
if title is None:
title = ""
plt.figure(figsize=[4, 5])
gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1])
ax0 = plt.subplot(gs[0])
# sample = Sample(pdfs, zgrid, ztrue)
pit_vals = integrate_pdfs_to_ztrue(pdfs, zgrid, ztrue)
pit_out_rate = np.mean((pit_vals < 0.005) | (pit_vals > 0.995))
qq = [np.linspace(0, 1, len(ztrue)), pit_vals[np.argsort(pit_vals)]]
if show_qq:
ax0.plot(
qq[0], qq[1], c="r", linestyle="-", linewidth=3, label=f"PIT$_{{out}}$:{pit_out_rate:.2f}"
)
ax0.plot([0, 1], [0, 1], color="k", linestyle="--", linewidth=2)
ax0.set_ylabel("Q$_{data}$", fontsize=18)
plt.ylim(-0.001, 1.001)
plt.xlim(-0.001, 1.001)
plt.title(title)
if show_pit:
# fzdata = Ensemble(interp, data=dict(xvals=zgrid, yvals=pdfs))
# pitobj = PIT(fzdata, ztrue)
# pit_vals = np.array(pitobj.pit_samps)
# pit_out_rate = pitobj.evaluate_PIT_outlier_rate()
try:
y_uni = float(len(pit_vals)) / float(bins)
except TypeError: # when bins is an array
y_uni = float(len(pit_vals)) / float(len(bins))
if not show_qq:
ax0.hist(pit_vals, bins=bins, alpha=0.7, label="LePHARE")
ax0.set_ylabel("Number")
ax0.hlines(y_uni, xmin=0, xmax=1, color="k")
plt.ylim(
0,
) # -0.001, 1.001)
else:
ax1 = ax0.twinx()
hist = ax1.hist(pit_vals, bins=bins, alpha=0.7)
ax1.set_ylim([0, 1.1 * np.max(hist[0][1:-1])])
ax1.set_ylabel("Number")
ax1.hlines(y_uni, xmin=0, xmax=1, color="k")
leg = ax0.legend(handlelength=0, handletextpad=0, fancybox=True, loc="upper left")
for item in leg.legend_handles:
item.set_visible(False)
if show_qq:
ax2 = plt.subplot(gs[1])
ax2.plot(
qq[0],
(qq[1] - qq[0]),
c="r",
linestyle="-",
linewidth=3,
)
plt.ylabel(r"$\Delta$Q", fontsize=18)
ax2.plot([0, 1], [0, 0], color="k", linestyle="--", linewidth=2)
plt.xlim(-0.001, 1.001)
plt.ylim(
np.min([-0.12, np.min(qq[1] - qq[0]) * 1.05]),
np.max([0.12, np.max(qq[1] - qq[0]) * 1.05]),
)
if show_pit:
if show_qq:
plt.xlabel("Q$_{theory}$ / PIT Value", fontsize=18)
else:
plt.xlabel("PIT Value", fontsize=18)
else:
if show_qq:
plt.xlabel("Q$_{theory}$", fontsize=18)
if savefig:
fig_filename = str("plot_pit_qq_" + "lephare.png")
plt.savefig(fig_filename)
else:
fig_filename = None
return fig_filename
[docs]
def dist_mass(self, nstep=10):
"""
Plot the distribution of stellar mass
This diagnostic shows the range and frequency of stellar mass values
derived by LePHARE for the sample.
Parameters
----------
nstep : int, optional
Number of bins to use in the histogram of mass values. Default is 10.
Notes
-----
- Uses the `MASS_BEST` and 'MASS_MED' column from the LePHARE output table.
- Only sources meeting the redshift and magnitude selection criteria are included.
See Also
--------
dist_model : Show distribution of best-fit SED model identifiers.
dist_filt : Show number of filters used in each fit.
dist_SFR : SFR distribution
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.dist_mass(nstep=15)
# Produces a histogram of MASS_BEST and MASS_MED values.
"""
## Create a figure with an array with nbRowZ*nbColZ subpanels
plt.clf()
f, axarr = plt.subplots(
self.nbRowZ, self.nbColZ, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
# No space between the figures
plt.subplots_adjust(hspace=0, wspace=0)
# label of the figure
f.text(0.5, 0.04, "log10(stellar mass)", ha="center", fontsize=15)
f.text(0.04, 0.5, "N", va="center", rotation="vertical", fontsize=15)
# Loop over the redshift bins
for rm in range(len(self.range_z) - 1):
# Define the subplots and pass the panel
zmin = self.range_z[rm]
zmax = self.range_z[rm + 1]
ax = axarr[int(rm // 2.0), int(rm % 2)]
# set the axis
ax.axis([6, 11.9, 0, 1.3])
# new condition within the redshift range
conda = self.cond & self.condgal & (self.zp > zmin) & (self.zp < zmax) & (self.massm > 0)
# check if some objects exist
if len(self.massm[conda]) > 1:
# Histogram of the mass
ax.hist(
self.massm[conda],
bins=nstep,
histtype="stepfilled",
density=1,
color="b",
label="median of the PDF",
)
# new condition within the redshift range
conda = self.cond & self.condgal & (self.zp > zmin) & (self.zp < zmax) & (self.massb > 0)
# check if some objects exist
if len(self.massb[conda]) > 1:
# Histogram of the mass
ax.hist(
self.massb[conda],
bins=nstep,
histtype="stepfilled",
density=1,
color="r",
alpha=0.5,
label=r"$minimum\; \chi^2$",
)
# labels
ax.annotate(
"$" + str(round(zmin, 2)) + " < z < " + str(round(zmax, 2)) + "$",
xy=(8, 1),
color="black",
fontsize=15,
)
# print the legend
if rm == 1:
ax.legend(prop={"size": 8})
return
[docs]
def dist_sfr(self, nstep=10):
"""
Plot the distribution of SFR
This diagnostic shows the range and frequency of SFR values
derived by LePHARE for the sample.
Parameters
----------
nstep : int, optional
Number of bins to use in the histogram of mass values. Default is 10.
Notes
-----
- Uses the `SFR_BEST` and 'SFR_MED' column from the LePHARE output table.
- Only sources meeting the redshift and magnitude selection criteria are included.
See Also
--------
dist_model : Show distribution of best-fit SED model identifiers.
dist_filt : Show number of filters used in each fit.
dist_mass : Mass distribution
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.dist_sfr(nstep=15)
# Produces a histogram of SFR_BEST and SFR_MED values.
"""
## Create a figure with an array with nbRowZ*nbColZ subpanels
plt.clf()
f, axarr = plt.subplots(
self.nbRowZ, self.nbColZ, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
# No space between the figures
plt.subplots_adjust(hspace=0, wspace=0)
# label of the figure
f.text(0.5, 0.04, "log10(SFR)", ha="center", fontsize=15)
f.text(0.04, 0.5, "N", va="center", rotation="vertical", fontsize=15)
# Loop over the redshift bins
for rm in range(len(self.range_z) - 1):
# Define the subplots and pass the panel
zmin = self.range_z[rm]
zmax = self.range_z[rm + 1]
ax = axarr[int(rm // 2.0), int(rm % 2)]
# set the axis
ax.axis([-2, 3, 0, 1.3])
# new condition within the redshift range
conda = self.cond & self.condgal & (self.zp > zmin) & (self.zp < zmax) & (self.sfrm > 0)
# check if some objects exist
if len(self.sfrm[conda]) > 1:
# Histogram of the SFR
ax.hist(
self.sfrm[conda],
bins=nstep,
histtype="stepfilled",
density=1,
color="b",
label="median of the PDF",
)
# new condition within the redshift range
conda = self.cond & self.condgal & (self.zp > zmin) & (self.zp < zmax) & (self.sfrb > 0)
# check if some objects exist
if len(self.sfrb[conda]) > 1:
# Histogram of the SFR
ax.hist(
self.sfrb[conda],
bins=nstep,
histtype="stepfilled",
density=1,
color="r",
alpha=0.5,
label=r"$minimum\; \chi^2$",
)
# labels
ax.annotate(
"$" + str(round(zmin, 2)) + " < z < " + str(round(zmax, 2)) + "$",
xy=(-1, 1),
color="black",
fontsize=15,
)
# print the legend
if rm == 1:
ax.legend(prop={"size": 8})
return
[docs]
def dist_ssfr(self, nstep=10):
"""
Plot the distribution of specific SFR
This diagnostic shows the range and frequency of sSFR values
derived by LePHARE for the sample.
Parameters
----------
nstep : int, optional
Number of bins to use in the histogram of mass values. Default is 10.
Notes
-----
- Uses the `SSFR_BEST` and 'SSFR_MED' column from the LePHARE output table.
- Only sources meeting the redshift and magnitude selection criteria are included.
See Also
--------
dist_model : Show distribution of best-fit SED model identifiers.
dist_filt : Show number of filters used in each fit.
dist_mass : Mass distribution
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.dist_ssfr(nstep=15)
# Produces a histogram of SSFR_BEST and SSFR_MED values.
"""
## Create a figure with an array with nbRowZ*nbColZ subpanels
plt.clf()
f, axarr = plt.subplots(
self.nbRowZ, self.nbColZ, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
# No space between the figures
plt.subplots_adjust(hspace=0, wspace=0)
# label of the figure
f.text(0.5, 0.04, "log10(sSFR)", ha="center", fontsize=15)
f.text(0.04, 0.5, "N", va="center", rotation="vertical", fontsize=15)
# Loop over the redshift bins
for rm in range(len(self.range_z) - 1):
# Define the subplots and pass the panel
zmin = self.range_z[rm]
zmax = self.range_z[rm + 1]
ax = axarr[int(rm // 2.0), int(rm % 2)]
# set the axis
ax.axis([-14, -7, 0, 1.3])
# new condition within the redshift range
conda = (
self.cond
& self.condgal
& (self.zp > zmin)
& (self.zp < zmax)
& (self.ssfrm > -20)
& (self.ssfrm < -6)
)
# check if some objects exist
if len(self.ssfrm[conda]) > 1:
# Histogram of the SFR
ax.hist(
self.ssfrm[conda],
bins=nstep,
histtype="stepfilled",
density=1,
color="b",
label="median of the PDF",
)
# new condition within the redshift range
conda = (
self.cond
& self.condgal
& (self.zp > zmin)
& (self.zp < zmax)
& (self.ssfrb > -20)
& (self.ssfrb < -6)
)
# check if some objects exist
if len(self.ssfrb[conda]) > 1:
# Histogram of the SFR
ax.hist(
self.ssfrb[conda],
bins=nstep,
histtype="stepfilled",
density=1,
color="r",
alpha=0.5,
label=r"$minimum\; \chi^2$",
)
# labels
ax.annotate(
"$" + str(round(zmin, 2)) + " < z < " + str(round(zmax, 2)) + "$",
xy=(-13.5, 1),
color="black",
fontsize=15,
)
# print the legend
if rm == 1:
ax.legend(prop={"size": 8})
return
[docs]
def mass_med_best(self):
"""
Plot stellar mass median of the PDF versus minimum χ²
This diagnostic compares the two mass estimates produced by
LePHARE: the median of the posterior probability distribution (`MASS_MED`)
and the best-fit value from the minimum χ² solution (`MASS_BEST`).
The comparison helps identify systematic biases between the two estimators
and assess the stability of mass determinations.
The figure is divided into subpanels based on magnitude bins defined in
`range_mag`. Each subpanel shows `MASS_MED` versus `MASS_BEST` with a 1:1
reference line and scatter indicating consistency between the two measures.
Notes
-----
- Only sources within the defined redshift and magnitude ranges are included.
- The plot can reveal populations where the PDF and best-fit solutions diverge,
such as multi-peaked PDFs or poor fits.
- Strong deviations from the 1:1 line indicate inconsistency between the
χ²-based and PDF-based estimates.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.mass_med_best()
# Produces a MASS_MED vs MASS_BEST comparison plot by magnitude bin.
"""
plt.clf()
f, axarr = plt.subplots(
self.nbRowM, self.nbColM, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
# No space between the figures
plt.subplots_adjust(hspace=0, wspace=0)
# label of the figure
f.text(0.5, 0.04, r"$mass \; minimum\; \chi^2$", ha="center", fontsize=15)
f.text(0.04, 0.5, r"$mass \; median\; PDF(z)$", va="center", rotation="vertical", fontsize=15)
# Loop over the magnitude bins
for rm in range(len(self.range_mag) - 1):
# Define the subplots and pass the panel
ax = axarr[int(rm // 2.0), int(rm % 2)]
# mag limits
magmin = self.range_mag[rm]
magmax = self.range_mag[rm + 1]
# set the axis
ax.axis([6, 11.9, 6, 11.9])
# new condition with the magnitude range
conda = self.cond & self.condgal & (self.mag > magmin) & (self.mag < magmax)
# Plot mass med versus best
ax.scatter(self.massb[conda], self.massm[conda], s=1, color="b", alpha=0.5, marker="s")
# Trace the limits 0.15(1+z)
x_zs = np.array([6, 15])
ax.plot(x_zs, x_zs, "r-")
# labels
ax.annotate(
"$" + str(round(magmin, 2)) + " < mag < " + str(round(magmax, 2)) + "$",
xy=(7, 11),
color="black",
fontsize=15,
)
return
[docs]
def sfr_med_best(self):
"""
Plot SFR median of the PDF versus minimum χ²
This diagnostic compares the two mass estimates produced by
LePHARE: the median of the posterior probability distribution (`SFR_MED`)
and the best-fit value from the minimum χ² solution (`SFR_BEST`).
The comparison helps identify systematic biases between the two estimators
and assess the stability of mass determinations.
The figure is divided into subpanels based on magnitude bins defined in
`range_mag`. Each subpanel shows `SFR_MED` versus `SFR_BEST` with a 1:1
reference line and scatter indicating consistency between the two measures.
Notes
-----
- Only sources within the defined redshift and magnitude ranges are included.
- The plot can reveal populations where the PDF and best-fit solutions diverge,
such as multi-peaked PDFs or poor fits.
- Strong deviations from the 1:1 line indicate inconsistency between the
χ²-based and PDF-based estimates.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.sfr_med_best()
# Produces a SFR_MED vs SFR_BEST comparison plot by magnitude bin.
"""
plt.clf()
f, axarr = plt.subplots(
self.nbRowM, self.nbColM, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
# No space between the figures
plt.subplots_adjust(hspace=0, wspace=0)
# label of the figure
f.text(0.5, 0.04, r"$SFR \; minimum\; \chi^2$", ha="center", fontsize=15)
f.text(0.04, 0.5, r"$SFR \; median\; PDF(z)$", va="center", rotation="vertical", fontsize=15)
# Loop over the magnitude bins
for rm in range(len(self.range_mag) - 1):
# Define the subplots and pass the panel
ax = axarr[int(rm // 2.0), int(rm % 2)]
# mag limits
magmin = self.range_mag[rm]
magmax = self.range_mag[rm + 1]
# set the axis
ax.axis([-2, 3.9, -2, 3.9])
# new condition with the magnitude range
conda = self.cond & self.condgal & (self.mag > magmin) & (self.mag < magmax)
# Plot SFR med versus best
ax.scatter(self.sfrb[conda], self.sfrm[conda], s=1, color="b", alpha=0.5, marker="s")
# Trace the limits 0.15(1+z)
x_zs = np.array([-6, 15])
ax.plot(x_zs, x_zs, "r-")
# labels
ax.annotate(
"$" + str(round(magmin, 2)) + " < mag < " + str(round(magmax, 2)) + "$",
xy=(-1, 3),
color="black",
fontsize=15,
)
return
[docs]
def mass_z(self):
"""
Plot stellar mass versus redshift for the sample.
This diagnostic shows the relation between stellar masses
and redshift, providing insight into mass evolution, selection effects, and completeness
limits of the survey.
Notes
-----
- Uses the stellar mass
- Only objects within the redshift and magnitude selection ranges are included.
- The figure can be divided into subpanels based on magnitude or redshift bins.
- Helps identify trends such as brightening or dimming with redshift
and the distribution of galaxies in the mass-redshift plane.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.mass_z()
# Produces plots of mass versus redshift
"""
plt.clf()
f, axarr = plt.subplots(
self.nbRowM, self.nbColM, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
plt.subplots_adjust(hspace=0, wspace=0)
f.text(0.5, 0.04, "$redshift$", ha="center", fontsize=15)
f.text(0.04, 0.5, "$log10(mass)$", va="center", rotation="vertical", fontsize=15)
for rm in range(len(self.range_mag) - 1):
ax = axarr[int(rm // 2.0), int(rm % 2)]
magmin = self.range_mag[rm]
magmax = self.range_mag[rm + 1]
ax.axis([self.z_min, self.z_max, 7, 11.9])
conda = self.cond & self.condgal & (self.mag > magmin) & (self.mag < magmax)
ax.scatter(self.zp[conda], self.massb[conda], s=1, color="b", alpha=0.2, marker="s")
ax.annotate(f"${magmin:.2f} < mag < {magmax:.2f}$", xy=(0.2, 11.5), color="black", fontsize=15)
return
[docs]
def sfr_z(self):
"""
Plot SFR versus redshift for the sample.
This diagnostic shows the relation between SFR
and redshift, providing insight into mass evolution, selection effects, and completeness
limits of the survey.
Notes
-----
- Uses the SFR
- Only objects within the redshift and magnitude selection ranges are included.
- The figure can be divided into subpanels based on magnitude or redshift bins.
- Helps identify trends such as brightening or dimming with redshift
and the distribution of galaxies in the mass-redshift plane.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.sfr_z()
# Produces plots of mass versus redshift
"""
plt.clf()
f, axarr = plt.subplots(
self.nbRowM, self.nbColM, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
plt.subplots_adjust(hspace=0, wspace=0)
f.text(0.5, 0.04, "$redshift$", ha="center", fontsize=15)
f.text(0.04, 0.5, "$log10(SFR)$", va="center", rotation="vertical", fontsize=15)
for rm in range(len(self.range_mag) - 1):
ax = axarr[int(rm // 2.0), int(rm % 2)]
magmin = self.range_mag[rm]
magmax = self.range_mag[rm + 1]
ax.axis([self.z_min, self.z_max, -3, 3.9])
conda = self.cond & self.condgal & (self.mag > magmin) & (self.mag < magmax)
ax.scatter(self.zp[conda], self.sfrb[conda], s=1, color="b", alpha=0.2, marker="s")
ax.annotate(f"${magmin:.2f} < mag < {magmax:.2f}$", xy=(0.2, 3.0), color="black", fontsize=15)
return
[docs]
def mass_sfr(self):
"""
Plot stellar mass versus SFR for the sample.
This diagnostic shows the main sequence
Notes
-----
- Uses the stellar mass and SFR
- Only objects within the redshift and magnitude selection ranges are included.
- The figure can be divided into subpanels based on magnitude or redshift bins.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.mass_sfr()
# Produces plots of mass versus SFR
"""
plt.clf()
f, axarr = plt.subplots(
self.nbRowM, self.nbColM, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
plt.subplots_adjust(hspace=0, wspace=0)
f.text(0.5, 0.04, "$log10(mass)$", ha="center", fontsize=15)
f.text(0.04, 0.5, "$log10(SFR)$", va="center", rotation="vertical", fontsize=15)
for rm in range(len(self.range_z) - 1):
ax = axarr[int(rm // 2.0), int(rm % 2)]
zmin = self.range_z[rm]
zmax = self.range_z[rm + 1]
ax.axis([7, 12, -3, 4])
# new condition with the redshift range and star-forming
conda = self.cond & (self.zp > zmin) & (self.zp < zmax) & (self.ssfrb > -11) & (self.ssfrb < 90)
ax.scatter(self.massb[conda], self.sfrb[conda], s=1, color="b", alpha=0.2, marker="s")
# new condition with the redshift range and quiescent
conda = self.cond & (self.zp > zmin) & (self.zp < zmax) & (self.ssfrb < -11)
ax.scatter(self.massb[conda], self.sfrb[conda], s=1, color="r", alpha=0.2, marker="s")
ax.annotate(f"${zmin:.2f} < z < {zmax:.2f}$", xy=(8, 3), color="black", fontsize=15)
return
[docs]
def lnuv_sfr(self):
"""
Plot dust-corrected L(NUV) versus SFR for the sample.
Notes
-----
- Uses the Lnuv and SFR
- Only objects within the redshift and magnitude selection ranges are included.
- The figure can be divided into subpanels based on magnitude or redshift bins.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.lnuv_sfr()
# Produces plots of L(NUV) versus SFR
"""
plt.clf()
f, axarr = plt.subplots(
self.nbRowM, self.nbColM, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
plt.subplots_adjust(hspace=0, wspace=0)
f.text(0.5, 0.04, "$log10(SFR)$", ha="center", fontsize=15)
f.text(0.04, 0.5, "$L(NUV)$", va="center", rotation="vertical", fontsize=15)
for rm in range(len(self.range_z) - 1):
ax = axarr[int(rm // 2.0), int(rm % 2)]
zmin = self.range_z[rm]
zmax = self.range_z[rm + 1]
ax.axis([-4, 3, 6, 11.9])
# new condition with the redshift range and star-forming
conda = self.cond & (self.zp > zmin) & (self.zp < zmax) & (self.sfrb > -90)
ax.scatter(self.sfrb[conda], self.lnuv[conda], s=1, color="b", alpha=0.2, marker="s")
ax.annotate(f"${zmin:.2f} < z < {zmax:.2f}$", xy=(-3, 10), color="black", fontsize=15)
return
[docs]
def masstolight_z(self):
"""
Plot mass-to-light ratio versus redshift for the sample.
This diagnostic shows the relation between mass-to-light ratio
and redshift, providing insight into SED fitting
Notes
-----
- Uses the L(K) and massb
- Only objects within the redshift and magnitude selection ranges are included.
- The figure can be divided into subpanels based on magnitude or redshift bins.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.masstolight_z()
# Produces plots of mass-to-light versus redshift
"""
plt.clf()
f, axarr = plt.subplots(
self.nbRowM, self.nbColM, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
plt.subplots_adjust(hspace=0, wspace=0)
f.text(0.5, 0.04, "$redshift$", ha="center", fontsize=15)
f.text(0.04, 0.5, r"$log10(M/L \, ratio)$", va="center", rotation="vertical", fontsize=15)
# Define the mass-to-light ratio
mlratio = (
self.massm
- self.lk
+ (0.4 * (51.605 - 5.14))
+ log10(3.0e18 * 2000.0 / pow(22000.0, 2.0))
- log10(3.826e33)
)
for rm in range(len(self.range_mag) - 1):
ax = axarr[int(rm // 2.0), int(rm % 2)]
magmin = self.range_mag[rm]
magmax = self.range_mag[rm + 1]
ax.axis([self.z_min, self.z_max, -1.4, 0.2])
conda = self.cond & self.condgal & (self.mag > magmin) & (self.mag < magmax)
ax.scatter(self.zp[conda], mlratio[conda], s=1, color="b", alpha=0.2, marker="s")
ax.annotate(f"${magmin:.2f} < mag < {magmax:.2f}$", xy=(0.2, 0), color="black", fontsize=15)
# Trace the limits
x = np.array([0, 6])
ax.plot(x, -0.27 * x - 0.05 - 0.24, "r-")
ax.plot(x, -0.18 * x - 0.05 - 0.24, "r-")
return
[docs]
def absmagk_mass(self):
"""
Plot M(Ks) versus mass for the sample.
Notes
-----
- Uses the M(Ks) and mass
- Only objects within the redshift and magnitude selection ranges are included.
- The figure can be divided into subpanels based on magnitude or redshift bins.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.absmagk_mass()
# Produces plots of M(K) versus mass
"""
plt.clf()
f, axarr = plt.subplots(
self.nbRowM, self.nbColM, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
plt.subplots_adjust(hspace=0, wspace=0)
f.text(0.5, 0.04, "$log10(mass)$", ha="center", fontsize=15)
f.text(0.04, 0.5, "$M(K)$", va="center", rotation="vertical", fontsize=15)
for rm in range(len(self.range_z) - 1):
ax = axarr[int(rm // 2.0), int(rm % 2)]
zmin = self.range_z[rm]
zmax = self.range_z[rm + 1]
ax.axis([6, 11.9, -24.9, -13])
# new condition with the magnitude range and star-forming
conda = self.cond & (self.zp > zmin) & (self.zp < zmax) & (self.ssfrm > -11) & (self.ssfrm < 0)
if len(self.massb[conda]) > 0:
ax.scatter(self.massb[conda], self.mabsk[conda], s=1, color="b", alpha=0.2, marker="s")
# new condition with the magnitude range and quiescent
conda = self.cond & (self.zp > zmin) & (self.zp < zmax) & (self.ssfrm < -11)
if len(self.massb[conda]) > 0:
ax.scatter(self.massb[conda], self.mabsk[conda], s=1, color="r", alpha=0.2, marker="s")
ax.annotate(f"${zmin:.2f} < z < {zmax:.2f}$", xy=(6.5, -23), color="black", fontsize=15)
return
[docs]
def absmagu_sfr(self):
"""
Plot M(U) versus SFR for the sample.
Notes
-----
- Uses the M(U) and SFR
- Only objects within the redshift and magnitude selection ranges are included.
- The figure can be divided into subpanels based on magnitude or redshift bins.
Examples
--------
>>> utils = lp.PlotUtils(t, sel_filt=3)
>>> utils.absmagu_sfr()
# Produces plots of M(U) versus SFR
"""
plt.clf()
f, axarr = plt.subplots(
self.nbRowM, self.nbColM, sharex=True, sharey=True, figsize=(12, 8), squeeze=False
)
plt.subplots_adjust(hspace=0, wspace=0)
f.text(0.5, 0.04, "$log10(SFR)$", ha="center", fontsize=15)
f.text(0.04, 0.5, "$M(U)$", va="center", rotation="vertical", fontsize=15)
for rm in range(len(self.range_z) - 1):
ax = axarr[int(rm // 2.0), int(rm % 2)]
zmin = self.range_z[rm]
zmax = self.range_z[rm + 1]
ax.axis([-4, 4, -24.9, -13])
# new condition with the redshift range and star-forming
conda = self.cond & (self.zp > zmin) & (self.zp < zmax) & (self.sfrb > -90)
condb = conda & (self.ebv <= 0.1)
ax.scatter(self.sfrb[condb], self.mabsu[condb], s=1, color="b", alpha=0.2, marker="s")
condb = conda & (self.ebv > 0.1) & (self.ebv < 0.3)
ax.scatter(self.sfrb[condb], self.mabsu[condb], s=1, color="g", alpha=0.2, marker="s")
condb = conda & (self.ebv >= 0.3)
ax.scatter(self.sfrb[condb], self.mabsu[condb], s=1, color="r", alpha=0.2, marker="s")
ax.annotate(f"${zmin:.2f} < z < {zmax:.2f}$", xy=(-3, -23), color="black", fontsize=15)
return
def integrate_pdfs_to_ztrue(pdfs, zgrid, ztrue):
"""
Integrate each PDF up to its corresponding true z value.
Parameters
----------
pdfs : np.ndarray, shape (m, n)
PDFs for m objects over n bins.
zgrid : np.ndarray, shape (n,)
Grid for the PDF bins.
ztrue : np.ndarray, shape (m,)
Value up to which to integrate each PDF.
Returns
-------
integrated : np.ndarray, shape (m,)
Integrated value of each PDF up to the corresponding ztrue.
"""
m, n = pdfs.shape
integrated = np.zeros(m)
for i in range(m):
# Interpolate PDF for this object
f = interp1d(zgrid, pdfs[i], kind="linear", bounds_error=False, fill_value=0.0)
# Determine the integration points: all zgrid points <= ztrue[i]
mask = zgrid <= ztrue[i]
if np.any(mask):
zvals = zgrid[mask]
pdfvals = f(zvals)
integrated[i] = trapezoid(pdfvals, zvals)
else:
integrated[i] = 0.0 # if ztrue[i] < zgrid[0]
return integrated