Spectrum descriptors for analysis¶
def load_bt_para(bt): thr = [3,2.0] if bt == ‘LR81’: adu_int = [ 97.0,123.0,False] roi = [310,348, 0,709]; fit = [284,304,358,378] axis_para = [50.72,498.21] else: print(‘Beamtime {} not available!’.format(bt)) return roi,fit,thr,axis_para,adu_int
Difference spectra¶
foo
First moment¶
bar
Note
Check observables (first moments) to established values to ensure that manual detector geometry works well.
FWHM¶
foobar
IAD¶
barfoo
Error bar estimation¶
Boot-strap¶
Fit function¶
########################
# smooth spectra and
# difference spectra
########################
def smooth_spec(xaxis,raw,bstep,xstep):
xmin = Emin - 2.5*bstep; xmax = Emax + 2.5*bstep
bins = np.zeros((int((xmax - xmin)/ bstep)))
for i in np.arange(len(bins)):
bins[i] = xmin + i*bstep
xes = np.zeros((len(bins))); inb = np.zeros((len(bins)))
for i in np.arange(len(xaxis)):
for j in np.arange(len(bins)):
if (xaxis[i] > (bins[j] - bstep/2)) & (xaxis[i] <= (bins[j] + bstep/2)):
xes[j] = raw[i] + xes[j]; inb[j] = inb[j] + 1
for i in np.arange(len(inb)):
if inb[i] == 0:
xes[i] = xes[i-1]; inb[i] = inb[i-1]
xes = xes / inb
spline = interp1d(bins, xes, kind='cubic')
xnew = np.arange(Emin-25*xstep,Emax+25*xstep,xstep)
return xnew,spline(xnew)
################
Emin,Emax = 6483,6497 # for area normalization
bstep = 0.75; xstep = 0.01 # for smoothing
plt.figure(figsize=(12,4))
################
plt.subplot(131)
plt.title('Raw')
imin,imax = 0,0
for i in np.arange(len(xaxis)):
if xaxis[i] < Emin: imin = i-1
elif xaxis[i] < Emax: imax = i+1
plt.plot(xaxis,spec0F/sum(spec0F[imin:imax])); plt.plot(xaxis,spec1F/sum(spec1F[imin:imax]))
plt.plot(xaxis,spec2F/sum(spec2F[imin:imax])); plt.plot(xaxis,spec3F/sum(spec3F[imin:imax]))
plt.plot(xaxis,3.0*(spec2F/sum(spec2F[imin:imax])-spec0F/sum(spec0F[imin:imax])))
plt.plot([Emin,Emax],[0,0],'k--')
plt.legend(('0F','1F','2F','3F','2F - 0F (x3)'),loc='upper left',fontsize=10)
Imax = np.max(spec0F/sum(spec0F[imin:imax]))
plt.xlim((Emin,Emax)); plt.ylim((-0.35*Imax,1.1*Imax))
plt.subplot(133)
plt.title('0F')
plt.plot(xaxis,(imax-imin)*spec0F/sum(spec0F[imin:imax]),'b')
################
plt.subplot(132)
plt.title('Smooth')
x0F,y0F = smooth_spec(xaxis,spec0F,bstep,xstep); x1F,y1F = smooth_spec(xaxis,spec1F,bstep,xstep)
x2F,y2F = smooth_spec(xaxis,spec2F,bstep,xstep); x3F,y3F = smooth_spec(xaxis,spec3F,bstep,xstep)
for i in np.arange(len(x0F)):
if x0F[i] < Emin: imin = i-1
elif x0F[i] < Emax: imax = i+1
plt.plot(x0F,y0F/sum(y0F[imin:imax])); plt.plot(x1F,y1F/sum(y1F[imin:imax]))
plt.plot(x2F,y2F/sum(y2F[imin:imax])); plt.plot(x3F,y3F/sum(y3F[imin:imax]))
plt.plot(x0F,3.0*(y2F/sum(y2F[imin:imax])-y0F/sum(y0F[imin:imax])))
plt.plot([Emin,Emax],[0,0],'k--')
Imax = np.max(y0F/sum(y0F[imin:imax]))
plt.xlim((Emin,Emax)); plt.ylim((-0.35*Imax,1.1*Imax))
plt.subplot(133)
plt.plot(x0F,(imax-imin)*y0F/sum(y0F[imin:imax]),'r')
plt.xlim((Emin,Emax))
################
plt.tight_layout()
plt.show()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-1-907923ce25f1> in <module>
23 Emin,Emax = 6483,6497 # for area normalization
24 bstep = 0.75; xstep = 0.01 # for smoothing
---> 25 plt.figure(figsize=(12,4))
26 ################
27 plt.subplot(131)
NameError: name 'plt' is not defined
########################
# track oxidation state
########################
def first_moment(xaxis,spec,hitfinder):
Emin,Emax = 6485,6495
for i in np.arange(len(xaxis)):
if xaxis[i] < Emin: imin = i
if xaxis[i] < Emax: imax = i + 1
int_counts = sum(spec[imin:imax])
if hitfinder: # approximate error estimate, hits
sigma = 2.29425*int_counts**(-0.48482)
else: # approximate error estimate, all
sigma = 3.03577*int_counts**(-0.501203)
xi = np.arange(min(xaxis),max(xaxis),0.01)
spline = interp1d(xaxis, spec, kind='cubic')
yi = spline(xi)
for i in np.arange(len(xi)):
if xi[i] < Emin: imin = i
if xi[i] < Emax: imax = i + 1
return sum(xi[imin:imax]*yi[imin:imax])/sum(yi[imin:imax]),sigma
def peak_max(xaxis,spec):
bstep = 0.75; xstep = 0.001
x,y = smooth_spec(xaxis,spec,bstep,xstep)
mind = np.where(y == np.amax(y))[0][0]
return x[mind]
def peak_width(xaxis,spec):
bstep = 0.75; xstep = 0.001
x,y = smooth_spec(xaxis,spec,bstep,xstep)
mind = np.where(y == np.amax(y))[0][0]
tmp_spec = np.abs(y/max(y)-0.5)
mleft = np.where(tmp_spec[:mind] == np.amin(tmp_spec[:mind]))[0][0]
mright = np.where(tmp_spec[mind:] == np.amin(tmp_spec[mind:]))[0][0] + mind
fwhm = x[mright] - x[mleft]
return fwhm
def int_abs_diff(xaxis,spec,spec_ref):
Emin,Emax = 6485,6495
for i in np.arange(len(xaxis)):
if xaxis[i] < Emin: imin = i
if xaxis[i] < Emax: imax = i + 1
nspec = spec / sum(spec[imin:imax])
nspec_ref = spec_ref / sum(spec_ref[imin:imax])
return sum(np.abs(nspec[imin:imax]-nspec_ref[imin:imax]))
########################
plt.figure(figsize=(10,6))
P = [0,1000,1050,1150,1250,1400,1730,3000,3050,3150,3250,4050,4230,5000]
states = [spec0F,spec1F,spec1F50,spec1F150,spec1F250,spec1F400,spec1F730,
spec2F,spec2F50,spec2F150,spec2F250,spec2F550,spec2F730,spec3F]
fm = []; pmax = []; fwhm = []; iad = []
for state in states:
fm.append(first_moment(xaxis,state,False)[0])
pmax.append(peak_max(xaxis,state))
fwhm.append(peak_width(xaxis,state))
iad.append(int_abs_diff(xaxis,state,states[0]))
plt.subplot(221)
plt.plot(P,fm,'bs--')
plt.xlim((-250,5250)); plt.xticks([0,1000,1500,3000,3500,5000],('0F','1F','+500','2F','+500','3F'))
plt.plot([0,1000,3000,5000],[6490.8234,6490.7626,6490.7369,6490.8053],'mo') # plot previous FM:s
plt.ylim((min(fm)-0.01,max(fm)+0.01))
plt.ylabel('First moment [eV]')
plt.subplot(222)
plt.plot(P,pmax,'bs--')
plt.xlim((-250,5250)); plt.xticks([0,1000,1500,3000,3500,5000],('0F','1F','+500','2F','+500','3F'))
plt.ylabel('Peak max [eV]')
plt.subplot(223)
plt.plot(P,fwhm,'bs--')
plt.xlim((-250,5250)); plt.xticks([0,1000,1500,3000,3500,5000],('0F','1F','+500','2F','+500','3F'))
plt.ylabel('FWHM [eV]')
plt.subplot(224)
plt.plot(P[1:],iad[1:],'bs--')
plt.xlim((-250,5250)); plt.xticks([0,1000,1500,3000,3500,5000],('0F','1F','+500','2F','+500','3F'))
plt.ylabel('Integrated absolute difference')
########################
plt.tight_layout()
##########################
# bootstrap error estimate
##########################
def fm_bootstrap(niter,runs,roi,fit,thr,adu_int,xaxis):
if thr == False:
x,y,adu,ndroplets,gmd = extract(runs)
else:
x,y,adu,ndroplets,rej = HitFinder(runs,thr,roi,fit,adu_int,False)
img,vmin,vmax = create_img(x,y,adu,adu_int,False)
spec = np.sum(remove_background(img, fit[0], fit[1], fit[2], fit[3])[:,roi[0]:roi[1]], axis=1)
interpolate_epix_gap(spec); tot_fm,sd = first_moment(xaxis,spec,True)
all_fm = []
for j in np.arange(niter):
nimages = len(ndroplets); indx = 0
x1 = []; y1 = []; adu1 = []; ndrop1 = []; x2 = []; y2 = []; adu2 = []; ndrop2 = []
for i in np.arange(nimages-1):
S = random.random()
if (S < 0.5):
x1.append(x[indx:indx + ndroplets[i]]); y1.append(y[indx:indx + ndroplets[i]])
adu1.append(adu[indx:indx + ndroplets[i]]); ndrop1.append(ndroplets[i])
else:
x2.append(x[indx:indx + ndroplets[i]]); y2.append(y[indx:indx + ndroplets[i]])
adu2.append(adu[indx:indx + ndroplets[i]]); ndrop2.append(ndroplets[i])
indx += ndroplets[i]
x1 = np.concatenate(x1); y1 = np.concatenate(y1); i1 = np.concatenate(adu1); n1 = ndrop1
x2 = np.concatenate(x2); y2 = np.concatenate(y2); i2 = np.concatenate(adu2); n2 = ndrop2
img,vmin,vmax = create_img(x1,y1,i1,adu_int,False)
spec = np.sum(remove_background(img, fit[0], fit[1], fit[2], fit[3])[:,roi[0]:roi[1]], axis=1)
interpolate_epix_gap(spec); fm,sigma = first_moment(xaxis,spec,True); all_fm.append(fm)
img,vmin,vmax = create_img(x2,y2,i2,adu_int,False)
spec = np.sum(remove_background(img, fit[0], fit[1], fit[2], fit[3])[:,roi[0]:roi[1]], axis=1)
interpolate_epix_gap(spec); fm,sigma = first_moment(xaxis,spec,True); all_fm.append(fm)
print(f'{j+1} / {niter}\r', end="")
avg_fm = np.mean(all_fm); sd = np.std(all_fm)
print('\nTotal / Mean +/- SD\n',np.around(tot_fm,4),' ',np.around(avg_fm,4),' ',np.around(sd,4))
return tot_fm,sd
def plot_fms(p,fm,sd,c):
plt.plot([1,2,3,4],fm,'s--',color=c)
for i in np.arange(len(fm)):
plt.plot([p[i],p[i]],[fm[i]-sd[i],fm[i]+sd[i]],'-',color=c)
plt.plot([p[i]-0.1,p[i]+0.1],[fm[i]-sd[i],fm[i]-sd[i]],'-',color=c)
plt.plot([p[i]-0.1,p[i]+0.1],[fm[i]+sd[i],fm[i]+sd[i]],'-',color=c)
return False
#########################
niter = 20
proper_fm, proper_sd = [],[]
fm,sd = fm_bootstrap(niter,runs_0F,roi,fit,thr,adu_xes,xaxis); proper_fm.append(fm); proper_sd.append(sd)
fm,sd = fm_bootstrap(niter,runs_1F,roi,fit,thr,adu_xes,xaxis); proper_fm.append(fm); proper_sd.append(sd)
fm,sd = fm_bootstrap(niter,runs_2F,roi,fit,thr,adu_xes,xaxis); proper_fm.append(fm); proper_sd.append(sd)
fm,sd = fm_bootstrap(niter,runs_3F,roi,fit,thr,adu_xes,xaxis); proper_fm.append(fm); proper_sd.append(sd)
plt.figure(figsize=(6,4))
#plot_fms([1,2,3,4],approx_fm,approx_sd,'b')
plot_fms([1,2,3,4],proper_fm,proper_sd,'r')
plt.legend(('Bootstrapping error',),loc='upper right',fontsize=9)
plt.xticks([1,2,3,4],('0F','1F','2F','3F'))
plt.ylim((proper_fm[2]-0.01,proper_fm[0]+0.01)); plt.xlim((0.8,4.2))
plt.tight_layout()
plt.show()