diff --git a/docs/Project.toml b/docs/Project.toml
index 33d3b88..fb61401 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,15 +1,19 @@
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 InitialMassFunctions = "e9251ff4-c148-4db3-bf46-89407750fae0"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
 [compat]
 DelimitedFiles = "1"
+Distributions = "0.25"
 Documenter = "1"
 InitialMassFunctions = "0.1"
 Plots = "1"
-PyPlot = "2"
 Printf = "<0.0.1, 1"
+PyPlot = "2"
+StatsBase = "0.34"
diff --git a/docs/src/fitting/fitting_intro.md b/docs/src/fitting/fitting_intro.md
index 7b38635..531e99b 100644
--- a/docs/src/fitting/fitting_intro.md
+++ b/docs/src/fitting/fitting_intro.md
@@ -40,14 +40,19 @@ Below we show a comparison of a smooth Hess diagram template constructed with [`
 
 ![Comparison of smooth Hess diagram template from `partial_cmd_smooth` and a Monte Carlo model made with `generate_stars_mass`.](figures/model_cmd.png)
 
-
 ```@example
 mv("../../../examples/templates/template_compare.svg", "figures/template_compare.svg") # hide
 mv("../../../examples/templates/sigma_distribution.svg", "figures/sigma_distribution.svg") # hide
 nothing # hide
 ```
 
-![Caption.](figures/template_compare.svg)
+A worked example comparing a sampled stellar population with a smooth Hess diagram template is available in `examples/templates/smooth_template.jl`. The output figure is shown below. A distance modulus of 25 mag is used for this example, with photometric error and completeness functions roughly based on those we observe in the JWST/NIRCAM data of WLM (see [Weisz et al. 2024](https://ui.adsabs.harvard.edu/abs/2024ApJS..271...47W)). 
+
+![Comparison of CMD-sampled population with smooth Hess diagram template.](figures/template_compare.svg)
+
+At left is a population of stars sampled from an SSP with the methods described in the section of the documentation on [simulating CMDs](@ref simulate). The points from the isochrone are colored orange. The next figure shows the binned Hess diagram derived from these data. The next figure shows our smooth Hess diagram template calculated for this SSP. The final figure at right shows the residual between the data and model in units of standard deviations. These are sometimes called Pearson residuals. Below we show the distribution of these residuals, which should be Gaussian with mean 0 and standard deviation 1 if the model were perfect. We show a Gaussian PDF with standard deviation 1 and mean equal to the observed mean of the residuals for comparison.
+
+![Distribution of data - model residuals, in units of standard deviations.](figures/sigma_distribution.svg)
 
 
 ```@docs
diff --git a/examples/templates/smooth_template.jl b/examples/templates/smooth_template.jl
index d3806ab..55170b1 100644
--- a/examples/templates/smooth_template.jl
+++ b/examples/templates/smooth_template.jl
@@ -124,7 +124,7 @@ plt.savefig(joinpath(@__DIR__,"template_compare.svg"), bbox_inches="tight")
 #################################
 # Distribution of σ discrepancies
 import Distributions: pdf, Normal, Poisson
-import StatsBase: mean
+import StatsBase: mean, std
 fig, ax1 = plt.subplots()
 hist1 = ax1.hist(filter(isfinite, signif), range=(-4,4), bins=25, density=true)
 ax1.set_xlim( extrema(hist1[2]) )
@@ -132,9 +132,14 @@ ax1.set_xlabel("Residual / Standard Deviation")
 ax1.set_ylabel("PDF")
 let xplot = first(hist1[2]):0.01:last(hist1[2])
     # mean_σ = mean(filter(Base.Fix1(<,-5), filter(isfinite,signif)))
-    mean_σ = mean(filter(isfinite,signif))
-    ax1.plot( xplot, pdf.(Normal(mean_σ,1.0),xplot))
+    tmpgood = filter(x->first(hist1[2]) <= x <= last(hist1[2]), filter(isfinite,signif))
+    mean_σ = mean(tmpgood)
+    std_σ = std(tmpgood)
+    ax1.plot( xplot, pdf.(Normal(0.0,1.0),xplot), label="Ideal")
+    ax1.plot( xplot, pdf.(Normal(mean_σ,std_σ),xplot), label="Observed", c="magenta")
     ax1.axvline(mean_σ, c="k", ls="--")
+    ax1.text(0.05,0.95,(@sprintf("Mean: %.2f\n\$\\sigma\$: %.2f", mean_σ, std_σ)),transform=ax1.transAxes, va="top", ha="left")
 end
+ax1.legend(loc="upper right")
 plt.savefig(joinpath(@__DIR__,"sigma_distribution.svg"), bbox_inches="tight")