diff --git a/XAS.py b/XAS.py index 2ff3c914555f848828a7075431e5ee6c8857af03..9170a263bf54a6e9ec60fde9bd3a3c709d43f41d 100644 --- a/XAS.py +++ b/XAS.py @@ -47,7 +47,6 @@ def absorption(T, Io, fluorescence=False): good = np.logical_and(np.isfinite(T), np.isfinite(Io)) T = T[good] Io = Io[good] - # return type of the structured array fdtype = [('muA', 'f8'), ('sigmaA', 'f8'), ('weights', 'f8'), ('muT', 'f8'), ('sigmaT', 'f8'), ('muIo', 'f8'), ('sigmaIo', 'f8'), @@ -64,21 +63,26 @@ def absorption(T, Io, fluorescence=False): weights = np.sum(Io) p = np.corrcoef(T, Io)[0,1] - - if fluorescence: - muA = muT/muIo - sigmaA = np.abs(muA)*(np.sqrt((sigmaT/muT)**2 + - (sigmaIo/muIo)**2 - 2*p*sigmaIo*sigmaT/(muIo*muT))) - else: - muA = -np.log(muT/muIo) - # from error propagation for correlated data - sigmaA = (np.sqrt((sigmaT/muT)**2 + (sigmaIo/muIo)**2 - - 2*p*sigmaIo*sigmaT/(muIo*muT))) - # direct calculation - #mask = (Io != 0) - #sigmaA = np.nanstd(-np.log(T[mask]/Io[mask])) + # weighted average of T/Io with Io as weights + muA = muT/muIo + + #Derivation of standard deviation + #1. using biased weighted sample variance: + #sigmaA = np.sqrt(np.average((T/Io - muA)**2, weights=Io)) + #2. using unbiased weighted sample variance (reliablility weights): + V2 = np.sum(Io**2) + sigmaA = np.sqrt(np.sum(Io*(T/Io - muA)**2) / (weights - V2/weights)) + + #3. using error propagation for correlated data: + #sigmaA = np.abs(muA)*(np.sqrt((sigmaT/muT)**2 + + # (sigmaIo/muIo)**2 - 2*p*sigmaIo*sigmaT/(muIo*muT))) + + if not fluorescence: + sigmaA = sigmaA/np.abs(muA) + muA = -np.log(muA) + return np.array([(muA, sigmaA, weights, muT, sigmaT, muIo, sigmaIo, p, counts)], dtype=fdtype) @@ -97,7 +101,6 @@ def binning(x, data, func, bins=100, bin_length=None): bins: the bins edges res: a structured array of binned data """ - if bin_length is not None: bin_start = np.amin(x) bin_end = np.amax(x) @@ -108,7 +111,7 @@ def binning(x, data, func, bins=100, bin_length=None): bins = np.linspace(bin_start, bin_end, bins) bin_centers = (bins[1:]+bins[:-1])/2 nb_bins = len(bin_centers) - + bin_idx = np.digitize(x, bins) dummy = func([]) res = np.empty((nb_bins), dtype=dummy.dtype) @@ -187,7 +190,7 @@ def xas(nrun, bins=None, Iokey='SCS_SA3', Itkey='MCP3apd', nrjkey='nrj', Iooffse dummy, nosample = binning(rundata['nrj'], rundata, whichIo, bins) muA = nosample['muA'] sterrA = nosample['sigmaA']/np.sqrt(nosample['counts']) - + bins_c = 0.5*(bins[1:] + bins[:-1]) if plot: