diff --git a/dispersiveXanes_alignment.py b/dispersiveXanes_alignment.py index cca56df..e6b84a7 100644 --- a/dispersiveXanes_alignment.py +++ b/dispersiveXanes_alignment.py @@ -18,6 +18,7 @@ import dispersiveXanes_utils as utils #defaultE = np.arange(1024)*0.12+7060 defaultE = (np.arange(1024)-512)*0.189+7123 defaultE = 7063.5+0.1295*np.arange(1024); # done on nov 30 2016, based on xppl3716:r77 +defaultE = 7064.5-6+0.1295*np.arange(1024); # done on jun 07 2017, based on reference spectrum from ESRF __i = np.arange(1024) __x = (__i-512)/512 diff --git a/dispersiveXanes_utils.py b/dispersiveXanes_utils.py index 10467d7..e90211e 100644 --- a/dispersiveXanes_utils.py +++ b/dispersiveXanes_utils.py @@ -47,6 +47,8 @@ def maskLowIntensity(p1,p2,threshold=0.03,squeeze=True): if squeeze: p1 = np.squeeze(p1); p2 = np.squeeze(p2) + p1.fill_value=np.nan + p2.fill_value=np.nan return p1,p2 def ratioOfAverage(p1,p2,threshold=0.03): diff --git a/sigma_vs_shot.py b/sigma_vs_shot.py new file mode 100644 index 0000000..fbddb6d --- /dev/null +++ b/sigma_vs_shot.py @@ -0,0 +1,43 @@ +import matplotlib.pyplot as plt +import numpy as np +import xppl37_spectra +import xanes_analyzeRun +from datastorage import DataStorage as ds +r=xppl37_spectra.xanes_analyzeRun.AnalyzeRun(81) +r.load() +if len(r.results.keys()) > 1: + p1 = np.vstack( ( r.results[c].p1 for c in range(0,10,2) ) ) + p2 = np.vstack( ( r.results[c].p2 for c in range(0,10,2) ) ) +else: + p1 = r.results[0].p1 + p2 = r.results[0].p2 +N = int(len(p1)/2) + +ref = slice(0,1000) +sam = slice(1000,1200) + +#ref = slice(200) +#sam = slice(200) + +ref = ds(p1=p1[ref],p2=p2[ref]) +sam = ds(p1=p1[sam],p2=p2[sam]) + +# "ratioOfAverage medianOfRatios" +data=xppl37_spectra.calcAbs(ref,sam,refKind="medianOfRatios") + +abs = data[-1] +idx = slice(200,800) +abs = abs[:,idx] + +fom = [] +fom1 = [] +x = np.arange(abs.shape[1]) +N = range(1,len(abs)) +#N = 2**np.arange(15) +for i in N: + av = np.nanmean(abs[:i],axis=0) + p = np.polyfit(x,av,5) + diff = av-np.polyval(p,x) + fom.append( np.nanstd(av) ) + fom1.append( np.nanstd(diff) ) + diff --git a/xanes_analyzeRun.py b/xanes_analyzeRun.py index 305185e..2f2345c 100644 --- a/xanes_analyzeRun.py +++ b/xanes_analyzeRun.py @@ -25,9 +25,9 @@ g_exp = "xppl3716" g_bml = g_exp[:3] x3py.config.updateBeamline(g_bml) - -g_folder_init = g_exp+"_init_pars/" -g_folder_out = g_exp+"_output/" +basedir = os.path.dirname(__file__) +g_folder_init = basedir + "/" +g_exp+"_init_pars/" +g_folder_out = basedir + "/" +g_exp+"_output/" g_folder_data = "/reg/d/psdm/"+g_bml+"/"+ g_exp +"/hdf5/" import socket diff --git a/xppl3716_init_pars/run0039_transform.npy b/xppl3716_init_pars/run0039_transform.npy new file mode 100644 index 0000000..218245e Binary files /dev/null and b/xppl3716_init_pars/run0039_transform.npy differ diff --git a/xppl3716_init_pars/run0080_transform.npy b/xppl3716_init_pars/run0080_transform.npy index bd9eea0..532c2cf 100644 Binary files a/xppl3716_init_pars/run0080_transform.npy and b/xppl3716_init_pars/run0080_transform.npy differ diff --git a/xppl3716_init_pars/run0081_transform.npy b/xppl3716_init_pars/run0081_transform.npy new file mode 100644 index 0000000..a9e609b Binary files /dev/null and b/xppl3716_init_pars/run0081_transform.npy differ diff --git a/xppl3716_init_pars/run0167_transform.npy b/xppl3716_init_pars/run0167_transform.npy new file mode 100644 index 0000000..162774e Binary files /dev/null and b/xppl3716_init_pars/run0167_transform.npy differ diff --git a/xppl37_spectra.py b/xppl37_spectra.py index 380fca6..b8af8e3 100644 --- a/xppl37_spectra.py +++ b/xppl37_spectra.py @@ -54,7 +54,8 @@ def calcSpectraForRun(run=82,calibs="all",realign=False,init="auto",alignCalib=0 elif isinstance(calibs,slice): calibsForOut = list(range(r.nCalib))[calibs] elif calibs == "all": - calibsForOut = r.results.keys() + calibsForOut = list(r.results.keys()) + calibsForOut.sort() else: calibsForOut = calibs # focused data have one single calibcycle ... @@ -62,9 +63,10 @@ def calcSpectraForRun(run=82,calibs="all",realign=False,init="auto",alignCalib=0 p1 = [r.results[calib].p1 for calib in calibsForOut] p2 = [r.results[calib].p2 for calib in calibsForOut] else: - idx = r.results[None].fom < 0.5 - p1 = [ r.results[None].p1[idx], r.results[None].p1[~idx] ] - p2 = [ r.results[None].p2[idx], r.results[None].p2[~idx] ] + calib = list(r.results.keys())[0] + idx = r.results[calib].fom < 0.5 + p1 = [ r.results[calib].p1[idx], r.results[calib].p1[~idx] ] + p2 = [ r.results[calib].p2[idx], r.results[calib].p2[~idx] ] return profile_ret( run =r, p1 =p1, p2=p2,calibs=calibsForOut) @@ -111,7 +113,11 @@ def calcSpectraForRefAndSample(run=82,refCalibs=slice(None,None,2),forceSpectraC def calcRef(r1,r2,calibs=None,threshold=0.05): """ r1 and r2 are list of 2d arrays (nShots,nPixels) for each calibcycle """ + if not isinstance(r1,(tuple,list)): r1 = (r1,) + if not isinstance(r2,(tuple,list)): r2 = (r2,) if calibs is None: calibs = list(range(len(r1))) + if not isinstance(calibs,(tuple,list)): calibs = (calibs,) + out = collections.OrderedDict() out["ratioOfAverage"] = dict() out["medianOfRatios"] = dict() @@ -153,21 +159,44 @@ 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",forceSpectraCalculation=False): +def calcAbs(ref,sample=None,threshold=0.05,refKind="medianOfRatios",merge_calibs=True): """ example of use - ratio = calcSampleAbs(82) - ratio = calcSampleAbs( (155,156) ) + ratio = calcAbsForRun(82) + ratio = calcAbsForRun( (155,156) ) + """ + if sample is None: sample=ref + temp = calcRef(ref.p1,ref.p2,threshold=threshold) + ref = temp[refKind]['all'] + if merge_calibs: + p1 = np.vstack( sample.p1 ) + p2 = np.vstack( sample.p2 ) + p1,p2 = utils.maskLowIntensity(p1,p2,threshold=threshold) + ratio = p2/p1 + ratio = ratio/ref + Abs = -np.log10(ratio) + else: + Abs = [] + p1 = [] + p2 = [] + for _p1,_p2 in zip(sample.p1,sample.p2): + _p1,_p2 = utils.maskLowIntensity(_p1,_p2,threshold=threshold) + ratio = _p2/_p1 + ratio = ratio/ref + _abs = -np.log10(ratio) + p1.append(_p1) + p2.append(_p2) + Abs.append(_abs) + return p1,p2,Abs + +def calcAbsForRun(run=82,refCalibs=slice(None,None,2),threshold=0.05,refKind="medianOfRatios",forceSpectraCalculation=False,merge_calibs=True): + """ example of use + ratio = calcAbsForRun(82) + ratio = calcAbsForRun( (155,156) ) """ 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 = utils.maskLowIntensity(p1,p2,threshold=threshold) - ratio = p2/p1 - ratio = ratio/ref - return p1,p2,-np.log10(ratio) + E = ref.run.E + p1,p2,Abs = calcAbs(ref,sample,threshold=threshold,refKind=refKind,merge_calibs=merge_calibs) + return E,p1,p2,Abs def showSpectra(run=82,shots=slice(5),calibs=0,averageEachCalib=False, normalization="auto",shifty=1,xlim=(7060,7180),showAv=True): @@ -209,7 +238,7 @@ def showAbs(run=82,shots=slice(5),normalization="auto",shifty=1,xlim=(7080,7180) filterShot = means that it filters out the filterShot*100 percentile """ E = alignment.defaultE - p1,p2,abs = calcSampleAbs(run=run,threshold=threshold) + _,p1,p2,abs = calcAbsForRun(run=run,threshold=threshold) p1_sum = p1.sum(-1) if filterShot>0: idx = p1_sum>np.percentile(p1_sum,filterShot*100) @@ -252,19 +281,23 @@ def showAbsWithSweep(run=(155,156),first=0,period=150,nSpectra=10,**kwards): shots = slice(first,first+period,int(period/nSpectra)) showAbs(run=run,shots=shots,**kwards) -def smoothSpectra(E,abs_spectra,res=0.5): +def smoothSpectra(E,abs_spectra,res=0.5,skip_close=5): from scipy import integrate + if isinstance(abs_spectra,np.ma.MaskedArray): abs_spectra = abs_spectra.filled() 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]) + for ispectrum,spectrum in enumerate(abs_spectra): + idx = np.isfinite(spectrum) Eclean = E[idx] + spectrum_clean = spectrum[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] + g = 1/np.sqrt(2*np.pi)/res*np.exp(-(Eclean-E[i])**2/2/res**2) + tointegrate = g*spectrum_clean out[ispectrum,i] = integrate.simps(tointegrate,x=Eclean) + out[ispectrum][~idx]=np.nan + temp = out[ispectrum].copy() + for i,o in enumerate(temp): + if np.any( np.isnan(temp[i-skip_close:i+skip_close] ) ): out[ispectrum][i] = np.nan return out def doLongCalc():