A Completely Unauthorized Users Guide
to EAZY-py

07 Understanding EAZY's output: z_ml, z_best, z_phot, and z_min_chisq

This is one of the most confusing aspects of EAZY-py — there are multiple "best redshift" concepts that are closely related but subtly different. Here are the output values as I understand them.

zml (Maximum Likelihood Redshift)

self.zml is the primary internal redshift estimate, and it's defined as the redshift at the maximum of lnp. However, there's some nuance here, as it is computed by eazy.evaluate_zml(), which:

1. Finds iz_best = argmax(lnp[i, :]) for each object
2. Fits a parabola around that peak for sub-grid-spacing precision
3. Returns the continuous redshift at the parabola peak as eazy.zml

It incorporates whatever priors were active when compute_lnp() was called:

self.compute_lnp(prior=True, beta_prior=True)
self.evaluate_zml(prior=True, beta_prior=True)
# Now self.zml reflects the prior-weighted maximum likelihood

z_best (zbest)

self.zbest is the redshift used for computing physical parameters (rest-frame colors, masses, etc.). It defaults to self.zml, but can be set to any user-specified value. If FIX_ZSPEC = True, zbest is set to z_spec for objects where it's available.

zout['z_phot']

In the output table, z_phot is the maximum-likelihood redshift used to compute everything in that table, effectively zml with whatever prior settings were passed to standard_output().

zout['z_raw_chi2']

This is the redshift that corresponds to the overall minimum χ² value across all of the redshifts in the redshift grid.

zout['z500'] (and z025, z160, z840, z975)

This redshift corresponding to the 50th percentile in the P(z) surface, derived from the output χ²(z), and a prior, if one is used. The other values are the 2.5, 16, 84, and 97.5th percentiles, for understanding uncertainty. You can specify these values with percentile_limits=[2.5, 16, 50, 84, 97.5] in self.standard_output().

Summary

09 Priors: Magnitude Prior & Beta Prior

The Apparent Magnitude Prior

The magnitude prior encodes the expectation that fainter sources are more likely to be at higher redshift. It is read from a file (default templates/prior_F160W_TAO.dat) that contains P(m, z) — the probability of a source being at redshift z given its apparent magnitude m.

# Recompute without prior to diagnose bias
ez.compute_lnp(prior=False, beta_prior=False)
ez.evaluate_zml(prior=False, beta_prior=False)

The Beta Prior

The UV slope beta prior (beta_prior=True) adds a soft constraint on the UV spectral slope, discouraging fits where the galaxy would need an unusual UV slope at the fitted redshift. The prior width is a sigmoid function of redshift:

• At z < z_split (default 4): tight, σ = 0.5
• At z > z_split: widens significantly toward σ = 20

Particularly useful for JWST high-z work where Lyman-break galaxies tend to have very blue UV slopes, and the beta prior helps avoid low-z degeneracies.

# Customize the beta prior shape
ez.prior_beta(width_params={
    'k': -5, 'z_split': 4,
    'sigma0': 20, 'sigma1': 0.5, 'center': -1.5
})

10 The Template Error Function (TEF)

The TEF addresses a fundamental problem: even perfect photometry gives poor χ² if the templates don't perfectly match your galaxies. The TEF adds wavelength-dependent "template uncertainty" to the error budget:

The TEF is read from TEMP_ERR_FILE and scaled by TEMP_ERR_A2. Increasing TEMP_ERR_A2 broadens all P(z) distributions, and could be used when you don't fully trust your templates.

11 Running standard_output: What You Get

self.standard_output() is the main output-generation function. It:

1. Refits template coefficients at zbest (or z_ml)
2. Computes redshift statistics (percentiles of P(z), risk function, etc.)
3. Computes rest-frame photometry in UBVJ and any other specified filters
4. Optionally computes SPS-derived physical parameters (stellar mass, SFR, Lv, Av)
5. Saves {MAIN_OUTPUT_FILE}.zout.fits and {MAIN_OUTPUT_FILE}.data.fits

zout, hdu = ez.standard_output(
    zbest=None,           # Use z_ml by default; pass array for custom redshifts
    prior=False,          # Apply magnitude prior
    beta_prior=False,     # Apply beta prior
    UBVJ=[153, 154, 155, 161],                   # Filter numbers for rest-frame U, B, V, J
    extra_rf_filters=[270, 274, 120, 121, ...],  # In case you want to calculate rest-frame photometry through these filters
    rf_pad_width=0.5,     # Padding (fraction of filter width) for rest-frame interp
    rf_max_err=2,         # Max allowed fractional error on rest-frame flux
    percentile_limits=[2.5, 16, 50, 84, 97.5],   # For z500, z160, etc.
    get_err=True,         # Compute rest-frame flux uncertainties via random draws
    save_fits=True,
    n_proc=0              # Parallelization for rest-frame computation
)

To reload products later:

from astropy.table import Table
zout = Table.read('photz.zout.fits')

# Or reload the full PhotoZ state
self = photoz.PhotoZ(params=params, load_products=True)
# self.zout, self.chi2_fit, self.fit_coeffs, self.ubvj all populated

NOTE: You will need to specify trhe rest-frame UBVJ filters, or else the code fails, so if you are defining your own filter file, either point these to arbitrary filter numbers (that exist!) or include the UBVJ filters from FILTER.RES.latest. I tend to just add on filters to the end of this file, as shown above.

12 The zout Table: Column-by-Column

Computing UVJ Colors

uv = -2.5 * np.log10(zout['restU'] / zout['restV'])
vj = -2.5 * np.log10(zout['restV'] / zout['restJ'])

uverr = 2.5 * np.sqrt((zout['restU_err']/zout['restU'])**2 +
                       (zout['restV_err']/zout['restV'])**2)
vjerr = 2.5 * np.sqrt((zout['restV_err']/zout['restV'])**2 +
                       (zout['restJ_err']/zout['restJ'])**2)

13 Iterative Zeropoint Corrections

There may be underlying differences in the fluxes between the photometric datapoints that can be ascertained with EAZY fits. The function iterate_zp_templates() iteratively adjusts per-filter flux zeropoints to minimize the template residuals as compared to the observed flues. After running, self.zp will contain the multiplicative correction factors per filter (close to 1.0 for well-calibrated data), which are normalized to the params['PRIOR_FILTER'] filter.

NITER = 3
NBIN  = min(self.NOBJ // 100, 180)

for i in range(NITER):
    sn   = self.fnu / self.efnu
    clip = (sn > 3).sum(axis=1) > 4   # objects with good multi-band SNR, you can adjust this to your liking

    self.iterate_zp_templates(
        idx=self.idx[clip],
        update_templates=False,
        update_zeropoints=True,
        iter=i,
        n_proc=4,
        error_residuals=(i > 0),   # rescale errors after first iteration
        NBIN=NBIN,
        get_spatial_offset=False
    )

The error_residuals=True option rescales per-filter error arrays so that the normalized χ² per filter is ~1, effectively correcting for over/underestimated uncertainties.

14 Template Sets: Which One to Use

Examine the templates after initialization:

import matplotlib.pyplot as plt
for i, temp in enumerate(self.templates):
    plt.plot(temp.wave / 1e4, temp.flux[0,:], label=temp.name, alpha=0.7)
plt.xlabel('Wavelength (µm)')
plt.xscale('log')
plt.yscale('log')
plt.legend(fontsize=6)

15 Visualizing Results

Individual SED Fits

# By object ID (from catalog)
fig, data = self.show_fit(
    id=12345,
    xlim=[0.3, 9],           # Wavelength range in µm
    show_components=True,    # Show individual template contributions
    logpz=True,              # Log scale for P(z)
    zr=[0, 7],               # Redshift range for P(z) panel
    show_prior=True,         # Overlay the magnitude prior
    show_stars=True,         # Overlay best-fit stellar templates
    delta_chi2_stars=-20,    # χ² threshold for flagging stars
    snr_thresh=2.0,          # SNR threshold for detections vs. limits
    showpz=0.6               # Fraction of figure height for P(z) panel
)

# By index rather than ID
fig, data = self.show_fit(id=42, id_is_idx=True)

# At a specific redshift (e.g., z_spec)
fig, data = self.show_fit(id=12345, zshow=1.234)

The data return value is a dict containing the SED, template fluxes, and P(z) arrays — useful for customizing the plot.

plot_fit_EAZY()

While the built-in function self.show_fit() works well, here is a python function I wrote which will allow you to plot the fit for individual sources with the minimum χ² value as well as the χ² surface, while also allowing the user to plot alternate redshift solutions and save the plot, or the SED fit and/or χ²(z) values.

Primarily I am including this to show the user how to go and grab information from self and zout, and specifically the power of self.show_fit() to extract information without making a plot.

        import numpy as np
        import matplotlib.pyplot as plt
        
        def plot_fit_EAZY(object_ID, self, zout, primary_z = -9999, alt_z = 0.0, h5file = False, savefile = None, savesed = False, savechisq = False):
	
	template_color = '#56B4E9'#'
	NIRC_photometry_color = '#D55E00'
	HST_photometry_color =  '#F786AA'
	chisq_surface_color = '#E69F00'
	zspec_color = '#332288'
	alternate_color = 'grey'

	# First, the redshift grid for the chisq surface. 
	output_zgrid = self.zgrid
	
	# the index of the source. 
	i = np.where(self.OBJID == object_ID)[0][0]

	# the chisq values. 
	if (h5file == False):
		output_chisq = self.chi2_fit[i,:]
	else:
		_data = self.get_object_data(i)
		z, fnu_i, efnu_i, ra_i, dec_i, output_chisq, zspec_i, ok_i = _data
	
	z_a_value = np.round(zout['z_raw_chi2'][i],2)
	
	if (primary_z < 0):
		data = self.show_fit(object_ID, show_fnu = True, zshow = z_a_value, get_spec = True)
	else:
		data = self.show_fit(object_ID, show_fnu = True, zshow = primary_z, get_spec = True)

	output_wavelength = data['templz']
	output_flux = data['templf']*1e3

	output_phot_wavelength = data['pivot']
	output_phot_template = data['model']*1e3

	output_phot_err_template = data['emodel']*1e3

	output_phot = data['fobs']*1e3

	output_phot_err = data['efobs']*1e3

	output_tef = data['tef']

	pos_flux_errors = np.where(output_phot_err > 0)[0]

	plt.figure(figsize = (12, 5))
	plt.subplot(121)
	if (primary_z < 0):
		plt.plot(output_wavelength/1e4, output_flux, color = template_color, lw = 3, zorder = 0, label = 'z = '+str(round(zout['z_raw_chi2'][i],2)))
	else:
		plt.plot(output_wavelength/1e4, output_flux, color = template_color, lw = 3, zorder = 0, label = 'z = '+str(round(primary_z,2)))
	
	plt.scatter(output_phot_wavelength/1e4, output_phot_template, marker = 's', s = 110, linewidth = 3, edgecolor = template_color, color = 'None', zorder = 5)
	plt.scatter(output_phot_wavelength[pos_flux_errors]/1e4, output_phot[pos_flux_errors], color = NIRC_photometry_color, s = 50, linewidth = 3, zorder = 10)
	plt.errorbar(output_phot_wavelength[pos_flux_errors]/1e4, output_phot[pos_flux_errors], yerr = output_phot_err[pos_flux_errors], color = 'black', ls = 'None', alpha = 0.4)
	plt.semilogy()
	plt.xlabel('Observed Wavelength ($\mu$m)')
	plt.ylabel('Flux (nJy)')
	plt.xlim(0.3, 5.0)
	pos_fluxes = np.where(output_phot > 0)[0]
	ymin = np.min(output_phot[pos_fluxes]) - (0.4 * np.min(output_phot[pos_fluxes]))
	if (ymin < 0.01):
		ymin = 0.01
	ymax = np.max(output_phot[pos_fluxes]) + (10.0 * np.max(output_phot[pos_fluxes]))
	if (ymax > (100.*np.median(output_phot[np.where(output_phot > 0)[0]]))):
		ymax = 100.*np.median(output_phot[np.where(output_phot > 0)[0]])
	print(ymin, ymax)
	plt.ylim(ymin, ymax)
	plt.title(self.OBJID[i])

	if (alt_z < 0):
		alt_z_value = alt_z * -1.0
		if (zout['z_raw_chi2'][i] > alt_z_value):
			zgrid_lt_7 = np.where(output_zgrid <= alt_z_value)[0]
			z_a_lt_7 = np.round(output_zgrid[zgrid_lt_7][np.argmin(output_chisq[zgrid_lt_7])],2)
			data_zlt_7 = self.show_fit(object_ID, show_fnu = True, zshow = z_a_lt_7, get_spec = True)

			delta_chisq = output_chisq[np.argmin(output_chisq[zgrid_lt_7])] - output_chisq[np.argmin(output_chisq)]
			print('Delta chisq = '+str(np.round(delta_chisq,2)))

			output_wavelength_zlt_7 = data_zlt_7['templz']
			output_flux_zlt_7 = data_zlt_7['templf']*1e3
			output_phot_template_zlt_7 = data_zlt_7['model']*1e3
		
			plt.plot(output_wavelength_zlt_7/1e4, output_flux_zlt_7, color = alternate_color, lw = 2, zorder = 0, label = 'z = '+str(round(z_a_lt_7,2)), alpha = 0.4)
		else:
			zgrid_gt_7 = np.where(output_zgrid > alt_z_value)[0]
			z_a_gt_7 = np.round(output_zgrid[zgrid_gt_7][np.argmin(output_chisq[zgrid_gt_7])],2)
			data_zgt_7 = self.show_fit(object_ID, show_fnu = True, zshow = z_a_gt_7, get_spec = True)

			delta_chisq = output_chisq[np.argmin(output_chisq[zgrid_gt_7])] - output_chisq[np.argmin(output_chisq)]
			print('Delta chisq = -'+str(np.round(delta_chisq,2)))

			output_wavelength_zgt_7 = data_zgt_7['templz']
			output_flux_zgt_7 = data_zgt_7['templf']*1e3
			output_phot_template_zgt_7 = data_zgt_7['model']*1e3
		
			plt.plot(output_wavelength_zgt_7/1e4, output_flux_zgt_7, color = alternate_color, lw = 2, zorder = 0, label = 'z = '+str(round(z_a_gt_7,2)), alpha = 0.4)
	if (alt_z > 0):
		alt_z_value = alt_z
		data_zalt = self.show_fit(object_ID, show_fnu = True, zshow = alt_z_value, get_spec = True)

		delta_chisq = output_chisq[np.where(output_zgrid == np.round(alt_z_value,2))[0]] - output_chisq[np.argmin(output_chisq)]
		print('Delta chisq = '+str(np.round(delta_chisq,2)))

		output_wavelength_zalt = data_zalt['templz']
		output_flux_zalt = data_zalt['templf']*1e3
		output_phot_template_zalt = data_zalt['model']*1e3
	
		plt.plot(output_wavelength_zalt/1e4, output_flux_zalt, color = alternate_color, lw = 2, zorder = 0, label = 'z = '+str(round(alt_z_value,2)), alpha = 0.4)
		

	plt.legend()

	plt.subplot(122)
	plt.plot(output_zgrid, output_chisq, color = chisq_surface_color, lw = 3, zorder = 10)
	if (primary_z < 0):
		plt.axvline(zout['z_raw_chi2'][i], lw = 3, zorder = 0)
	else:
		plt.axvline(primary_z, lw = 3, zorder = 0)

	if (alt_z < 0):
		if (zout['z_raw_chi2'][i] > 7.0):
			plt.axvline(z_a_lt_7, color = alternate_color, lw = 3, alpha = 0.4, zorder = 0)
		else:
			plt.axvline(z_a_gt_7, color = alternate_color, lw = 3, alpha = 0.4, zorder = 0)
	if (alt_z > 0):
		plt.axvline(alt_z_value, color = alternate_color, lw = 3, alpha = 0.4, zorder = 0)
		
	plt.xlabel('Redshift')
	plt.ylabel('$\chi^2$')

	plt.tight_layout()
	if (savefile == None):
		plt.show()
	else:
		savefile_name = str(object_ID)+'_wave_sed.png'
		plt.savefig(savefile_name, dpi = 300)

	if (savesed == True):
		np.savetxt(str(object_ID)+'_wave_sed.dat', np.c_[output_wavelength/1e4, output_flux])
		np.savetxt(str(object_ID)+'_wave_sed_photometry.dat', np.c_[output_phot_wavelength/1e4, output_phot_template])
	
	if (savechisq == True):
		np.savetxt(str(object_ID)+'_z_chisq.dat', np.c_[output_zgrid, output_chisq])
		

z_phot vs z_spec Comparison

If you happen to have spectroscopic redshifts for your sources, you could use the built-in function self.zphot_zspec():

fig = self.zphot_zspec()   # Produces a standard comparison plot

Interactive Dashboard

Similarly, there's a Plotly/Dash browser app that you can run which will allow you to explore the source fits interactively:

from eazy import visualization
app = visualization.EazyExplorer(photoz=self, zout=zout)
app.make_dash_app(server_mode='external', port=8050)
# Opens a Plotly/Dash browser app for interactive exploration

16 Common Pitfalls & Tips

for i, filt in enumerate(self.RES.filters):
    if 'f160w' in filt.name.lower():
        print(i+1, filt.name)   # +1 because 1-indexed

import pickle
with open("tempfilt_example.pkl","wb") as f:
      pickle.dump(self.tempfilt,f)

with open("tempfilt_example.pkl","rb") as inp:
      test_tempfilt = pickle.load(inp)

self = eazy.photoz.PhotoZ(params=params,
                   param_file=None,
                   translate_file='zphot.translate',
                   tempfilt = test_tempfilt, # <-- right here! 
                   zeropoint_file=None,
                   load_prior=False,
                   load_products=False,
                   n_proc=0)

unfitted = (self.chi2_fit == 0).all(axis=1)
print(f"{unfitted.sum()} objects not fit")

chi2_best = self.chi2_fit[np.arange(self.NOBJ), ez.izbest]
chinu     = chi2_best / np.maximum(self.nusefilt - 1, 1)
bad_fits  = chinu > 5   # Flag very poor fits

17 Key Internal Attributes Reference