diff --git a/src/toolbox_scs/detectors/hrixs.py b/src/toolbox_scs/detectors/hrixs.py
index 89218a10d8b28352d1bbb00f1b5d358d1bd32c01..31f6984f3b653e470d53e4130c9c19a682d192c4 100644
--- a/src/toolbox_scs/detectors/hrixs.py
+++ b/src/toolbox_scs/detectors/hrixs.py
@@ -305,6 +305,10 @@ class hRIXS:
     
     # Early energy calibration from runs 31-36 (to estimate mono drift)
     PIX2EV_POLY = [6.31196512e-02, 6.05502748e+02]
+    # Width of fit window to calculate hv (estimated mono ± this)
+    HV_FIT_DELTA = 30
+    # Gaussian fit parameters for elastic line finding
+    HV_FIT_SIGMA = 2
 
     ENERGY_INTERCEPT = 0
     ENERGY_SLOPE = 1
@@ -413,7 +417,8 @@ class hRIXS:
         return self.CURVE_A, self.CURVE_B
 
     
-    def centroid(self, data, bins=None, return_hits=False):
+    def centroid(self, data, bins=None, return_hits=False, fit_elastic=False,
+                 hv=None, fit_delta=None):
         #*************************************************************
         # Carry out hit finding on data and bin them in grid
         # Allows for 
@@ -436,10 +441,30 @@ class hRIXS:
         hit_y = []
         hits = []
         ret = np.zeros((len(data["hRIXS_det"]), bins))
+        #*************************************************************
+        # If the 'actual' photon energy is desired, we'll need to know
+        # where to start looking for the elastic line
+        #*************************************************************
+        if fit_elastic:
+            elastic_pixel = np.zeros(len(data["hRIXS_det"]))
+            found_hv = np.zeros(len(data["hRIXS_det"]))
+            if fit_delta == None:
+                fit_delta = self.HV_FIT_DELTA
+            if hv == None:
+                # ?? We'll make a guess later
+                hv = -1 * np.ones(len(data["hRIXS_det"]))
+            else:
+                # hv is given as one value or a train-vector
+                if type(hv) == type([]):
+                    if len(hv) != len(data["hRIXS_det"]):
+                        raise ValueError("Size of hv vector doesn't match trainID")
+                else:
+                    hv = hv * np.ones(len(data["hRIXS_det"]))
+                
         #*************************************************************
         # Handle each Aq image separately
         #*************************************************************
-        for image, r in zip(data["hRIXS_det"], ret):
+        for i, (image, r) in enumerate(zip(data["hRIXS_det"], ret)):
             use_image = image.to_numpy()
             #*************************************************************
             # Treat background by optionally 
@@ -493,6 +518,25 @@ class hRIXS:
                 rc[:, 0], bins=bins,
                 range=(0, self.Y_RANGE.stop - self.Y_RANGE.start))
             r[:] = hy
+            #*************************************************************
+            # Return the found elastic peak energy, if requested
+            #*************************************************************
+            if fit_elastic:
+                if hv[i] == -1:
+                    # Assume elastic line to be the strongest peak
+                    gx = np.where(hy == max(hy))[0][0]
+                else:
+                    # Use given guess for where it is
+                    gx = np.where(hx + self.Y_RANGE.start > self.energy2pixel(hv[i]))[0][0]
+                popt, pcov = curve_fit(gauss1d, hx[gx-fit_delta:gx+fit_delta],
+                                       hy[gx-fit_delta:gx+fit_delta],
+                                      p0 = [max(hx[gx-fit_delta:gx+fit_delta]),
+                                           hx[gx], self.HV_FIT_SIGMA, 0],
+                                      bounds=([0, 0, 0, -np.Inf],[np.Inf, np.Inf, 100, np.Inf]))
+                
+                elastic_pixel[i] = popt[1]+ self.Y_RANGE.start
+                found_hv[i] = self.pixel2energy(popt[1]+ self.Y_RANGE.start)
+                                       
         #*************************************************************
         # Setup and assing a linear energy grid
         #*************************************************************
@@ -507,6 +551,12 @@ class hRIXS:
                         xhits=(("trainId"), hit_x),
                         yhits=(("trainId"), hit_y))
         #**********************************************
+        # If hv was fitted, return it
+        #**********************************************
+        if fit_elastic:
+            data = data.assign(fitted_hv = (("trainId"), found_hv),
+                              elastic_pixel = (("trainId"), elastic_pixel))
+        #**********************************************
         # Always assign the spectrum to data
         #**********************************************
         data = data.assign(spectrum=(("trainId", "energy"), ret))