diff --git a/xppl37_spectra.py b/xppl37_spectra.py index a171db5..d079b5c 100644 --- a/xppl37_spectra.py +++ b/xppl37_spectra.py @@ -3,6 +3,7 @@ import matplotlib.pyplot as plt import copy import argparse import collections +import alignment import xanes_analyzeRun 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"] +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 = run for initial alignment, use auto is you want to use same run alignCalib = calibcycle for alignment @@ -32,12 +35,13 @@ def calcProfilesForRun(run=82,calibs="all",realign=False,init="auto",alignCalib= r.saveTransform() # next line is used to force calculations in case of realignment fname = 'auto' if not realign else "thisfiledoesnotexists" + if force: fname = "thisfiledoesnotexists" if len(r.results) == 0: try: r.load(fname) print("Loading previously saved results") 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) # cannot take the output from r.results because it might have been calculated 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) -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): - refRun = xanes_analyzeRun.AnalyzeRun(run) - sampleRun = refRun - if isinstance(refCalibs,slice): refCalibs = list(range(refRun.nCalib))[refCalibs] - if isinstance(refCalibs,int): refCalibs = (refCalibs,) + run = xanes_analyzeRun.AnalyzeRun(run) + if isinstance(refCalibs,slice): refCalibs = list(range(run.nCalib))[refCalibs] + if isinstance(refCalibs,int): refCalibs = [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)): refRun = xanes_analyzeRun.AnalyzeRun(run[0]) sampleRun = xanes_analyzeRun.AnalyzeRun(run[1],initAlign=run[0]) refCalibs = [0,] sampleCalibs = [0,] - if refRun == sampleRun: - # otherwise same does not save them both - temp = calcProfilesForRun(refRun,calibs=sampleCalibs+refCalibs); - ref = calcProfilesForRun(refRun,calibs=refCalibs) - sample = calcProfilesForRun(sampleRun,calibs=sampleCalibs) + ref = calcSpectraForRun(refRun,calibs=refCalibs,\ + force=forceSpectraCalculation) + sample = calcSpectraForRun(sampleRun,calibs=sampleCalibs,\ + force=forceSpectraCalculation) return ref,sample 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): """ 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) kind_of_av = list(refs.keys()) 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() 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 ratio = calcSampleAbs(82) 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) ref = temp[refKind]['all'] p1 = np.vstack( sample.p1 ) p2 = np.vstack( sample.p2 ) 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 = ratio/ref - return ratio - - - + return p1,p2,-np.log10(ratio) -def main(run,refCalib=0,force=False): - pass - #r = calcProfiles(run,refCalibs=refCalib,force=force) +def showSpectra(run=82,shots=slice(5),calibs=0,averageEachCalib=False, + normalization="auto",shifty=1,xlim=(7060,7180),showAv=True): + """ 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__": - main(args.run,force=args.force) + pass + #main(args.run,force=args.force)