Sorting for sample hits¶
Sort for hits: test different schemes and compare hit rate, projections, final spectra
sort away weak pulses (gmd)
sort away wak hits (total xes signal)
can sort with Crystal hits (not shown here)
sorr with our scheme
Sorting using XRD¶
barfoo
Sorting by photon count¶
Two thresholds
def HitFinder(runs,thr,roi,fit,adu_int,print_individual):
x = []; y = []; adu = []; ndroplets = []
rejx = []; rejy = []; rejadu = []; rejndroplets = []
thrrat = thr[1] * (roi[1]-roi[0])/(fit[1]-fit[0]+fit[3]-fit[2])
tot_hits,tot_imgs = 0,0
for run in runs:
tmp_x,tmp_y,tmp_adu,tmp_ndroplets,tmp_gmd = load_run(run)
nimages = len(tmp_ndroplets)
nhit = 0; nmiss = 0; indx = 0
for i in np.arange(nimages-1):
nroi_tmp = 0; n_tmp = 0
for j in np.arange(indx,indx + tmp_ndroplets[i]):
if ((tmp_adu[j] >= adu_int[0]) & (tmp_adu[j] <= adu_int[1])):
if((tmp_y[j] >= roi[0]) & (tmp_y[j] <= roi[1])): nroi_tmp += 1
elif((tmp_y[j] >= fit[0]) & (tmp_y[j] <= fit[1])): n_tmp += 1
elif((tmp_y[j] >= fit[2]) & (tmp_y[j] <= fit[3])): n_tmp += 1
if n_tmp == 0.0: rat = float(nroi_tmp)
else: rat = float(nroi_tmp)/float(n_tmp)
if (rat >= thrrat) and (nroi_tmp >= thr[0]):
x.append(tmp_x[indx:indx + tmp_ndroplets[i]]); y.append(tmp_y[indx:indx + tmp_ndroplets[i]])
adu.append(tmp_adu[indx:indx + tmp_ndroplets[i]]); ndroplets.append(tmp_ndroplets[i])
nhit += 1
else:
rejx.append(tmp_x[indx:indx + tmp_ndroplets[i]]); rejy.append(tmp_y[indx:indx + tmp_ndroplets[i]])
rejadu.append(tmp_adu[indx:indx + tmp_ndroplets[i]]); rejndroplets.append(tmp_ndroplets[i])
nmiss += 1
indx += tmp_ndroplets[i]
if print_individual == True: print(' Hit rate of run {:d}: {:.4f}'.format(run, float(nhit)/float(nhit+nmiss)))
tot_hits += nhit; tot_imgs += nhit + nmiss
print('Total hit rate: {:.4f}, with totally {:d} runs'.format((float(tot_hits)/float(tot_imgs)),tot_imgs))
x = np.concatenate(x); y = np.concatenate(y); adu = np.concatenate(adu)
if len(rejx) > 0:
rejx = np.concatenate(rejx); rejy = np.concatenate(rejy)
rejadu = np.concatenate(rejadu); rejndroplets = rejndroplets
else: rejx = x; rejy = y; rejadu = 0.01*adu; rejndroplets = ndroplets
return x,y,adu,ndroplets,[rejx,rejy,rejadu,rejndroplets]
def HitFinder(runs,thr,roi,fit,adu_int,print_individual,bt):
# note: does not rotate the image before identifying hits
x = []; y = []; adu = []; ndroplets = []; gmd = []
rejx = []; rejy = []; rejadu = []; rejndroplets = []; rejgmd = []
thrrat = thr[1] * (roi[1]-roi[0])/(fit[1]-fit[0]+fit[3]-fit[2]); nrat = 0.0
hitphot,hitrat,rejphot,rejrat = [],[],[],[]
tot_hits,tot_imgs = 0,0
for run in runs:
tmp_x,tmp_y,tmp_adu,tmp_ndroplets,tmp_gmd = load_runs([run],bt)
nimages = len(tmp_ndroplets)
nhit = 0; nmiss = 0; indx = 0
for i in np.arange(nimages-1):
nroi_tmp = 0; n_tmp = 0
for j in np.arange(indx,indx + tmp_ndroplets[i]):
if ((tmp_adu[j] >= adu_int[0]) & (tmp_adu[j] <= adu_int[1])):
if((tmp_y[j] >= roi[0]) & (tmp_y[j] <= roi[1])): nroi_tmp += 1
elif((tmp_y[j] >= fit[0]) & (tmp_y[j] <= fit[1])): n_tmp += 1
elif((tmp_y[j] >= fit[2]) & (tmp_y[j] <= fit[3])): n_tmp += 1
if n_tmp == 0.0: rat = float(nroi_tmp); nrat += 1
else: rat = float(nroi_tmp)/float(n_tmp)
if (rat >= thrrat) and (nroi_tmp >= thr[0]):
hitphot.append(nroi_tmp); hitrat.append(rat*(fit[1]-fit[0]+fit[3]-fit[2])/(roi[1]-roi[0]))
x.append(tmp_x[indx:indx + tmp_ndroplets[i]]); y.append(tmp_y[indx:indx + tmp_ndroplets[i]])
adu.append(tmp_adu[indx:indx + tmp_ndroplets[i]]); ndroplets.append(tmp_ndroplets[i]); gmd.append(tmp_gmd[i])
nhit += 1
else:
rejphot.append(nroi_tmp); rejrat.append(rat*(fit[1]-fit[0]+fit[3]-fit[2])/(roi[1]-roi[0]))
rejx.append(tmp_x[indx:indx + tmp_ndroplets[i]]); rejy.append(tmp_y[indx:indx + tmp_ndroplets[i]])
rejadu.append(tmp_adu[indx:indx + tmp_ndroplets[i]]); rejndroplets.append(tmp_ndroplets[i]); rejgmd.append(tmp_gmd[i])
nmiss += 1
indx += tmp_ndroplets[i]
if print_individual == True: print(' Hit rate of run {:d}: {:.4f}'.format(run, float(nhit)/float(nhit+nmiss)))
tot_hits += nhit; tot_imgs += nhit + nmiss
print 'Hit rate: {:.4f}, including {:d} / {:d},'.format((float(tot_hits)/float(tot_imgs)),tot_hits,tot_imgs),'runs',runs[0],'-',runs[len(runs)-1]
x = np.concatenate(x); y = np.concatenate(y); adu = np.concatenate(adu)
if len(rejx) > 0:
rejx = np.concatenate(rejx); rejy = np.concatenate(rejy)
rejadu = np.concatenate(rejadu); rejndroplets = rejndroplets
else: rejx = x; rejy = y; rejadu = 0.01*adu; rejndroplets = ndroplets; rejgmd = gmd
return x,y,adu,ndroplets,gmd,[rejx,rejy,rejadu,rejndroplets,rejgmd],[hitphot,hitrat,rejphot,rejrat]
def compare_hits(xaxis,runs,thr,roi,fit,adu_xes,print_rate,plot_spec):
allx,ally,alladu,allndroplets,gmd = extract(runs)
x,y,adu,ndroplets,rej = HitFinder(runs,thr,roi,fit,adu_xes,print_rate)
img,vmin,vmax = create_img(x,y,adu,adu_xes,False)
allimg,vmin,vmax = create_img(allx,ally,alladu,adu_xes,False)
rejimg,vmin,vmax = create_img(rej[0],rej[1],rej[2],adu_xes,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)
allspec = np.sum(remove_background(allimg, fit[0], fit[1], fit[2], fit[3])[:,roi[0]:roi[1]], axis=1)
interpolate_epix_gap(allspec)
rejspec = np.sum(remove_background(rejimg, fit[0], fit[1], fit[2], fit[3])[:,roi[0]:roi[1]], axis=1)
interpolate_epix_gap(rejspec)
if plot_spec == True:
plt.plot(xaxis,spec) ; plt.plot(xaxis,allspec)
plt.plot(xaxis,rejspec); plt.plot(xaxis,allspec-spec)
plt.ylim((0,1.05*max(allspec))); plt.xlim((6480,6500))
return spec
##################
thr = [3,2.0] # hit thresholds: total hits in ROI and ratio
plt.figure(figsize=(8,6))
plt.subplot(221); plt.title('0F'); spec0F = compare_hits(xaxis,runs_0F,thr,roi,fit,adu_xes,True,True)
plt.subplot(222); plt.title('1F'); spec1F = compare_hits(xaxis,runs_1F,thr,roi,fit,adu_xes,True,True)
plt.subplot(223); plt.title('2F'); spec2F = compare_hits(xaxis,runs_2F,thr,roi,fit,adu_xes,True,True)
plt.subplot(224); plt.title('3F'); spec3F = compare_hits(xaxis,runs_3F,thr,roi,fit,adu_xes,True,True)
plt.subplot(221); plt.legend(('Hits','All','Misses','All - hits'),fontsize=9,loc='upper left')
plt.tight_layout()
plt.show()
def proj_spectrum(runs,roi,fit,adu_int,axis_para,thr,plot_proj,bt):
if axis_para == False:
xaxis = False
else:
xaxis = energy_axis(axis_para[0],axis_para[1],bt)
if thr == False:
x,y,adu,ndroplets,gmd = load_runs(runs,bt)
print 'Looking at data from beamtime',bt,'with',len(ndroplets),'shots.'
img,vmin,vmax = create_img(runs,x,y,adu,adu_int,adu_int[2],bt)
if bt == 'LV30':
img[:400,:] = 0.0
elif bt == 'LS10':
img[550:,:] = 0.0; img[:,:150] = 0.0; img[:,250:] = 0.0
cor_spec = np.sum(remove_background(img, fit[0], fit[1], fit[2], fit[3])[:,roi[0]:roi[1]], axis=1)
interpolate_epix_gap(cor_spec,bt)
if plot_proj:
tmp_spat = np.sum(img, axis=0)
cor_spat = np.sum(remove_background(img, fit[0], fit[1], fit[2], fit[3]), axis=0)
tmp_spec = np.sum(img[:,roi[0]:roi[1]], axis=1)
else:
x,y,adu,ndroplets,gmd,rej,hitinfo = HitFinder(runs,thr,roi,fit,adu_int,False,bt)
img,vmin,vmax = create_img(runs,x,y,adu,adu_int,adu_int[2],bt)
if bt == 'LV30':
img[:400,:] = 0.0
elif bt == 'LS10':
img[550:,:] = 0.0; img[:,:150] = 0.0; img[:,250:] = 0.0
cor_spec = np.sum(remove_background(img, fit[0], fit[1], fit[2], fit[3])[:,roi[0]:roi[1]], axis=1)
interpolate_epix_gap(cor_spec,bt)
if plot_proj:
x_all,y_all,adu_all,ndroplets_all,gmd_all = load_runs(runs,bt)
img_all,vmin,vmax = create_img(runs,x_all,y_all,adu_all,adu_int,adu_int[2],bt)
cor_spat = np.sum(img, axis=0)
tmp_spat = np.sum(img_all, axis=0)
tmp_spec = np.sum(remove_background(img_all, fit[0], fit[1], fit[2], fit[3])[:,roi[0]:roi[1]], axis=1)
interpolate_epix_gap(tmp_spec,bt)
if plot_proj:
fig = plt.figure(figsize=(12,3))
################
plt.subplot(141); plt.title('Detector image',fontsize=9)
plt.imshow(img, origin='lower',vmin=vmin, vmax=vmax,aspect='auto')
plt.plot([roi[0],roi[0]],[roi[2],roi[3]],'w-'); plt.plot([roi[1],roi[1]],[roi[2],roi[3]],'w-')
plt.plot([fit[0],fit[0]],[roi[2],roi[3]],'w--');plt.plot([fit[1],fit[1]],[roi[2],roi[3]],'w--')
plt.plot([fit[2],fit[2]],[roi[2],roi[3]],'w--');plt.plot([fit[3],fit[3]],[roi[2],roi[3]],'w--')
plt.xlim((fit[0]-20,fit[3]+20))
plt.ylim((roi[2],roi[3]))
plt.xticks(fontsize=8); plt.yticks(fontsize=9)
################
plt.subplot(142); plt.title('Spatial projection',fontsize=9)
plt.plot(tmp_spat/len(ndroplets))
plt.plot(cor_spat/len(ndroplets),'r')
plt.plot((tmp_spat-cor_spat)/len(ndroplets),'k')
ylims = [0,0.8*max(tmp_spat)/len(ndroplets)]
plt.plot([roi[0],roi[0]],ylims,'g'); plt.plot([roi[1],roi[1]],ylims,'g')
plt.plot([fit[0],fit[0]],ylims,'g--');plt.plot([fit[1],fit[1]],ylims,'g--')
plt.plot([fit[2],fit[2]],ylims,'g--');plt.plot([fit[3],fit[3]],ylims,'g--')
plt.xlim((fit[0]-20,fit[3]+20))
plt.ylim((0,1.05*max(tmp_spat)/len(ndroplets)))
plt.xticks(fontsize=9); plt.yticks(fontsize=9)
plt.ylabel('Signal / no. shots',fontsize=9)
################
plt.subplot(143); plt.title('Spectrum',fontsize=9)
if axis_para == False:
plt.plot(tmp_spec/len(ndroplets))
plt.plot(cor_spec/len(ndroplets),'r')
plt.plot((tmp_spec-cor_spec)/len(ndroplets),'k')
else:
plt.plot(xaxis,tmp_spec/len(ndroplets))
plt.plot(xaxis,cor_spec/len(ndroplets),'r')
plt.plot(xaxis,(tmp_spec-cor_spec)/len(ndroplets),'k')
plt.xlim(6480,6500)
if thr == False: plt.legend(('Raw','Corr.','Bkgrd'),fontsize=9,loc='upper left')
else: plt.legend(('All','Hits','Diff.'),fontsize=9,loc='upper left')
plt.ylabel('Signal / no. shots',fontsize=9)
plt.ylim((0,1.05*max(tmp_spec)/len(ndroplets)))
plt.xticks(fontsize=9); plt.yticks(fontsize=9)
################
plt.subplot(144); plt.title('Final spectrum',fontsize=9)
if axis_para == False:
plt.plot(cor_spec,'r')
else:
plt.plot(xaxis,cor_spec,'r')
plt.xlim(6480,6500)
plt.ylabel('Total intensity [counts]',fontsize=9)
plt.ylim((0,1.05*max(cor_spec)))
plt.xticks(fontsize=8); plt.yticks(fontsize=9)
################
plt.tight_layout()
return xaxis,cor_spec
File "<ipython-input-1-f37cf2ebdaf2>", line 71
print 'Hit rate: {:.4f}, including {:d} / {:d},'.format((float(tot_hits)/float(tot_imgs)),tot_hits,tot_imgs),'runs',runs[0],'-',runs[len(runs)-1]
^
SyntaxError: invalid syntax
def first_moment(xaxis,spec,hitfinder):
#Emin,Emax = 6485,6495 #; print 'Thomas'
Emin,Emax = 6485.5,6495.5
#Emin,Emax = 6486.5,6496.5 ; print 'Uwe'
#Emin,Emax = 6486,6500 ; print 'Clemens'
#Emin,Emax = 6483,6500 ; print 'Clemens2'
#print 'Energy interval',Emin,'-',Emax
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 fm_bootstrapdata(niter,runs,x,y,adu,ndroplets,gmd,roi,fit,thr,adu_int,axis_para,bt):
xaxis = energy_axis(axis_para[0],axis_para[1],bt)
img,vmin,vmax = create_img(runs,x,y,adu,adu_int,False,bt)
spec = np.sum(remove_background(img, fit[0], fit[1], fit[2], fit[3])[:,roi[0]:roi[1]], axis=1)
interpolate_epix_gap(spec,bt); 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(runs,x1,y1,i1,adu_int,False,bt)
spec = np.sum(remove_background(img, fit[0], fit[1], fit[2], fit[3])[:,roi[0]:roi[1]], axis=1)
interpolate_epix_gap(spec,bt); fm,sigma = first_moment(xaxis,spec,True); all_fm.append(fm)
img,vmin,vmax = create_img(runs,x2,y2,i2,adu_int,False,bt)
spec = np.sum(remove_background(img, fit[0], fit[1], fit[2], fit[3])[:,roi[0]:roi[1]], axis=1)
interpolate_epix_gap(spec,bt); 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')
#print(np.around(tot_fm,4),' ',np.around(avg_fm,4),' ',np.around(sd,4))
return tot_fm,sd
def fm_bootstrap(niter,runs,roi,fit,thr,adu_int,axis_para,bt):
xaxis = energy_axis(axis_para[0],axis_para[1],bt)
if thr == False:
x,y,adu,ndroplets,gmd = load_runs(runs,bt)
else:
x,y,adu,ndroplets,gmd,rej,hitinfo = HitFinder(runs,thr,roi,fit,adu_int,False,bt)
img,vmin,vmax = create_img(runs,x,y,adu,adu_int,False,bt)
spec = np.sum(remove_background(img, fit[0], fit[1], fit[2], fit[3])[:,roi[0]:roi[1]], axis=1)
interpolate_epix_gap(spec,bt); 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(runs,x1,y1,i1,adu_int,False,bt)
spec = np.sum(remove_background(img, fit[0], fit[1], fit[2], fit[3])[:,roi[0]:roi[1]], axis=1)
interpolate_epix_gap(spec,bt); fm,sigma = first_moment(xaxis,spec,True); all_fm.append(fm)
img,vmin,vmax = create_img(runs,x2,y2,i2,adu_int,False,bt)
spec = np.sum(remove_background(img, fit[0], fit[1], fit[2], fit[3])[:,roi[0]:roi[1]], axis=1)
interpolate_epix_gap(spec,bt); 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')
print(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
def smooth_spec(xaxis,raw,bstep,xstep,Emin,Emax):
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)
def find_indx(x,Emin,Emax):
for i in np.arange(len(x)):
if x[i] < Emin: imin = i
if x[i] < Emax: imax = i + 1
return imin,imax
def intensity_spec(runs,hitthr,intthr,roi,fit,adu_int,axis_para,bt):
xaxis = energy_axis(axis_para[0],axis_para[1],bt)
if hitthr == False:
x,y,adu,ndrop,gmd = load_runs(runs,bt)
hitinfo = False
else:
x,y,adu,ndrop,gmd,rej,hitinfo = HitFinder(runs,hitthr,roi,fit,adu_int,False,bt)
x0,y0,adu0,ndrop0,gmd0,hits0 = sort_intensity(x,y,adu,ndrop,gmd,intthr,hitinfo)
print('Total number of shots:',len(gmd0))
img,vmin,vmax = create_img(runs,x0,y0,adu0,adu_int,adu_int[2],bt)
spec = np.sum(remove_background(img, fit[0], fit[1], fit[2], fit[3])[:,roi[0]:roi[1]], axis=1)
raw_spec = np.sum(img[:,roi[0]:roi[1]], axis=1)
spat = np.sum(img, axis=0)
interpolate_epix_gap(spec,bt)
if hitthr != False:
gmd_com = np.sum(np.array(hits0)*np.array(gmd0))/np.sum(np.array(hits0))
else:
gmd_com = False
fm, sd = first_moment(xaxis,spec,True)
return spec,raw_spec,spat,fm,sd,gmd0,gmd_com,hits0
def sort_intensity(x,y,adu,ndrop,gmd,thr,hitinfo):
x0,y0,adu0,ndrop0,gmd0,indx = [],[],[],[],[],0
if hitinfo == False:
hits0 = False
for i in np.arange(len(gmd)):
if (gmd[i] >= thr[0]) and (gmd[i] <= thr[1]):
x0.append(x[indx:indx + ndrop[i]]); y0.append(y[indx:indx + ndrop[i]])
adu0.append(adu[indx:indx + ndrop[i]]); ndrop0.append(ndrop[i]); gmd0.append(gmd[i])
indx += ndrop[i]
else:
hits0 = []
for i in np.arange(len(gmd)):
if (gmd[i] >= thr[0]) and (gmd[i] <= thr[1]):
x0.append(x[indx:indx + ndrop[i]]); y0.append(y[indx:indx + ndrop[i]])
adu0.append(adu[indx:indx + ndrop[i]]); ndrop0.append(ndrop[i]); gmd0.append(gmd[i]); hits0.append(hitinfo[0][i])
indx += ndrop[i]
x0 = np.concatenate(x0); y0 = np.concatenate(y0); adu0 = np.concatenate(adu0)
return x0,y0,adu0,ndrop0,gmd0,hits0
def intensity_fm(niter,runs,hitthr,intthr,roi,fit,adu_int,axis_para,bt):
xaxis = energy_axis(axis_para[0],axis_para[1],bt)
if hitthr == False:
x,y,adu,ndrop,gmd = xest.load_runs(runs,bt)
hitinfo = False
else:
x,y,adu,ndrop,gmd,rej,hitinfo = HitFinder(runs,hitthr,roi,fit,adu_int,False,bt)
x0,y0,adu0,ndrop0,gmd0,hits0 = sort_intensity(x,y,adu,ndrop,gmd,intthr,hitinfo)
#
img,vmin,vmax = create_img(runs,x0,y0,adu0,adu_int,adu_int[2],bt)
spec = np.sum(remove_background(img, fit[0], fit[1], fit[2], fit[3])[:,roi[0]:roi[1]], axis=1)
interpolate_epix_gap(spec,bt)
if hitthr != False:
gmd_com = np.sum(np.array(hits0)*np.array(gmd0))/np.sum(np.array(hits0))
else:
gmd_com = False
fm, sd = first_moment(xaxis,spec,True)
print('Center-of-mass intensity:',np.around(gmd_com,4))
print('First moment and approx. error',np.around(fm,4),np.around(sd,4))
tot_fm,sd = fm_bootstrapdata(niter,runs,x0,y0,adu0,ndrop0,gmd0,roi,fit,thr,adu_int,axis_para,bt)
return tot_fm,sd,gmd_com
def proj_spectrum(runs,roi,fit,adu_int,axis_para,thr,plot_proj,bt):
xaxis = xest.energy_axis(axis_para[0],axis_para[1],bt)
x,y,adu,ndroplets,gmd,rej,hitinfo = xest.HitFinder(runs,thr,roi,fit,adu_int,False,bt)
img,vmin,vmax = create_img(runs,x,y,adu,adu_int,adu_int[2],bt)
cor_spec = np.sum(xest.remove_background(img, fit[0], fit[1], fit[2], fit[3])[:,roi[0]:roi[1]], axis=1)
xest.interpolate_epix_gap(cor_spec,bt)
return xaxis,cor_spec
def plot_fms(p,fm,sd,c):
for i in np.arange(len(fm)):
plt.plot(p[i],fm[i],'s',color=c)
plt.plot([p[i],p[i]],[fm[i]-sd[i],fm[i]+sd[i]],'-',color=c)
plt.plot([p[i]-40,p[i]+40],[fm[i]-sd[i],fm[i]-sd[i]],'-',color=c)
plt.plot([p[i]-40,p[i]+40],[fm[i]+sd[i],fm[i]+sd[i]],'-',color=c)
return False
# Spectra, check of different runs
bt = 'LV68'; roi,fit,thr,axis_para,adu_int = xest.load_bt_para(bt)
roi5 = [495,565,0,709]; fit5 = [roi5[0]-60,roi5[0]-25,roi5[1]+20,roi5[1]+55]; roi,fit = roi5,fit5
runs = [244]; specA = xest.proj_spectrum(runs,roi,fit,adu_int,axis_para,thr,False,bt)
runs = [268]; specB = xest.proj_spectrum(runs,roi,fit,adu_int,axis_para,thr,False,bt)
imin,imax = xest.find_indx(specA[0],6484,6496)
plt.figure(figsize=(6,3))
tmp = specA; plt.plot(tmp[0],tmp[1]/sum(tmp[1][imin:imax]),label='244')
tmp = specB; plt.plot(tmp[0],tmp[1]/sum(tmp[1][imin:imax]),label='268')
plt.legend(loc='upper left',fontsize=8); plt.xlim((6480,6500)); plt.tight_layout(); plt.show()
import scipy.signal
def count_photons(xaxis,spec):
imin,imax = xest.find_indx(xaxis,6485,6495)
return sum(spec[imin:imax])
def smooth(y1,y2,imin,imax):
y1i = scipy.signal.savgol_filter(y1, 51, 2)
y2i = scipy.signal.savgol_filter(y2, 51, 2)
return y1i/sum(y1i[imin:imax])-y2i/sum(y2i[imin:imax])
def first_moment(xaxis,spec,hitfinder,Emin,Emax):
print 'Energy interval',Emin,'-',Emax
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