Alignment routine runs, but needs a careful validation & debugging.

src/toolbox_scs/detectors/hrixs.py

@@ -493,7 +493,7 @@ class hRIXS:
         return data
 
     
-    def align_readouts(self,data,method,start,stop):
+    def align_readouts(self,data,method,start,stop,fit_tol=None):
         import scipy as sp
         from scipy import optimize 
         #********************************************
@@ -507,8 +507,14 @@ class hRIXS:
         # RETURNS LINE CENTER POSITIONS 
         # (and in future perhaps shifted spectra)
         #********************************************
+        #********************************************
+        #********************************************
+        # PART 1: find shifts
+        #********************************************
+        #********************************************
         searchinds = (data.energy >=start)*(data.energy <=stop)
         peak_posis = []
+        chan_shift = []
         #********************************************
         # Simple maximum alignment
         #********************************************
@@ -516,11 +522,16 @@ class hRIXS:
             #********************************************
             # Find the max for each of the spectra
             #********************************************
-            for spec in data.spectrum:
+            for ind,spec in enumerate(data.spectrum):
                 x          = data.energy.to_numpy()[searchinds]
                 y          = spec.to_numpy()[searchinds]
                 maxipos    = np.argmax(y)
                 peak_posis.append(x[maxipos])
+                if ind==0:
+                    maxipos0 = maxipos
+                    chan_shift.append(0)
+                else:
+                    chan_shift.append(maxipos-maxipos0)
         #********************************************
         # Alignment based on autocorrelation
         # this is a relative alignment method
@@ -536,6 +547,7 @@ class hRIXS:
                     y0       = spec.to_numpy()[searchinds]
                     maxipos0 = np.argmax(spec.to_numpy()[searchinds])
                     peak_posis.append(x0[maxipos0])
+                    chan_shift.append(0)
                 else:        
                     x        = data.energy.to_numpy()[searchinds]
                     y        = spec.to_numpy()[searchinds]
@@ -544,7 +556,10 @@ class hRIXS:
                     maxpos   = np.argmax(corr)
                     shift    = maxpos-corr_len
                     peak_posis.append(x[maxipos0+shift])
+                    chan_shift.append(shift)
         elif method.lower() == 'gauss_fit':
+            if fit_tol == None:
+                raise Exception('Gauss fit requires a tolerance value.')
             #********************************************
             # Define needed functions
             #********************************************
@@ -576,7 +591,7 @@ class hRIXS:
                 #********************************************
                 # Fit by minimizing least squares error
                 #********************************************
-                p          = optimize.minimize(Cost,p0,args=(x,y),bounds=bnds,method='L-BFGS-B',tol=1e-6,
+                p          = optimize.minimize(Cost,p0,args=(x,y),bounds=bnds,method='L-BFGS-B',tol=fit_tol,
                                             options={'disp':0,'maxiter':1000000})
                 if p.success:
                     peak_posis.append(p.x[1])
@@ -587,7 +602,45 @@ class hRIXS:
                     raise Exception('align_readouts(): can not fit a gaussian to the data.')
         else: 
             raise Exception('align_readouts() did not recognize the method.')
-        return peak_posis    
+        #********************************************
+        #********************************************
+        # PART 2: shift aquisitions
+        # For max_value and autocorrelation,
+        # the shift is integer bins and no rebinning
+        # is needed
+        # FOR gaussian fit, rebinning the given
+        # hit lists
+        #********************************************
+        #********************************************
+        if method.lower() == 'max_value' or method.lower() == 'autocorrelation':
+            #********************************************
+            # Since peak_posis[0] matches with one grid
+            # point, the shift in energy scale will be
+            # an integer multiple of grid spacing. Thus
+            # we can subtract and interpolation will 
+            # not cause a problem.
+            #********************************************
+            energy_aligned=data.energy.to_numpy()-peak_posis[0]
+            ret = []
+            for position,spec in zip(peak_posis,data.spectrum):
+                grid = energy_aligned
+                x    = data.energy.to_numpy() - position
+                y    = spec.to_numpy()
+                ret.append(np.interp(grid,x,y))
+            spectrum_aligned = np.array(ret)
+        elif method.lower() == 'gauss_fit':
+            e0  = peak_posis[0]
+            energy_aligned = np.linspace(self.Y_RANGE.start, self.Y_RANGE.stop, self.BINS)-e0
+            ret = []
+            for position,hits in zip(peak_posis,data.xhits):
+                spec = np.histogram(hits.to_numpy()-position,bins=self.BINS,
+                                    range=(0-e0, self.Y_RANGE.stop - self.Y_RANGE.start-e0))[0]
+                ret.append(spec)
+            spectrum_aligned = np.array(ret)
+        else:
+            raise Exception('align_readouts() did not recognize the method.')
+        
+        return energy_aligned,spectrum_aligned   
     
     
     def integrate(self, data):