From b53354da9b4b3e3b51670f000f47c61ceee0d381 Mon Sep 17 00:00:00 2001 From: Marco Cammarata Date: Fri, 9 Jun 2017 16:42:35 +0200 Subject: [PATCH] end of marion visit (main figures ready --- dispersiveXanes_alignment.py | 1 + dispersiveXanes_utils.py | 2 + sigma_vs_shot.py | 43 ++++++++++++ xanes_analyzeRun.py | 6 +- xppl3716_init_pars/run0039_transform.npy | Bin 0 -> 1703 bytes xppl3716_init_pars/run0080_transform.npy | Bin 647 -> 3417 bytes xppl3716_init_pars/run0081_transform.npy | Bin 0 -> 3416 bytes xppl3716_init_pars/run0167_transform.npy | Bin 0 -> 3413 bytes xppl37_spectra.py | 81 ++++++++++++++++------- 9 files changed, 106 insertions(+), 27 deletions(-) create mode 100644 sigma_vs_shot.py create mode 100644 xppl3716_init_pars/run0039_transform.npy create mode 100644 xppl3716_init_pars/run0081_transform.npy create mode 100644 xppl3716_init_pars/run0167_transform.npy 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 0000000000000000000000000000000000000000..218245e330a68d7c12f8b673a441abd45af0c0c1 GIT binary patch literal 1703 zcmZuxTWl0n7@nzxmZ1wMRZ$TaL3V3PWv5;cl^GFxWXGMtDmo^Z$?WXVj_l6#pP5Z| zg(hfZyZc7O7a$})LsCtQ2~m>`F)=X_jlLKXFNjYKMlkxoll9!XGqj$EcJ}<2^IiUP zPU!~=``-|Sap7pD?1XlZ*_F;5o5+l$Gl%^k3M|j`gR(=m$E`}};QesQns)H-?9fPh zCVl&oj`!GJtvcNpwf(>ut=1}$YXyPTkj%ibeJ_lHnjJ|Xuq??dvnYs~Z~}U4)Aq}b z9ZHTzJ|wuq5CuUfL9ekA|F2k)>*ZiY6IK?XZx;HCuqvyGL%dlzYE0uA$>_nk6FXtG zCTQX;3=|=i)mW!vutwMQ`~qI|kETShw$y$T*3H13S)&&h2ELnvyKK`e)m)6;GfgS1 zxHfLQzVNYFsQ2izFj$1U7n5bUXKH=H5ewT%FCY=i5mMXO%2M-~fLRX?a zkEJpB!6s}r*5L#v2>ifwC#_l-%Aq?64lehj%K+6%(k0SLWm10XkqjY5$W^_kr4V98eB*XTWK_Jl`N}qNbZ7r9D zkbf+eh4{sYaPm8l`tF8*E=`X7k~p1Wu0swfS5A|41}T4!(bWW*=5gtnVPi=7;O*bg z+~CnfjPrgRDVM+a{hL!i&Qa%iU3Va!K*~))9Qbs8o~lR)g;S0dz?1DRXWr?lHf6`A zjai!`>~`u@3%1|2=V{cLyLjOGiA@)g!atg@i#oz5^%Ex~?Ws`|)L97#&*~Eo#)0_mk(ZbwE+RGm^R?WOzJ7XMV%A=XpOHnD{9$^v72ozc97M{E i(pNw2PJL#qz_q3 zd6lm3EOivKNu`qIP(r0s`COr>N~vN90Z%HFYx8yyY!0kK zRVtCnx202sP&y}z5QfBqfk2=YL*v8w-jOV3a}5}#F`Tf^YK$O^tZNQVwaVIx^*xP4(fk)pooK8GehE|_*7@T zyUpE6>8_?I4#}{GhY}9cm_#@{KAQK_O664|+uBi5Vbo+MWF}hRvhmLL)7(u)i{PQ z%?gegAQ-`Pk()uN*O*B-R->UWMrDM3wgpn|MrX@xXZ@5YW(n#z!fcJ>2`9uutaov? zJ7GiH$S^N>tqnNQ3*H>UTn=7lK*~65^CWCi83CLmdQT=qG@^w0@nOs?WYS56Q>-Nt zK4+zKt(x56Y)i1C3j}g1;R_l~gwrafnCLd(i;SCV^#}1K5k8&Jtnp<+i!T|H^>i6A z|1mc%a4(**h`rtu!9qcvK{!)`2xs|_#^*A=Vv&-6ce8UzjnlIvdukMC3-GIib2Q?F zMaI5td$Lq0G!(M!_?m@Dv6?ID?1uP81S@G64VSWQ^tv&zxgbTp?R5 z=tAx&uiyvjYvh_r57S>zW#7O6AXG8DL-4MYBG2m*#H5xxA{KQI% zk+aO@?uy`Ak=aPNPGb|{dX1kBGAAt6`m>k`kQiSXLYvukIoH}=@#*KW!8#-jzlfPp zW&+TcoEz9JX&vj9>G>_zR*6hMuX?!ee<>z!CETd7jj;Xm0?-q|4w2kR_?5;@gk2gp z4=}*Eb_*LIZrHM^xN)o1EQ810HP5-qy@eM$Blxwb-A(w7#%~F?nWnIne7nW?DlLFL zB7X;Aug32PcbX=v1ktewndI&b&aHW8T@z1{P3)*ExVs2NCi8_hS|}s2aCK@VJ1VAUvt@6ya$r*}oc}5t)Ay{-yCO;W>@xO=)cnSfF{Kt#Q}Ioo(khYur^4yeQz82rp~ALU`3!XOS}EyYZS2TJ5f%>D*~0#)VZ!ykJt}V>ez!)|aM|Dk_Tf9J5lpojY1Sd8eT5!YqxgBSYXLP{_^e@! zL5*dq)i|@RX6nm{*Dz{d=K5UqeFmxhWQ$epZ{XvZYBj;E4`5oEj@}?t2QtsjQrpO8 z#7Kv$iL&SU>(7HaNceoh!F(rm2ve;NHS5Ee4oaxhf||rTpQYBJbd^^%Ti@o=bI$HRBUxrog?M?2UYG-XCs;0;mA5zPAQguwVnrha^Ff~mq z)75l}n#L@Sh>yjy3BBroIv~ zAq_Lns;ZQbe88w#vcYPOGtk*gwL0FcPheWfn)fxM8kzT5Z$OFkiM$sp#fP$`DxYmo zbND%xNVJx+{BIHq)$ delta 387 zcmca9)y}$skFlPGfq|hYKhvmCC^f(U=jZhw2*89lLrI}>QfCHpL+Qo?bi85@uR65Ap5f4Z|3m8zZ*_EY5I($)?a2Ru=?<7 zdzA!Lp~DS_mnu5|jlcp9Z?G>^$zW|OR1GRrD=Ab@0tZQEUP)?Rab`(np@zEyGsmT^ kyhVB;(}GBCskNbfulv65`@Zk{rsv#uGnpj#=KIdO@80wO=iL9Ed*7Rb z(VN?Dyfx%(a(32qr1D8sv#z>k*VdZ)>Y8miRZxj+Jf}KR(!D9snNM+lKAq@J@pElm zeRXg3*el#umCP2qx=RhooJuux6*~)=gi?u8IIdF3TsB`&#bhB2C-4eqJAzS!TAe+p zO2(79j#M%qPGzML#-v!tah!Hcj!orzXQGhFHe!m#RKi*Pm`0dh+Zw8~WF3W4H)F!F zD)!ylj2W#?Yp5SH3A1Wj1EyV=9gRk}Jja(Yt!@Ex+6QjK++Lhr8=K6){2hsIoHH<> zac*onKZ^!=)}jmZdU0NDY#R5eT!B-`WpVxs%r>7)bS6|7^V2~#+nTW;$mRmV!jWt? z;KH%lI0anPo}Lt&Au*P_>+8HB@2Gbgi_$#j#e^jq)r5;|`p71`f55@gi7od3_OFAv{+S8a0{-D<`DgD9}~W*i44l=>!p0 zmov{aV@;6xRfJ~Ae4=D)C39Jpkqhf2@Or`qjg5pBW^1-C^Q;eCafuMF*0_dnt;Th+ z**s>Z8E7iWy@;X9tS%vJ63k}8^%^fGY%w({E1l{=s}(?|verB39a!r3clpO%+#r}X zLR4cb;l>FXG)b77qK3=rOBGzTQ~llEUhm{i|0x$Q5w4dKZq|4i;TEG?N~xTRnvtdB?cIUn2yagf5LN zA!i+DjdOI-Eq)*|XAc zja+eJ+;>K)4hsGKgg0wEKsZ#Ux0U`N>mu_HvaMh{Y}qPyVgyHo>@9?cG~P;h*i@c~ zT%P44XQGomar8`dj|$ymgtuu765d`hAv4E!jGvGr$I#;fdMDutjdu}-Mz?JV@3sKZ zc!zhlcbC7r6!(W*yhkv{2*)*^B%G+&YDR_7!1ofK)_5P` z{jn*`92C!~idbHsKUlrQ+c&?R3+Dp@d4}*ojb{lTGF2C-8vutbP*%^0t?-9kZ{J)O zrv>(5!bdbbfGno2h0nF4~Q>`&zA{b(fBIiYZ_ml zP()_;e#7QBrii{N;BOJWt??bgcPkddmi+fDN)}`1Jnz_M|Is_W`(1orfIlGoP~%60 zACHsByfuHq63OaY!GqKMrnSR&G?#zU#UI7zPlP{f{DtsWjlY>u>bU; z)oUX7hrs?x_?O1N3ICb!jxoyqxAHw}i}-&plsi5?{{hq%5 zx;9rWmJtrIgzu!P*=lvMXU9=neTN{8Rq!%@MU!Mw#a_gI&f} zs~Xd;Wjpq)FR4228Ae?^{NUv`&&ymPCoE|(GJujf0d6>PO?FzqYYR_vHy z=c+5YmzhRB&sZsQsYdB%o+g7`$yTdXroCE5MiC<)bZU+CGvX@a)y!6_wWht!wAZsO znlF5n%QUGC{Fsc#+lv|gi;|7U!}-okG8I-ExxZ~ksI93g+OJypft1U31lXt|{ItV( Xod~s+#OrF~bq(8MqM)wjdwaw4u$nmN literal 0 HcmV?d00001 diff --git a/xppl3716_init_pars/run0167_transform.npy b/xppl3716_init_pars/run0167_transform.npy new file mode 100644 index 0000000000000000000000000000000000000000..162774eb027d1d6467b671194c4faf1bc7e1e286 GIT binary patch literal 3413 zcmZ`+33yc16`lzMVIV9T7f>N8l7JYUK_y1W3|M0^&49mg?zE3%IQ)F0naLw@AYO8>PL&;(`lTtV*JMKVF7)w0|&Lu2k2eMNV!Am#G(vt50T4cH9gcTa+ z5mqKbtQ6W@Pk&j3aX#}dv&KQRitz=63pK)oh%vTar89Y~V&g?@EXnjPXJ-r-3FubF^3^GP;c$_YD~E{T;|2JhOky*s!eBy#3WlbejzFg#|Y~*))Ovg zwzbxz(c0i;VsZtcQ{zg)%QZG6=JK3LfR}A??y79WRHf{7W(Zw2OmTRN`;0s0KHs>- zJs8C+gsqzp*VstdR1;}L$gYkXE}LJnT$3EXMtELH*sSp?!mCYYGD;Ox(v;{B5*9Vv z#Nm!^are5r`l7g2fUhBJ(ReN4I$uf+u!7gwBuG)Yx4FmsKO_ZxJt3vhLr7aw@hoE2 z;cDp>lMJCxV=E!6ku&-7wyoMLH(2*1mCljH&TjXytDPgxL9XuC3qe0&KqF5mSVQRo zzhM-EVh4iKC=yBUP0roU*!a&K;Sk;^OY9`vtnntoE^9ik#x#mIi{);D*4RUMi^hn_mbaDF zY~5lRr1Ju?-iC;CpL=r$-YVc*341keBivq%#1u{h`?w}2m#XQbx3OA?d}qZS$4%6z z(BDaTyT&^RcU9?ajlYw#BGG%k`q;)So3CJFxLe5XA>6BRAK_i5?xcsRF~0kFbkYs? zACGRo&^C9@d3i)zD6@?sNe~UluYv~&VH_hL+%mx;SL-Y><0-?YCJ`FI#JISUiBQY zAlZlOw-q)!cewZ87{!MK@C@Nujpqo@o1*j5>%op%m^6c|rJbj^686oG;+P;mO!$a~ zLvW|=yL7`xdEccIa`{pwU(A&%_?XX&^$pH}*h1%)`WQYgT%RC(QsYyEPftA}hVU82 z$vG(nJ?I{tbV7Vqz@H<0UgHadFIr1k*V$)YCwxSFNo>AM_=?6?318FrdQBmj{{4pa zZ&D$BQ^4OMd|Tr?gzx$`#FqW{EJ`|dd!uvkN_YP@=gug;FTfuVeyH&y!jGq<$hsxDTy%M?3H{;rlt^7aG4L{K}Um+aSNTz?!}vi{du|{w?8m z8owv}!CFe1R+-B^9K|2S=1+t_Yy5@qSB<}!inWK*-z`Q;XntX-JL*2Wu|0-=2R2AnZ8c3M>b06~%<5Tc zHN%WgFyomlP3hQ!W4s*IES~TM74GEH+I7wFil_?73{%eLJE=J=wQ4ZqxhxHt{Vb@G z&Etu9ZCUYwqUQ6IUB_QpYJuo^l@s|+>Liw0HJb6sEGJ&c5p@bb*D?gusiJ3w(+p}M zORW}}@#!qhmz?=-sCsz@tNGlUs5x(Q)tRFAzLd-hoO}Jrt?vHLO;NR2rg)7dd?(e! zQmeDf_-vLF6*8pG;b%)C2S&_&!rALQ!l&Ee)iJeHCYk1316{^at7bE9VLAD%uc+m$ zo8+hWZK76)e!@>*yQ%YJiWgYPcT(rG)T-5tFJS5G7O!p9g{)<&@$WBo*jy?sdgh53 z>?)R8U1Y`=%gh8~{EJOpB6>!wHdbvcwYt=d+s*hgmSyvKuL`+{TEma&WU{B6n>f(D;2S1RCvK$A7Rg9l*FjfbG-Ab(18LRay%c+vOobNrY FF9HEmGgbfq literal 0 HcmV?d00001 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():