feat: unfolded spectrum (wip with Fergus's edits)

src/plots-recipes.jl

@@ -224,11 +224,20 @@ end
 
     result = r.args[1] isa FittingResult ? r.args[1][1] : r.args[1]
     data = get_dataset(result)
-    x = plotting_domain(data)
+    model = get_model(result)
+
+    x_plot = plotting_domain(data)
     y = invoke_result(result, result.u)
+    x_domain = make_output_domain(ContiguouslyBinned(), data)
+    # TODO: pass the result.u to invokemodel (currently not possible without knowing frozen paramters)
+    update_model!(model, result)
+    y_unfolded = invokemodel(x_domain, model)[data.data.data_mask]
+
+    unfolded_spectrum = @. x_plot^0 * result.objective * (y_unfolded / y)
+    unfolded_std = @. sqrt(result.variance) * x_plot^0 * result.objective * (y_unfolded / y)
 
-    y_unfolded = @. x^2 * result.objective / y # <-- allow power of x to be a parameter
-    # TODO: multiply this by the unfolded model
+    unfolded_spectrum = @. y_unfolded 
+    unfolded_std = @. unfolded_spectrum / 100
 
     ylabel --> "Unfolded spectrum [E F_E]" # <-- check units
     xlabel --> "Energy (keV)"
@@ -238,23 +247,15 @@ end
         label = make_label(data)
     end
 
-    @series begin
-        linestyle --> :dash
-        seriestype --> :hline
-        label --> false
-        color --> modelcolor
-        [1.0]
-    end
-
     @series begin
         markerstrokecolor --> datacolor
         label --> label
         seriestype --> :scatter
         markershape --> :none
         markersize --> 0.5
-        yerror --> @. sqrt(result.variance) * x^2 * result.objective / y
+        yerror --> unfolded_std
         xerror --> SpectralFitting.bin_widths(data) ./ 2
-        x, y_unfolded
+        x_plot, unfolded_spectrum
     end
 end