diff --git a/src/toolbox_scs/detectors/hrixs.py b/src/toolbox_scs/detectors/hrixs.py
index ca800e8595514d13243b7ff3452692ce86531031..96a0a8ed13c24ab59cf86f8ed94d6ed88072770a 100644
--- a/src/toolbox_scs/detectors/hrixs.py
+++ b/src/toolbox_scs/detectors/hrixs.py
@@ -493,28 +493,73 @@ class hRIXS:
return data
- def align_readouts(self,data,method,start,stop,fit_tol=None):
+ def align_readouts(self,data_list,start,stop,method='max_value',ImageGrouping=1,fit_tol=None):
+ import xarray as xa
import scipy as sp
from scipy import optimize
#********************************************
- # aligns spectra in a given data xarray
+ # aligns spectra in a given list of data from
+ #
+ # data = h.from_run(runind,2776)
+ # OR
+ # data = h.centroid(data,return_hits=True)
+ #
+ # One can choose to handle
+ # Grouping = 'run'
+ # Grouping = integer (1 means each image separately)
# METHOD
# -max_value
# -autocorrelation
# -gauss_fit
# start and stop are values of data.energy
# that define the range of these operations
- # RETURNS LINE CENTER POSITIONS
- # (and in future perhaps shifted spectra)
+ # RETURNS:
+ # energy grid and shifted spectra
#********************************************
+ # Accept a list of data xarrays. I needed
+ # to write it this way as xarray was hard
+ #********************************************
+ energy = data_list[0].energy.to_numpy()
+ if type(ImageGrouping)==type(1):
+ #********************************************
+ # Group data to bunches of a given number of
+ # Images
+ #********************************************
+ data = xa.concat(data_list,dim='trainId')
+ spectra = []
+ xhits = []
+ Ngroups = data.spectrum.shape[0]//ImageGrouping
+ # First the groups with full length
+ for groupind in range(Ngroups):
+ spectra.append(np.sum(data.spectrum[groupind*ImageGrouping:(groupind+1)*ImageGrouping].to_numpy(),axis=0))
+ xhits.append(np.hstack(data.xhits[groupind*ImageGrouping:(groupind+1)*ImageGrouping].to_numpy()))
+ # if the leftover tail is more than a half of full length, add it too
+ if data.spectrum.shape[0]%ImageGrouping > ImageGrouping/2:
+ spectra.append(np.sum(data.spectrum[Ngroups*ImageGrouping::].to_numpy(),axis=0))
+ xhits.append(np.hstack(xhits[Ngroups*ImageGrouping::].to_numpy()))
+ elif ImageGrouping.lower() == 'run' or ImageGrouping.lower() == 'runs':
+ #********************************************
+ # Group data by runs
+ #********************************************
+ spectra = []
+ xhits = []
+ for rundata in data_list:
+ runspectrum = np.zeros(energy.shape)
+ runxhits = []
+ for spec,hits in zip(rundata.spectrum,rundata.xhits):
+ runspectrum += spec.to_numpy()
+ runxhits.append(hits.to_numpy())
+ spectra.append(runspectrum)
+ xhits.append(np.hstack(runxhits))
+ else:
+ raise Exception('align_readouts() needs a reasonable grouping argument')
#********************************************
#********************************************
- # PART 1: find shifts
+ # PART 1: find peak positions
#********************************************
#********************************************
- searchinds = (data.energy >=start)*(data.energy <=stop)
+ searchinds = (energy >=start)*(energy <=stop)
peak_posis = []
- chan_shift = []
#********************************************
# Simple maximum alignment
#********************************************
@@ -522,16 +567,13 @@ class hRIXS:
#********************************************
# Find the max for each of the spectra
#********************************************
- for ind,spec in enumerate(data.spectrum):
- x = data.energy.to_numpy()[searchinds]
- y = spec.to_numpy()[searchinds]
+ for ind,spec in enumerate(spectra):
+ x = energy[searchinds]
+ y = spec[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
@@ -541,22 +583,20 @@ class hRIXS:
#********************************************
# Find the max for each of the spectra
#********************************************
- for ind,spec in enumerate(data.spectrum):
+ for ind,spec in enumerate(spectra):
if ind == 0:
- x0 = data.energy.to_numpy()[searchinds]
- y0 = spec.to_numpy()[searchinds]
- maxipos0 = np.argmax(spec.to_numpy()[searchinds])
+ x0 = energy[searchinds]
+ y0 = spec[searchinds]
+ maxipos0 = np.argmax(spec[searchinds])
peak_posis.append(x0[maxipos0])
- chan_shift.append(0)
else:
- x = data.energy.to_numpy()[searchinds]
- y = spec.to_numpy()[searchinds]
+ x = energy[searchinds]
+ y = spec[searchinds]
corr_len = np.sum(searchinds)
corr = sp.signal.correlate(y,y0, mode='full')
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.')
@@ -565,18 +605,19 @@ class hRIXS:
#********************************************
def Gauss(grid,x0,sigma):
# Returns a normalized bell curve
- # with center at x0 and sigma
+ # with center at x0 and std sigma
# on grid
return 1.0/(sigma*np.sqrt(2.0*np.pi))*np.exp(-0.5*(grid-x0)**2/sigma**2)
def Cost(p,grid,spec):
+ # Returns squared error of spec and a bell curve presented on grid
return np.sum(np.square(p[0]*Gauss(grid,p[1],p[2])-spec))
#********************************************
# Find the max for each of the spectra
#********************************************
- for spec in data.spectrum:
- x = data.energy.to_numpy()[searchinds]
- y = spec.to_numpy()[searchinds]
+ for spec in spectra:
+ x = energy[searchinds]
+ y = spec[searchinds]
#********************************************
# Initial Guess and bounds
#********************************************
@@ -596,6 +637,9 @@ class hRIXS:
if p.success:
peak_posis.append(p.x[1])
else:
+ #********************************************
+ # If fitting fails plot why and quit
+ #********************************************
plt.figure()
plt.plot(x,y,'.')
plt.plot(x,p.x[0]*Gauss(x,p.x[1],p.x[2]))
@@ -620,27 +664,27 @@ class hRIXS:
# we can subtract and interpolation will
# not cause a problem.
#********************************************
- energy_aligned=data.energy.to_numpy()-peak_posis[0]
+ energy_aligned=energy-peak_posis[0]
ret = []
- for position,spec in zip(peak_posis,data.spectrum):
+ for position,spec in zip(peak_posis,spectra):
grid = energy_aligned
- x = data.energy.to_numpy() - position
- y = spec.to_numpy()
+ x = energy - position
+ y = spec
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,
+ for position,hits in zip(peak_posis,xhits):
+ spectrum = np.histogram(hits-position,bins=self.BINS,
range=(0-e0, self.Y_RANGE.stop - self.Y_RANGE.start-e0))[0]
- ret.append(spec)
+ ret.append(spectrum)
spectrum_aligned = np.array(ret)
else:
raise Exception('align_readouts() did not recognize the method.')
- return energy_aligned,spectrum_aligned
+ return energy_aligned,spectrum_aligned
def integrate(self, data):