mcammara
/
dispersiveXanes
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

improved calculations (calcSpectraForRefAndSample) and plotting (showSpectra and showAbs)

This commit is contained in:
Marco Cammarata 2016-12-01 11:16:46 +01:00
parent 66dd58a94b
commit 527f88b50c
1 changed files with 129 additions and 24 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

153
xppl37_spectra.py
View File

@ -3,6 +3,7 @@ import matplotlib.pyplot as plt
import copy import copy
import argparse import argparse
import collections import collections
import alignment
import xanes_analyzeRun import xanes_analyzeRun
parser = argparse.ArgumentParser(description='Process argv') parser = argparse.ArgumentParser(description='Process argv')
@ -16,8 +17,10 @@ profile_ret = collections.namedtuple("profile_ret",["run","p1","p2","calibs"])
nice_colors = ["#1b9e77", "#d95f02", "#7570b3"] nice_colors = ["#1b9e77", "#d95f02", "#7570b3"]
gradual_colors = ['#014636', '#016c59', '#02818a', '#3690c0', '#67a9cf', '#a6bddb', '#d0d1e6']#, '#ece2f0']
def calcProfilesForRun(run=82,calibs="all",realign=False,init="auto",alignCalib=0):
def calcSpectraForRun(run=82,calibs="all",realign=False,init="auto",alignCalib=0,force=False):
""" init and alignCalib are used only if realignment is performed: """ init and alignCalib are used only if realignment is performed:
init = run for initial alignment, use auto is you want to use same run init = run for initial alignment, use auto is you want to use same run
alignCalib = calibcycle for alignment alignCalib = calibcycle for alignment
@ -32,12 +35,13 @@ def calcProfilesForRun(run=82,calibs="all",realign=False,init="auto",alignCalib=
r.saveTransform() r.saveTransform()
# next line is used to force calculations in case of realignment # next line is used to force calculations in case of realignment
fname = 'auto' if not realign else "thisfiledoesnotexists" fname = 'auto' if not realign else "thisfiledoesnotexists"
if force: fname = "thisfiledoesnotexists"
if len(r.results) == 0: if len(r.results) == 0:
try: try:
r.load(fname) r.load(fname)
print("Loading previously saved results") print("Loading previously saved results")
except FileNotFoundError: except FileNotFoundError:
r.analyzeScan(shots=slice(600),calibs=calibs,nImagesToFit=0,nSaveImg=4) r.analyzeScan(calibs=calibs,nImagesToFit=0,nSaveImg=4)
r.save(overwrite=True) r.save(overwrite=True)
# cannot take the output from r.results because it might have been calculated for # cannot take the output from r.results because it might have been calculated for
# a bigger range than asked for. # a bigger range than asked for.
@ -54,23 +58,41 @@ def calcProfilesForRun(run=82,calibs="all",realign=False,init="auto",alignCalib=
return profile_ret( run =r, p1 =p1, p2=p2,calibs=calibsForOut) return profile_ret( run =r, p1 =p1, p2=p2,calibs=calibsForOut)
def calcProfilesForRefAndSample(run=82,refCalibs=0,force=False): def calcSpectraForRefAndSample(run=82,refCalibs=slice(None,None,2),forceSpectraCalculation=False):
""" Function to calculate the Spectra with and without (ref) the sample.
It can analyze two kinds of runs:
* Single run with that alternates IN and OUT (like run 82)
in this case use something like:
calcSpectraForRefAndSample(82,refCalibs=slice(None,None,2)
or
calcSpectraForRefAndSample(82,refCalibs=(0,2,4,6))
* Multiple runs (one with reference and one with sample) like run 155 and 156
in this case use something like:
calcSpectraForRefAndSample( run=(155,156) )
where the first is the reference run and the second the one with the sample.
in this second case the refCalibs does not play a role
use forceSpectraCalculation = True, to re-read the images
"""
if isinstance(run,int): if isinstance(run,int):
refRun = xanes_analyzeRun.AnalyzeRun(run) run = xanes_analyzeRun.AnalyzeRun(run)
sampleRun = refRun if isinstance(refCalibs,slice): refCalibs = list(range(run.nCalib))[refCalibs]
if isinstance(refCalibs,slice): refCalibs = list(range(refRun.nCalib))[refCalibs] if isinstance(refCalibs,int): refCalibs = [refCalibs,]
if isinstance(refCalibs,int): refCalibs = (refCalibs,)
sampleCalibs = [c+1 for c in refCalibs] sampleCalibs = [c+1 for c in refCalibs]
# need a single call (for sample and ref) to save all calibcycles
calcSpectraForRun(run,calibs=refCalibs+sampleCalibs,\
force=forceSpectraCalculation);
ref = calcSpectraForRun(run,calibs=refCalibs)
sample = calcSpectraForRun(run,calibs=sampleCalibs)
elif isinstance(run,(list,tuple)): elif isinstance(run,(list,tuple)):
refRun = xanes_analyzeRun.AnalyzeRun(run[0]) refRun = xanes_analyzeRun.AnalyzeRun(run[0])
sampleRun = xanes_analyzeRun.AnalyzeRun(run[1],initAlign=run[0]) sampleRun = xanes_analyzeRun.AnalyzeRun(run[1],initAlign=run[0])
refCalibs = [0,] refCalibs = [0,]
sampleCalibs = [0,] sampleCalibs = [0,]
if refRun == sampleRun: ref = calcSpectraForRun(refRun,calibs=refCalibs,\
# otherwise same does not save them both force=forceSpectraCalculation)
temp = calcProfilesForRun(refRun,calibs=sampleCalibs+refCalibs); sample = calcSpectraForRun(sampleRun,calibs=sampleCalibs,\
ref = calcProfilesForRun(refRun,calibs=refCalibs) force=forceSpectraCalculation)
sample = calcProfilesForRun(sampleRun,calibs=sampleCalibs)
return ref,sample return ref,sample
def calcRef(r1,r2,calibs=None,threshold=0.05): def calcRef(r1,r2,calibs=None,threshold=0.05):
@ -96,7 +118,7 @@ def calcRef(r1,r2,calibs=None,threshold=0.05):
def showDifferentRefs(run=82,refCalibs=slice(None,None,2),threshold=0.05): def showDifferentRefs(run=82,refCalibs=slice(None,None,2),threshold=0.05):
""" example plots showing how stable are the different ways of taking spectra """ """ example plots showing how stable are the different ways of taking spectra """
prof = calcProfilesForRun(run,calibs=refCalibs) prof = calcSpectraForRun(run,calibs=refCalibs)
refs = calcRef(prof.p1,prof.p2,calibs=prof.calibs) refs = calcRef(prof.p1,prof.p2,calibs=prof.calibs)
kind_of_av = list(refs.keys()) kind_of_av = list(refs.keys())
fig,ax=plt.subplots(len(kind_of_av)+1,1,sharex=True,sharey=True) fig,ax=plt.subplots(len(kind_of_av)+1,1,sharex=True,sharey=True)
@ -117,29 +139,112 @@ def showDifferentRefs(run=82,refCalibs=slice(None,None,2),threshold=0.05):
ax[-1].legend() ax[-1].legend()
for a in ax: a.grid() for a in ax: a.grid()
def calcSampleAbs(run=82,refCalibs=slice(None,None,2),threshold=0.05,refKind="medianOfRatios"): def calcSampleAbs(run=82,refCalibs=slice(None,None,2),threshold=0.05,refKind="medianOfRatios",forceSpectraCalculation=False):
""" example of use """ example of use
ratio = calcSampleAbs(82) ratio = calcSampleAbs(82)
ratio = calcSampleAbs( (155,156) ) ratio = calcSampleAbs( (155,156) )
""" """
ref,sample = calcProfilesForRefAndSample(run,refCalibs=refCalibs) ref,sample = calcSpectraForRefAndSample(run,refCalibs=refCalibs,forceSpectraCalculation=forceSpectraCalculation)
temp = calcRef(ref.p1,ref.p2,calibs=ref.calibs,threshold=threshold) temp = calcRef(ref.p1,ref.p2,calibs=ref.calibs,threshold=threshold)
ref = temp[refKind]['all'] ref = temp[refKind]['all']
p1 = np.vstack( sample.p1 ) p1 = np.vstack( sample.p1 )
p2 = np.vstack( sample.p2 ) p2 = np.vstack( sample.p2 )
print(p1.shape) print(p1.shape)
p1,p2 = xanes_analyzeRun.maskLowIntensity(p1,p2,threshold=0.1) p1,p2 = xanes_analyzeRun.maskLowIntensity(p1,p2,threshold=threshold)
ratio = p2/p1 ratio = p2/p1
ratio = ratio/ref ratio = ratio/ref
return ratio return p1,p2,-np.log10(ratio)
def main(run,refCalib=0,force=False): def showSpectra(run=82,shots=slice(5),calibs=0,averageEachCalib=False,
pass normalization="auto",shifty=1,xlim=(7060,7180),showAv=True):
#r = calcProfiles(run,refCalibs=refCalib,force=force) """ averageEachCalib: if True, plot only one (averaged) spectrum per calibcycle
normalization: if "auto", the max of the spectra that will be plotted will be used
"""
r = xanes_analyzeRun.AnalyzeRun(run=run)
r.load()
calibsSaved = list(r.results.keys()); calibsSaved.sort()
res = [ r.results[c].p2 for c in calibsSaved ]
if isinstance(calibs,slice): res=res[calibs]
if isinstance(calibs,int): res=[res[calibs],]
avCalibs = [ np.nanmedian(spectra,axis=0) for spectra in res ]
if averageEachCalib:
res = avCalibs
showAv = False; # it does not make sense to plot it twice !
fig,ax = plt.subplots(len(res),1,sharex=True,sharey=True,squeeze=False)
if normalization == "auto": normalization = np.nanmax( [temp[shots] for temp in res] )
for (av,spectra,a) in zip(avCalibs,res,ax[:,0]):
spectra_norm = spectra[shots]/normalization
av_norm = av/normalization
for i,spectrum in enumerate(spectra_norm):
color = gradual_colors[i%len(gradual_colors)]
a.axhline(i*shifty,ls='--',color=color)
if showAv: a.fill_between(r.E,i*shifty,av_norm+i*shifty,color='#d95f0e',alpha=0.4)
print(i)
a.plot(r.E,spectrum+i*shifty,color=color,lw=2)
ax[0][0].set_xlim(*xlim)
ax[0][0].set_title("Run %d"%run)
if not averageEachCalib: ax[0][0].set_ylim(0,shifty*(len(spectra_norm)))
def showAbs(run=82,shots=slice(5),normalization="auto",shifty=1,xlim=(7080,7180),showAv=True,smoothWidth=0):
""" normalization: if "auto", the max of the spectra that will be plotted will be used
"""
E = alignment.defaultE
p1,p2,abs = calcSampleAbs(run=run,threshold=0.01)
p1_av = np.nanmedian(p1,axis=0)
p2_av = np.nanmedian(p2,axis=0)
abs_av = np.nanmedian(abs,axis=0)
p1 = p1[shots]; p2=p2[shots]; abs = abs[shots]
if smoothWidth > 0: abs = smoothSpectra(E,abs,res=smoothWidth)
fig,ax = plt.subplots(1,2,sharex=True,sharey=True,squeeze=False)
ax = ax[0]
if normalization == "auto": normalization = np.nanmax( p1 )
for ishot,(s1,s2,a) in enumerate(zip(p1,p2,abs)):
s1_norm = s1/normalization
s2_norm = s2/normalization
color = gradual_colors[ishot%len(gradual_colors)]
ax[0].axhline(ishot*shifty,ls='--',color=color)
ax[1].axhline(ishot*shifty,ls='--',color=color)
if showAv:
if ishot == 0: ax[0].fill_between(E,ishot*shifty,p1_av/normalization+ishot*shifty,color='#d95f0e',alpha=0.6)
ax[1].plot(E,abs_av+ishot*shifty,color='#d95f0e',lw=2,zorder=20)
ax[0].plot(E,s1_norm+ishot*shifty,ls = '-' ,color='0.8',lw=2)
ax[0].plot(E,s2_norm+ishot*shifty,ls = '-' ,color='0.3',lw=2)
ax[1].plot(E,a+ishot*shifty,color=color,lw=2)
ax[0].set_xlim(*xlim)
ax[0].set_title("Run %d"%run)
ax[1].set_ylabel("Sample Absorption")
ax[0].set_ylabel("Normalized Spectra")
ax[0].set_ylim(0,shifty*(p1.shape[0]))
def smoothSpectra(E,abs_spectra,res=0.5):
from scipy import integrate
if abs_spectra.ndim == 1: abs_spectra=abs_spectra[np.newaxis,:]
out = np.empty_like(abs_spectra)
for ispectrum in range(abs_spectra.shape[0]):
idx = np.isfinite(abs_spectra[ispectrum])
Eclean = E[idx]
for i in range(len(E)):
g = 1/np.sqrt(2*np.pi)/res*np.exp(-(E-E[i])**2/2/res**2)
tointegrate = g*abs_spectra[ispectrum]
# filter out nans
tointegrate = tointegrate[idx]
out[ispectrum,i] = integrate.simps(tointegrate,x=Eclean)
return out
def doLongCalc():
#calcSpectraForRefAndSample(82,forceSpectraCalculation=True)
# scanning requires a lower level call
r = xanes_analyzeRun.AnalyzeRun(84)
r.analyzeScan(nShotsPerCalib="all",calibs="all",nSaveImg=2,calibsToFit='even',nImagesToFit=3)
r.save(overwrite=True)
#calcSpectraForRefAndSample(84,forceSpectraCalculation=False)
calcSpectraForRefAndSample(96,forceSpectraCalculation=True)
calcSpectraForRefAndSample((155,156),forceSpectraCalculation=True)
#r = calcSpectra(run,refCalibs=refCalib,force=force)
if __name__ == "__main__": if __name__ == "__main__":
main(args.run,force=args.force) pass
#main(args.run,force=args.force)