Post-processing and plotting results#

A RetrievalResult already exposes the retrieved field, spectrum and their intensities and phases. croak.processing goes further: it turns a result into a display-ready ProcessedResult with oversampled profiles, fitted dispersion coefficients, marginals, residuals and a spectral-edge quality flag, everything the figures need. Read these alongside the trace error, which measures agreement with the model rather than with the pulse.

pr = croak.process_result(result, measured=td.trace, Iomega_meas=td.Iomega)
pr.fwhm_retr, pr.fwhm_tl       # retrieved and transform-limited FWHM (s)
pr.gdd_fs2, pr.tod_fs3, pr.fod_fs4   # fitted dispersion (fs², fs³, fs⁴)

The spectral counterparts of fwhm_retr are two small helpers rather than fields, since they apply to any intensity-vs-wavelength curve — the retrieved spectrum, the measured one, or a known truth’s:

from croak.processing import peak_wavelength, spectral_fwhm

peak_wavelength(pr.wavelength, pr.Ilam)   # m
spectral_fwhm(pr.wavelength, pr.Ilam)     # m

Pass the wavelength density (Ilam, what the spectrum panel draws) rather than the frequency density Iw: the peak of one is not the peak of the other. Both helpers accept an ascending or descending axis, drop the unphysical bins where the grid’s 2πc/(ω+ω₀) axis crosses zero, and return nan when the quantity is undefined (an empty axis, or a spectrum whose half maximum the grid never reaches). spectral_fwhm takes an optional level= to measure a width at another fraction of the peak.

Absolute power and peak power#

A retrieval fixes the pulse shape but not its absolute scale, so by default the temporal intensity is normalised to a peak of 1 (“a.u.”). If you measured the pulse energy E (e.g. with a power meter), pass it as energy= (in joules) and process_result turns the normalised intensity into absolute instantaneous power. With the intensity normalised to unit peak, instantaneous power is P(t) = E · I(t) / ∫I(t)dt, so the peak power is E / ∫I dt — the energy divided by the effective duration:

pr = croak.process_result(result, measured=td.trace, energy=100e-6)  # 100 µJ
pr.peak_power       # retrieved-pulse peak power (W)
pr.peak_power_tl    # transform-limited peak power (W) at the same energy

For a fixed spectral amplitude and energy, a flat spectral phase gives the maximum peak power, so the transform-limited pulse sets the ceiling. Supply the energy of the same beam the retrieval describes (an on-axis field describes the centre of the beam, not its total power). Energy is conserved under the dispersion-stage (phase-only) compensation, so compressing the pulse raises its peak power towards peak_power_tl at fixed E. With energy unset (the default), peak_power/peak_power_tl are None and every plot stays in normalised units.

plot_temporal() (and plot_retrieval() via its own energy= argument) automatically switch the temporal axis to absolute power with an SI prefix — chosen by si_power() — when a peak power is available.

A known ground truth scales the same way through peak_power(), which applies the same E / ∫(I/I_peak) dt definition over the truth’s own time axis:

truth.peak_power(energy)      # W, or None when the energy is unknown

What process_result computes#

process_result() fills a ProcessedResult with:

  • Temporal — peak-centred, oversampled intensity It_retr and phase phi_t (linear ramp removed), the transform-limited reference It_tl, and the FWHMs fwhm_retr / fwhm_tl. When a pulse energy is supplied, also the absolute peak_power / peak_power_tl (W) — see below.

  • Spectral (wavelength domain) — frequency density Iw and wavelength density Ilam (including the \(|d\omega/d\lambda|\) Jacobian), phase phi_w with a quartic fit phi_w_fit, and the dispersion coefficients gdd_fs2 / tod_fs3 / fod_fs4 from that fit. If you pass Iomega_meas, the measured spectrum Iw_meas is included for comparison.

  • Trace — the retrieved trace, the measured trace and their marginals and weighted residual, plus the convergence history errors.

  • Quality — edge_energy, the fraction of the retrieved pulse energy sitting in the unmeasured grid-edge bins (see below).

Helpers are also exposed individually: marginals(), residuals(), oversample(), spectrogram() and shiftnorm().

Is the spectrum real#

A low trace error does not by itself mean the retrieved spectrum is trustworthy. The outermost TAPER_COLLAR_BINS (10) frequency bins at each end of the grid carry almost no measurement, because the edge taper takes them to zero. The field there is nearly free, and because intensity cannot go below zero it drifts upward: spectral energy appears at the band edges that was never in the trace.

edge_energy_fraction() measures it, and ProcessedResult.edge_energy carries it. plot_retrieval() prints it in the spectrum panel, in red above croak.plotting.EDGE_ENERGY_WARN (1 %):

pr = croak.process_result(result, measured=td.trace)
if pr.edge_energy > 0.01:
    print(f"{pr.edge_energy:.1%} of the pulse energy is off the measured band")

Rules of thumb from our deep-UV TG-FROG measurements: under 0.1 % is a sound retrieval; a few percent means part of the reported spectrum, and therefore of the reported duration, is an artefact. On one dataset, 15.6 % edge energy inflated the reported FWHM by 14 % against the same retrieval with the edge bins removed. A low value rules out this failure only; it does not certify the rest of the field.

The causes, in order of importance:

  1. R_omega=True used without reg_spectrum or phase_only: by far the biggest contributor, and it lowers the reported error while doing it. See per-frequency scaling.

  2. A lam_min/lam_max band tight enough that the taper cuts live signal; load_and_clean warns about this at load time.

  3. Noise fitted into signal-free bins on a noisy trace when no measured spectrum is supplied as a target. The COPRA warm-up of warm-lbfgs has no spectral penalty even when one is; see the thickness-uncertainty example.

Note the metric counts bins, so it is only meaningful on a grid with room for two non-overlapping collars (≥ 21 points); it returns 0.0 below that.

Plotting#

croak.plotting is built on Matplotlib without global pyplot state, so the figures embed cleanly (the GUI uses the same functions). The headline composite is the twelve-panel plot_retrieval():

fig = croak.plot_retrieval(result, measured=td.trace, Iomega_meas=td.Iomega,
                          lam_min=td.lam_min, lam_max=td.lam_max)
fig.savefig("retrieval.png", dpi=110)

Its annotations are the numbers you judge a retrieval by, so they are printed finely enough to compare runs:

  • the retrieved-trace panel carries the FROG error as a percentage to two decimals (0.71%), above the trace dimensions;

  • the spectrum panel reports the fitted GDD (fs²) and TOD (fs³) to one decimal, matching the dashed polynomial phase fit drawn over the retrieved phase;

  • the temporal panel’s legend gives the intensity FWHM of every curve — retrieved R, transform-limited TL, and (when a ground truth is passed) the truth overlay — so the three durations read side by side.

Pass truth= (a TruthPulse) to draw a known pulse beneath the retrieved one. The synthetic and simulated GUI workflows do this automatically, but it is an ordinary library feature — build one with TruthPulse.from_spectrum whenever you know the answer:

truth = croak.TruthPulse.from_spectrum(grid, ew, omega0)
fig = croak.plot_retrieval(result, measured=trace, truth=truth)

SHG: settle the direction of time first#

An SHG-FROG trace is exactly invariant under \(E(t) \to E^*(-t)\), so a pulse and its mirror image fit the measurement equally well and the solver returns whichever branch it happened to land in. Overlay a truth without accounting for that and a perfect retrieval can look like a failure — the two curves mirrored about \(t=0\).

When the true pulse is known, resolve_time_direction() picks the matching branch. Apply it to the RetrievalResult before anything derives from it, so the plot, any ProcessedResult and the reported dispersion all describe the same pulse:

result, flipped = croak.resolve_time_direction(result, truth)
pr = croak.process_result(result, measured=trace)

It compares temporal intensity profiles by their best normalised cross-correlation over all lags, which keeps the decision independent of the absolute phase and the time origin (both ambiguous too), and reliable on an imperfect retrieval: only the relative ranking of the two branches has to be right. A time-symmetric intensity, such as a chirped Gaussian’s, scores alike on both branches and is left unchanged. It is a no-op for PG and SD, which fix the direction of time.

Reversal is applied by conjugating the spectrum, which is exact: the stored trace, mu and error remain valid. Note that the spectral phase changes sign, so the retrieved GDD and TOD only become comparable with the truth once the direction is settled. The wizard does this automatically for synthetic and simulated sessions, where a truth is available.

Framing the axes#

Two independent controls, because the panels use two different abscissae:

  • lam_min / lam_max (m) bound the wavelength axes — the spectrum and spectrogram panels. On the measured-data path take them from the cleaned TraceData, as above.

  • flim — (f_min, f_max) in Hz — bounds the frequency axis of the trace and residual panels.

They are separate because the trace sits at the signal frequency, which for SHG is twice the pulse’s. flim matters when the retrieval grid is much wider than the signal, which is typical when retrieving directly on a synthesis grid rather than one cropped by load_and_clean(); without it the trace is a thin stripe across a mostly empty axis:

f0 = omega0 / (2 * np.pi)
fig = croak.plot_retrieval(result, measured=trace,
                           flim=(2 * f0 - 0.3e15, 2 * f0 + 0.3e15))  # SHG

Panels and pages#

The twelve panels are independent functions over one data bundle, so any subset can be drawn. retrieval_plot_data() computes everything the panels share once — the processed result, axis bounds, colour scales, and the spectrogram on first use; plot_retrieval_page() lays a PageSpec of panel keys out as a figure; and each draw_<key>(ax, data) draws one panel into an axis of your own:

data = croak.plotting.retrieval_plot_data(result, measured=td.trace,
                                          Iomega_meas=td.Iomega,
                                          lam_min=td.lam_min, lam_max=td.lam_max)
fig = croak.plotting.plot_retrieval_page(data, croak.plotting.RETRIEVAL_PAGES["pulse"])

ax = Figure().add_subplot()
croak.plotting.draw_spectral(ax, data)   # one panel, in your own figure

RETRIEVAL_PANELS maps the keys — measured_log, retrieved_log, measured_lin, retrieved_lin, temporal, spectral, convergence, residual, spectral_filter, freq_marginal, delay_marginal, spectrogram — to their titles and draw functions; RETRIEVAL_PAGES holds the three 2×2 pages the wizard shows as tabs (Traces, Pulse, Diagnostics), and RETRIEVAL_OVERVIEW the 3×4 layout plot_retrieval draws. The spectrum panel is drawn against frequency by default; pass spectral_axis="wavelength" to plot_retrieval or retrieval_plot_data (or axis="wavelength" to plot_spectral) for the spectrometer’s axis — both are unit-peak densities in their own variable, so the wavelength view carries the \(1/\lambda^2\) Jacobian and the frequency view does not. Colourbars shared by a pair of trace panels are drawn whenever both panels are on the page. Every draw_* returns the artists it created (for example TemporalArtists), which is what a live view updates in place instead of redrawing.

Other composites and single-axis plotters:

Function

Shows

plot_retrieval()

12-panel overview (trace, residual, time, spectrum, marginals, convergence)

plot_retrieval_page()

one page of retrieval panels (a PageSpec)

retrieval_plot_data()

the shared inputs of the retrieval panels

draw_measured_log, …, draw_spectrogram

one retrieval panel each (see RETRIEVAL_PANELS)

plot_frog_filter()

6-panel preprocessing before/after

plot_simulated_trace()

retrieved vs measured trace

plot_residual

a single residual image

plot_temporal, plot_spectral

intensity + phase in one axis (with a truth overlay); the spectrum against axis="frequency" (default) or "wavelength"

plot_convergence

trace error vs iteration

plot_marginal, plot_spectrogram

a marginal / Gabor spectrogram

The single-axis plotters take an existing Matplotlib axis, so you can compose your own figures.

Live views (GUI support)#

For custom interactive applications, see the live-view contracts in the developer’s guide.

Windowing the retrieved field#

To trim ringing or out-of-band energy from the result before plotting (without re-running the retrieval), apply croak.processing.post_filter() — a reversible temporal/spectral Planck taper that returns a new result. Pass the measurement windows td.lamm_lims / td.tau_lims (λm/τm) as lam_lims / tau_lims. The stored trace and error still describe the unfiltered retrieval. See Grids and filtering.

Residual dispersion#

process_result reports the fitted GDD/TOD/FOD. To explore compressing the pulse — applying extra dispersion and watching the peak power — see the dispersion-tuning tutorial and croak.dispersion.

Saving#

Persist a result (and its processed form) to HDF5 with croak.save.save_result(); see Saving and loading.