diff --git a/pes_to_spec/__init__.py b/pes_to_spec/__init__.py
index 02af2a4e342acba42b23495baf4b64f079fd1a1c..d5e0f8e872f3dd0799f2d1bed563b7fe249070b5 100644
--- a/pes_to_spec/__init__.py
+++ b/pes_to_spec/__init__.py
@@ -2,4 +2,4 @@
 Estimate high-resolution photon spectrometer data from low-resolution non-invasive measurements.
 """
 
-VERSION = "0.3.7"
+VERSION = "0.3.8"
diff --git a/pes_to_spec/model.py b/pes_to_spec/model.py
index 44d7a1487602e9189125c4841ddf74f4cee70cb1..422a6021940de4970486f699e976cfe0660d3c73 100644
--- a/pes_to_spec/model.py
+++ b/pes_to_spec/model.py
@@ -79,14 +79,14 @@ def deconv(y: np.ndarray, yhat: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.
 
 def fit_gaussian(x: np.ndarray, y: np.ndarray) -> lmfit.ModelResult:
     """Fit Gaussian and return the fit result."""
-    def gaussian(x, amplitude, centre, sigma):
-        return amplitude * np.exp(-0.5 * (x - centre)**2 / (sigma**2))
+    def gaussian(x, amplitude, centre, log_sigma):
+        return amplitude * np.exp(-0.5 * ((x - centre)**2)*np.exp(-2*log_sigma))
     gmodel = lmfit.Model(gaussian)
-    result = gmodel.fit(y, x=x, centre=0.0, amplitude=np.amax(y), sigma=1.0)
+    result = gmodel.fit(y, x=x, centre=0.0, amplitude=np.amax(y), log_sigma=0.0)
     return result
 
 def get_resolution(y: np.ndarray, y_hat: np.ndarray, e: np.ndarray,
-                   e_center: Optional[float]=None, e_width: Optional[float]=None) -> Tuple[np.ndarray, np.ndarray, np.ndarray, lmfit.ModelResult]:
+                   e_center: Optional[float]=None, e_width: Optional[float]=None) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, lmfit.ModelResult]:
     """
     Given the true y and the predicted y, together with the energy axis e,
     estimate the impulse response of the system and return the Gaussian fit to it.
@@ -99,7 +99,7 @@ def get_resolution(y: np.ndarray, y_hat: np.ndarray, e: np.ndarray,
       e_center: If given the energy value, for which to probe the resolution.
       e_width: The width of the energy neighbourhood to probe if e_center is given.
 
-    Returns: The centered energy axis, the impulse response, the transfer function and the fit result.
+    Returns: The centered energy axis, the impulse response, the transfer function, the grating spectral density and the fit result.
     """
     e_range = e[-1] - e[0]
     e_axis = np.linspace(-0.5*e_range, 0.5*e_range, len(e))
@@ -112,7 +112,7 @@ def get_resolution(y: np.ndarray, y_hat: np.ndarray, e: np.ndarray,
         y_sel = y_sel*f
         y_hat_sel = y_hat_sel*f
     h, H, S = deconv(y_sel, y_hat_sel)
-    return e_axis, h, H, fit_gaussian(e_axis, np.absolute(h))
+    return e_axis, h, H, S, fit_gaussian(e_axis, np.absolute(h))
 
 class PromptNotFoundError(Exception):
     """
@@ -796,7 +796,7 @@ class Model(TransformerMixin, BaseEstimator):
                               "n_obs",
                               "pca_threshold",
                               "high_res_fwhm",
-                              "resolution_per_energy", "resolution_per_energy_unc"
+                              "resolution_per_energy", "resolution_per_energy_unc",
                               "resolution_energy_bins"
                               ]
 
@@ -1014,10 +1014,10 @@ class Model(TransformerMixin, BaseEstimator):
         e = high_res_photon_energy[0,:] if len(high_res_photon_energy.shape) == 2 else high_res_photon_energy
 
         # get average resolution
-        e_axis, h, H, result = get_resolution(y,
-                                              y_hat,
-                                              e
-                                              )
+        e_axis, h, H, S, result = get_resolution(y,
+                                                 y_hat,
+                                                 e
+                                                 )
 
         de = e[1] - e[0]
         E = np.fft.fftfreq(len(e), de)
@@ -1063,23 +1063,23 @@ class Model(TransformerMixin, BaseEstimator):
         #if self.resolution < 0:
         #    print("Warning: Resolution calculation failed. The model can still be used, but this is probably a red flag!")
         #
-        self.resolution = result.best_values["sigma"]*2.355
+        self.resolution = np.exp(result.best_values["log_sigma"])*2.355
         print("Resolution:", self.resolution)
 
         # estimate resolution using power spectra
         e_min = e.min()
         e_max = e.max()
         e_probe = np.linspace(e_min, e_max, 5)
-        e_width = (e_max - e_min)/(len(e_probe[key])-1)
+        e_width = (e_max - e_min)/(len(e_probe)-1)
         width = list()
         width_unc = list()
         for e_ in e_probe:
-            e_axis_centered, h, H, result = get_resolution(y,
-                                                           y_hat,
-                                                           e_axis,
-                                                           e_center=e_,
-                                                           e_width=e_width)
-            width += [result.best_values["sigma"]*2.355]
+            _, _, _, _, result = get_resolution(y,
+                                                y_hat,
+                                                e,
+                                                e_center=e_,
+                                                e_width=e_width)
+            width += [np.exp(result.best_values["log_sigma"])*2.355]
             width_unc += [np.sqrt(result.covar[2,2])*2.355]
         width = np.array(width)
         width_unc = np.array(width_unc)