Skip to content

GitLab

Commit 477cc82a authored by Gonzalez, Miguel's avatar Gonzalez, Miguel
Browse files

Corrected error in formula for QENS_IL2 in interface. Added extra output for... 

Corrected error in formula for QENS_IL2 in interface. Added extra output for Salsa, taking care of the case where the xtra information is not available (Issue 25)
parent 4237df39
Hide whitespace changes
Inline Side-by-side
lamp_mac/str_fit.pro
View file @ 477cc82a
......@@ -1743,7 +1743,7 @@ case uv[0] of
Pexplain[6] = 'HWHM (meV) of localized component.'
end
;QENS_IL2
85: begin & formu ='Y = QENS_IL1 (X, P1, P2, P3, P4, P5, P6, P7, P8) ;QENS Ionic liquids (+ two localized components)'
85: begin & formu ='Y = QENS_IL2 (X, P1, P2, P3, P4, P5, P6, P7, P8) ;QENS Ionic liquids (+ two localized components)'
npp = 8 ;Number of parameters
editf = 0 ;Not editable formula
fixpi = [0,0,0,0,0,0,0,0] ;Fixed params (definitively)
......@@ -7863,7 +7863,7 @@ for ispec=prems,derch,fol do begin
DD2k=di_conv(DD2k,wres) ;*Wf to keep a good position
endif else DD2k=(DD2k-Wf)
DD2k[idk]=0. ;only what was higher Wf
;stop
BB =di_keeb(DD2K>0.,0,(DD2p-Wf)/2.,idd0,shf) & BB[0:wlim-1]=0 & BB[n-wlim:n-1]=0 ;get Bands in BB
idx =where(BB[idd0] gt 0.) ; ""
if idx[0] ge 0 then APx=idd0[idx] else begin & bid=max(DD2G,APx) & BB[APx]=DD2G[APx] & endelse
......@@ -7948,7 +7948,7 @@ for ispec=prems,derch,fol do begin
AH2 =atmp[ind] & AP2=atmp[ind+1] & WD2=abs(atmp[ind+2]) & WDs=s[ind+2] & st=sort(AP2)
APx =round((AP2-Xa)/Xunit)
AH2 =AH2[ st] & AP2=AP2[ st] & WD2=WD2[st] & APx=APx[st] & WDs=WDs[st]
;stop
idk =where((AH2 gt 0) and (AP2 ge Xa) and (AP2 le Xb));height>0 and position in margin
if idk[0] ge 0 then begin
AH2 =AH2[idk] & AP2=AP2[idk] & WD2=WD2[idk] & APx=APx[idk] & WDs=WDs[idk] & nek=n_elements(AP2)
......@@ -8578,7 +8578,7 @@ end4:
;**************************************************
;**************************************************
if (modl eq 1) or (modl eq 2) then if not sep2 then begin ;return a one peak parameters set from a double peak
if (modl eq 1) or (modl eq 2) then if not sep2 then begin ;return a one peak parameter set from a double peak
for i=0,nspec-1 do begin & fct=fit_pars[i,1+0,*] * fit_pars[i,1+2,*] & tfct=total(fct)
for j=0,npr do begin & fit_pars[ i,j,0]=total(fit_pars[ i,j,*]*fct,3)/tfct
fit_sigmas[i,j,0]=total(fit_sigmas[i,j,*]*fct,3)/tfct & endfor
......@@ -8606,97 +8606,154 @@ endif
if swww ne 5 then begin & fit_pars=float(fit_pars) & fit_sigmas=float(fit_sigmas) & endif
;Where the data come from (D20 ID31 SALSA Unknown)
ptxt=datp.par_txt & xtr =[0.] & ytr =[0.] & ztr=[0.] & omg=[0.] & tth=[0.] & xro=[0.] & yro=[0.] & zro=[0.] & inst='?'
moni=reform([datp.N]) & time=reform([datp.N])
ptxt=datp.par_txt & xtr =[0.] & ytr =[0.] & ztr=[0.] & omg=[0.] & tth=[0.]
xro=[0.] & yro=[0.] & zro=[0.] & inst='?'
moni=reform([datp.N]) & time=reform([datp.N])
;D20
if n_elements(ptxt) gt 18 then if strpos(ptxt[17],'Tr.1') eq 0 then begin & idx=17 & ido=14 & idt=13 & inst='D20'
if (size(datp.pv))[1] eq n_elements(datp.p) then begin & xtr=datp.pv[idx,*] & ytr=datp.pv[idx+1,*];& ztr=datp.pv(idx+2,*)
omg=datp.pv[ido,*] & tth=datp.pv[idt,*]
endif else begin & xtr=datp.p[idx] & ytr=datp.p[idx+1] & omg=datp.p[ ido] & tth=datp.p[ idt] & endelse & endif
if n_elements(ptxt) gt 18 then if strpos(ptxt[17],'Tr.1') eq 0 then begin
idx=17 & ido=14 & idt=13 & inst='D20'
if (size(datp.pv))[1] eq n_elements(datp.p) then begin
xtr=datp.pv[idx,*] & ytr=datp.pv[idx+1,*] ;& ztr=datp.pv(idx+2,*)
omg=datp.pv[ido,*] & tth=datp.pv[idt,*]
endif else begin
xtr=datp.p[idx] & ytr=datp.p[idx+1]
omg=datp.p[ ido] & tth=datp.p[ idt]
endelse
endif
;ID31
if strpos(ptxt[0],'#F ') eq 0 then begin & idx=0 & ido=0 & idt=0
if strpos(ptxt[0],'#F ') eq 0 then begin
idx=0 & ido=0 & idt=0
for i=0,n_elements(ptxt)-1 do if strpos(ptxt[i],'X Translation') eq 0 then idx=i else $
if strpos(ptxt[i],'Omega=') eq 0 then ido=i else if strpos(ptxt[i],'2_theta=') eq 0 then idt=i else $
if strpos(ptxt[i],'Chi=') eq 0 then idc=i & if idx gt 0 then inst ='ID31'
if idx gt 0 then if (size(datp.pv))[1] eq n_elements(datp.p) then begin & xtr =datp.pv[idx,*] & ytr =datp.pv[idx+1,*] & ztr=datp.pv[idx+2,*]
omg =datp.pv[ido,*] & tth =datp.pv[idt,*]
endif else begin & xtr=datp.p[idx] & ytr=datp.p[idx+1] & ztr =datp.p[idx+2] & omg =datp.p[ ido] & tth=datp.p[ idt] & endelse & endif
if strpos(ptxt[i],'Omega=') eq 0 then ido=i else $
if strpos(ptxt[i],'2_theta=') eq 0 then idt=i else $
if strpos(ptxt[i],'Chi=') eq 0 then idc=i & if idx gt 0 then inst ='ID31'
if idx gt 0 then if (size(datp.pv))[1] eq n_elements(datp.p) then begin
xtr =datp.pv[idx,*] & ytr =datp.pv[idx+1,*] & ztr=datp.pv[idx+2,*]
omg =datp.pv[ido,*] & tth =datp.pv[idt,*]
endif else begin
xtr=datp.p[idx] & ytr=datp.p[idx+1] & ztr =datp.p[idx+2]
omg =datp.p[ ido] & tth=datp.p[ idt]
endelse
endif
;SALSA
if strpos(ptxt[0],'2*Th') eq 0 then begin & idx=0 & ido=1 & idt=0
for i=0,n_elements(ptxt)-1 do if strpos(ptxt[i],'X_translation') eq 0 then idx=i & if idx gt 0 then inst='SALSA'
if idx gt 0 then if (size(datp.pv))[1] eq n_elements(datp.p) then begin & xtr =datp.pv[idx,*] & ytr =datp.pv[idx+1,*] & ztr=datp.pv[idx+2,*]
xro =datp.pv[idx+3,*] & yro =datp.pv[idx+4,*] & zro=datp.pv[idx+5,*]
omg =datp.pv[ido,*] & tth =datp.pv[idt,*]
endif else if (((size(datp.pv))[2] eq vspec) and ((size(datp.pv))[1] ge 10)) $
or ( (vspec eq 1) and (n_elements(datp.pv) ge 10) and (n_elements(datp.pv) lt n_elements(datp.p)) ) $
then begin & idx =5 & ido=4 & idt=3 & datp.pv=round(datp.pv*10.)/10000.
xtr =datp.pv[idx,*] & ytr =datp.pv[idx+1,*] & ztr=datp.pv[idx+2,*]
xro =datp.pv[idx+3,*] & yro =datp.pv[idx+4,*] & zro=datp.pv[idx+5,*]
omg =datp.pv[ido,*] & tth =datp.pv[idt,*]
time=datp.pv[0,*]*1000 & moni=datp.pv[1,*]*1000
endif else begin & xtr=datp.p[idx] & ytr=datp.p[idx+1] & ztr =datp.p[ idx+2] & omg =datp.p[ ido] & tth=datp.p[ idt]
xro =datp.p[ idx+3] & yro =datp.p[ idx+4] & zro=datp.p[ idx+5] & endelse & endif
nxtr =n_elements(xtr) -1 & nytr =n_elements(ytr) -1 & nztr=n_elements(ztr)-1
nomg =n_elements(omg) -1 & ntth =n_elements(tth) -1 & nxro=n_elements(xro)-1 & nyro=n_elements(yro)-1 & nzro=n_elements(zro)-1
nmoni=n_elements(moni)-1 & ntime=n_elements(time)-1
if strpos(ptxt[0],'2*Th') eq 0 then begin
idx=0 & ido=1 & idt=0
for i=0,n_elements(ptxt)-1 do begin
if strpos(ptxt[i],'X_translation') eq 0 then idx=i
if idx gt 0 then inst='SALSA'
endfor
if idx gt 0 then if (size(datp.pv))[1] eq n_elements(datp.p) then begin
xtr =datp.pv[idx,*] & ytr =datp.pv[idx+1,*] & ztr=datp.pv[idx+2,*]
xro =datp.pv[idx+3,*] & yro =datp.pv[idx+4,*] & zro=datp.pv[idx+5,*]
omg =datp.pv[ido,*] & tth =datp.pv[idt,*]
endif else begin
if (((size(datp.pv))[2] eq vspec) and ((size(datp.pv))[1] ge 10)) $
or ( (vspec eq 1) and (n_elements(datp.pv) ge 10) and (n_elements(datp.pv) lt n_elements(datp.p)) ) then begin
idx =5 & ido=4 & idt=3
datp.pv=round(datp.pv*10.)/10000.
xtr = datp.pv[idx,*] & ytr = datp.pv[idx+1,*] & ztr = datp.pv[idx+2,*]
xro = datp.pv[idx+3,*] & yro = datp.pv[idx+4,*] & zro = datp.pv[idx+5,*]
omg = datp.pv[ido,*] & tth = datp.pv[idt,*]
time = datp.pv[0,*]*1000 & moni = datp.pv[1,*]*1000
chi = datp.pv[20,*] & phi = datp.pv[21,*]
angl = datp.pv[12,*]
endif else begin
xtr = datp.p[idx] & ytr = datp.p[idx+1] & ztr = datp.p[ idx+2]
xro = datp.p[ idx+3] & yro = datp.p[ idx+4] & zro = datp.p[ idx+5]
omg = datp.p[ ido] & tth = datp.p[ idt]
endelse
endelse
endif
nxtr =n_elements(xtr) -1 & nytr =n_elements(ytr) -1 & nztr=n_elements(ztr)-1
nomg =n_elements(omg) -1 & ntth =n_elements(tth) -1 & nxro=n_elements(xro)-1 & nyro=n_elements(yro)-1 & nzro=n_elements(zro)-1
nmoni=n_elements(moni)-1 & ntime=n_elements(time)-1
;WRITE OUTPUT FILE
;*****************
title=title+' Fit_function:'
if modl eq 5 then for i =0,nk-1 do title=title+' ,'+UserModl.name[i] else $
if modl eq 4 then title =title+'FREE' else $
if frac eq 0. then title =title+'Gauss' else $
if frac eq 1. then title =title+'Lorentz' else title=title+'Pseudo_Voigt'
if sb gt 0 then bktxt ='(+a1.shapeBgd) '
if sb le 0 then bktxt ='(+a1.x +a2.x^2) '
ChiText = 'Mean CHI square = '+ string(round(ChiMean*100.)/100. ,format='(F6.2)')
title=title+' Fit_function:'
if modl eq 5 then for i =0,nk-1 do title=title+' ,'+UserModl.name[i] else $
if modl eq 4 then title =title+'FREE' else $
if frac eq 0. then title =title+'Gauss' else $
if frac eq 1. then title =title+'Lorentz' else title=title+'Pseudo_Voigt'
if sb gt 0 then bktxt ='(+a1.shapeBgd) '
if sb le 0 then bktxt ='(+a1.x +a2.x^2) '
if n_elements(datp.y) eq sw[2] then YY=datp.y else YY=indgen(sw[2])+1
ChiText = 'Mean CHI square = '+ string(round(ChiMean*100.)/100. ,format='(F6.2)')
if n_elements(datp.y) eq sw[2] then YY=datp.y else YY=indgen(sw[2])+1
if n_elements(file_out) ne 1 then file_out=''
if file_out gt ' ' then begin
on_ioerror,miswrt & u=-1
explain=['']
if keyword_set(append) then openw,u,file_out,/append,/get_lun else openw,u,file_out,/get_lun
if keyword_set(append) then begin & printf,u,'' & printf,u,'***New fit' & printf,u,'***New fit' & printf,u,''
endif else if inst ne 'SALSA' then explain=['',$
'Sigma=s return errors for all fitted parameters [nspectra , nparams , nfunct]',$
'Res_dfit=d return filtered fit (no background, deconvoluted) [nx , nspectra]' ,$
'Fit=f return fitted spectra (+background, convoluted) [nx , nspectra]' ,$
'Fout=n return fitted spectra in Wrks number N [nx , nspectra]' ,$
'Convfunc=c return convoluted functions (no background) [nx , nspectra , nfunct]' ,$
'out_subF=r return deconvoluted functions indexed by out_subi=[index] [nx , nspectra , nindex]']
on_ioerror,miswrt & u=-1
explain=['']
if keyword_set(append) then openw,u,file_out,/append,/get_lun else openw,u,file_out,/get_lun
if keyword_set(append) then begin & printf,u,'' & printf,u,'***New fit' & printf,u,'***New fit' & printf,u,''
endif else if inst ne 'SALSA' then explain=['',$
'Sigma=s return errors for all fitted parameters [nspectra , nparams , nfunct]',$
'Res_dfit=d return filtered fit (no background, deconvoluted) [nx , nspectra]' ,$
'Fit=f return fitted spectra (+background, convoluted) [nx , nspectra]' ,$
'Fout=n return fitted spectra in Wrks number N [nx , nspectra]' ,$
'Convfunc=c return convoluted functions (no background) [nx , nspectra , nfunct]' ,$
'out_subF=r return deconvoluted functions indexed by out_subi=[index] [nx , nspectra , nindex]']
if swww eq 5 then begin & fit_pars=float(fit_pars) & fit_sigmas=float(fit_sigmas) & endif
if swww eq 5 then begin & fit_pars=float(fit_pars) & fit_sigmas=float(fit_sigmas) & endif
if oklab then widget_control,ilab,bad_id=ii,set_value=strmid('-> '+file_out+' ',0,20)
printf,u,title
if datp.other_tit gt '' then printf,u,datp.other_tit
if datp.history gt '' then printf,u,datp.history
printf,u,ChiText
if n_elements(incmd) eq 1 then printf,u,'Call: '+incmd else printf,u,' '
printf,u,' '
if inst ne '?' then begin
printf,u,'Spectra pseudo_V e height e center e fWmH e background e ' +bktxt+ ' a2 e fit_area e '+$
' raw_int Xtransl. Ytransl. Ztransl. Omega 2_theta'+$
' Xrotation Yrotation Zrotation Monitor Time'
for i=0,nspec-1 do for k=0,nk-1 do $
printf,u,YY[frst+i],' ',(fit_pars[i,0,k]>0),(fit_sigmas[i,0,k]),$
(fit_pars[i,1,k]) ,(fit_sigmas[i,1,k]),$
(fit_pars[i,2,k]) ,(fit_sigmas[i,2,k]),$
(fit_pars[i,3,k]) ,(fit_sigmas[i,3,k]),$
(fit_pars[i,4,k]) ,(fit_sigmas[i,4,k]),$
(fit_pars[i,5,k]) ,(fit_sigmas[i,5,k]),$
(fit_pars[i,6,k]) ,(fit_sigmas[i,6,k]),$
(fit_pars[i,9,k]) ,(fit_sigmas[i,9,k]),$
(fit_pars[i,7,k]),$
xtr[(frst+i)<nxtr],ytr[(frst+i)<nytr],ztr[(frst+i)<nztr],omg[(frst+i)<nomg],$
tth[(frst+i)<ntth],xro[(frst+i)<nxro],yro[(frst+i)<nyro],zro[(frst+i)<nzro],$
long(moni[(frst+i)<nmoni]),long(time[(frst+i)<ntime]),$
format='(F7.2,A,F5.3,4X,F5.3,23G,4X,I8,1X,I8)'
if oklab then widget_control,ilab,bad_id=ii,set_value=strmid('-> '+file_out+' ',0,20)
printf,u,title
if datp.other_tit gt '' then printf,u,datp.other_tit
if datp.history gt '' then printf,u,datp.history
printf,u,ChiText
if n_elements(incmd) eq 1 then printf,u,'Call: '+incmd else printf,u,' '
printf,u,' '
chi = 1.1111
phi = 2.2222
angl = 3.3333
if inst ne '?' then begin
printf, u, format='("Spectra", 2x, "pseudo_V", 1x, "e", 11x, "height", 9x, "e", $)'
printf, u, format='(14x, "center", 7x, "e", 16x, "fWmH", 11x, "e", 14x, "background", $)'
printf, u, format='(6x, "e", 14x, A, 2x, "a2", 13x, "e", 13x, "fit_area", 7x, "e", $)', bktxt
printf, u, format='(14x, "raw_int", 9x, "Xtransl.", 7x, "Ytransl.", 7x, "Ztransl.", $)
printf, u, format='(7x, "Omega", 10x, "2_theta", 6x, "Xrotation", 4x, "Yrotation", $)
printf, u, format='(7x, "Zrotation", 5x, "Monitor", 3x, "Time", $)
;extra fields added for Salsa
if n_elements(chi) gt 0 then printf, u, format='(14x, "Chi", $)'
if n_elements(phi) gt 0 then printf, u, format='(12x, "Phi", $)'
if n_elements(angl) gt 0 then printf, u, format='(10x, "Anglesx1000", $)'
printf, u, format=''
for i=0,nspec-1 do begin
for k=0,nk-1 do begin
printf, u, YY[frst+i], ' ', $
(fit_pars[i,0,k]>0), (fit_sigmas[i,0,k]), $
(fit_pars[i,1,k]) , (fit_sigmas[i,1,k]), $
(fit_pars[i,2,k]) , (fit_sigmas[i,2,k]), $
(fit_pars[i,3,k]) , (fit_sigmas[i,3,k]), $
(fit_pars[i,4,k]) , (fit_sigmas[i,4,k]), $
(fit_pars[i,5,k]) , (fit_sigmas[i,5,k]), $
(fit_pars[i,6,k]) , (fit_sigmas[i,6,k]), $
(fit_pars[i,9,k]) , (fit_sigmas[i,9,k]), $
(fit_pars[i,7,k]), $
xtr[(frst+i)<nxtr], ytr[(frst+i)<nytr], ztr[(frst+i)<nztr], $
omg[(frst+i)<nomg], tth[(frst+i)<ntth], $
xro[(frst+i)<nxro], yro[(frst+i)<nyro], zro[(frst+i)<nzro],$
long(moni[(frst+i)<nmoni]), $
long(time[(frst+i)<ntime]), $
format='(F7.2,A,F5.3,4X,F5.3,23G,4X,I8,1X,I8,$)'
;print extra info (for Salsa) only if available
if n_elements(chi) gt 0 then printf, u, chi[(frst+i)<nxtr], format='(1X, G, $)'
if n_elements(phi) gt 0 then printf, u, phi[(frst+i)<nxtr], format='(1X, G, $)'
if n_elements(angl) gt 0 then printf, u, angl[(frst+i)<nxtr], format='(1X, G, $)'
printf, u, format=''
endfor
endfor
endif else begin
printf,u,'Spectra pseudo_V e P1.height e P2.center e P3.fWmH e background e ' +bktxt+ ' a2 e fit_area e '+$
' raw_int'
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment