Commit 97acd415 authored by Gonzalez, Miguel's avatar Gonzalez, Miguel

Added to git and to runtime d17_2dnorm routines by Thomas

parent e7ff1f0d
;Modified 21.10.18 TS: Added option of background mask for mainbeam to d17_2dnorm_runs()
;Modified 20.10.18 TS: Removed foreground normalization of main beam and corrected a mistake in error calcs (d17_2dnorm_win, s17_2dnorm_runs(), d17_2dnrom_mask)
function d17_2dnorm,win,mainbeam,lrange=lrange
;+
; LAMP function to normalize a workspace by a main beam measurement
......@@ -15,7 +18,7 @@ function d17_2dnorm,win,mainbeam,lrange=lrange
; e.g. w3=d17_mbnorm(w2,w1,lrange=[2,20])
;
; CALLS TO:
; Last modified: 8.11.13 PG
; Last modified: 20.10.18 TS: Removed foreground normalization of main beam and corrected a mistake in error calcs
;-
; *** Collect the workspace information
......@@ -43,7 +46,8 @@ endelse
lmb=mbdatp.x
xout=windatp.x
yout=windatp.y
;print, n_elements(mainbeam)
;print, n_elements(lmb)
; *** Interpolate the wavelengths of the main beam on to the wavelengths for the data
mb=interpol(mainbeam,lmb,lwin)
emb=interpol(mbdatp.e,lmb,lwin)
......@@ -74,7 +78,7 @@ endfor
; *** Carry out the normalization
wout=win/mbdiv
eout=sqrt((windatp.e/mbdiv)^2+(ediv*wout/mbdiv)^2)
eout=sqrt((windatp.e/mbdiv)^2+(ediv*wout/mbdiv^2)^2)
; *** Check for NaN and infinity
wnan=where(finite(wout) ne 1)
......@@ -99,7 +103,7 @@ if n_elements(lrange) eq 2 then begin
endelse
endif else print,'lrange not set. Output is not truncated'
print, '2D_norm is done'
; *** Modifify the 'datp'
MOD_DATP,woutdatp,'x',xout
MOD_DATP,woutdatp,'y',yout
......@@ -108,3 +112,9 @@ GIVE_DATP,woutdatp
return,wout
end
This diff is collapsed.
This diff is collapsed.
;Modified 20.10.18 TS: Removed foreground normalization of main beam and corrected a mistake in error calcs
; ************************************************************************************
; *************** Function to normalize a 2D measurement with direct beam (main beam) *********************************
; ************************************************************************************
function d17_2dnorm_win,win,mainbeam, lrange=lrange,norm=norm,mask=mask
;+ NOTE: using the win version with multiple read runnumbers does not add the counttimes!! The normalization will be off!
; LAMP function to normalize a workspace by a main beam measurement using two workspaces loaded in lamp
; Hacked from d17_2Dnorm
; The workspaces do not need to be treated, just loaded, for example with rdrun('')
; NECESSARY VARIABLES:
; win=the input workspace
; mainbeam=the workspace with the main beam data
;
; OPTIONAL VARIABLES:
; lrange = a two-element array with the maximum and minimum values of lambda
; norm = a string to define the normalization of the data (default is time)
; mask = define the foreground for the main beam measurement. Only box allowed at the moment without restricting the wavelength. See "maskme" for further details.
;
;
; e.g. w3=d17_mbnorm(w2,w1,lrange=[2,20])
;
; CALLS TO:
; Last modified: 13.12.16 TS
;-
d17_errmsg,1,'//'
d17_errmsg,1,'2D normalization of given workspace'
d17_errmsg,1,'//'
; *** Collect the workspace information
TAKE_DATP,windatp
;PG: error calc:
if (total(windatp.e) eq 0) then MOD_DATP,windatp, 'E' , sqrt(win)
ezero=where(windatp.e eq 0)
if (ezero[0] ne -1) then windatp.e[ezero]=1
TAKE_DATP,mbdatp,/THIRD
;PG: error calc:
if (total(mbdatp.e) eq 0) then MOD_DATP,mbdatp, 'E' , sqrt(mainbeam)
ezero=where(mbdatp.e eq 0)
if (ezero[0] ne -1) then mbdatp.e[ezero]=1
;TS Figure out if normalization string is given
if n_elements(norm) ne 1 then norm='time'
if (norm ne 'time') then begin
if (norm ne 'mon') then begin
if (norm ne 'none') then begin
d17_errmsg,1,'ERROR: You must enter time, mon or none for normalization'
return,0
endif else begin
print, 'no normalization'
endelse
endif else begin
d17_errmsg,1,'ERROR: Currently the MON is not read'
return,0
endelse
endif else begin
print, 'Time Normalization'
endelse
;Normalize win by time or
if (norm ne 'none') then begin
time_win=windatp.p[1]
endif else begin
time_win=1
endelse
;Normalize MB by time or
if (norm ne 'none') then begin
time_mb=mbdatp.p[1]
endif else begin
time_mb=1
endelse
;;mon_all=[mon_all,tdatp.p[2]]
;
;
print, 'Normalized WIN by time:',string(time_win)
print, 'Normalized MB by time:',string(time_mb)
d17_errmsg,1,'CountTimes: data ' + strtrim(string(time_win),1) +'; mainbeam '+ strtrim(string(time_mb),1)
;;TS Determine if the masks are correctly defined and find out the limits::
if n_elements(mask) ne 1 then begin
d17_errmsg,1,'No mainbeam mask is applied!!!!!!'
d17_errmsg,1,'Setting ForegroundWidth to 1!!!!!!'
mask='none'
foregroundwidth=1
endif else begin
bra=strpos(mask,'[')
ket=strpos(mask,']')
limstr=strmid(mask,bra+1,ket-bra)
flims=0
comma=0
while comma ne -1 do begin
comma=strpos(limstr,',')
if comma ne -1 then begin
flims=[flims,fix(strmid(limstr,0,comma))]
endif else flims=[flims,fix(strmid(limstr,0,strlen(limstr)-1))]
limstr=strmid(limstr,comma+1,strlen(limstr)-comma+1)
endwhile
flims=flims[1:n_elements(flims)-1]
foregroundwidth=flims[1]-flims[0]
print,'THESE ARE THE foreground LIMITS ', string(flims)
endelse
;if n_elements(background) ge 1 then begin
; d17_errmsg,1,'Sorry, Backround subtraction is not implemented in this version'
;
; print, 'Background subtraction not implemented in this version'
;endif
woutdatp=windatp
winsize=size(win,/DIMENSIONS)
;Check if win was given in 1d or 2d
if n_elements(winsize) eq 1 then begin
cflag=0
lwin=windatp.x
endif else begin
cflag=1
lwin=windatp.y
endelse
;Check if mainbeam was given in 1d or 2d
winsizemb=size(mainbeam,/DIMENSIONS)
if n_elements(winsizemb) eq 1 then begin
cflagmb=0
lmb=mbdatp.x
endif else begin
cflagmb=1
lmb=mbdatp.y
endelse
; *** Mask the MainBeam workspace if required
if cflagmb eq 1 then begin
mainbeam=maskme(mask,mainbeam,datp=mbdatp)
;*** Check dimensions of the workspace. If it is 1D, continue else integrate horizontally
;;** NOTE: This makes the Y-axis the new X-axis
w_size=size(mainbeam)
if w_size[0] eq 2 then begin
mainbeam=total(mainbeam,1)
endif
;recalculate the errors for the ROI
emainb=sqrt(mainbeam)
;emainb[where (mainbeam eq 0.0)] =1.0 ;Not done at this stage
endif else begin
emainb=sqrt(mainbeam)
;emainb[where (mainbeam eq 0.0)] =1.0
endelse
xout=windatp.x
yout=windatp.y
;Normalize data and mainbeam by time and foregroundwidth
win=win/time_win;
ein=windatp.e/time_win;
mainbeam=mainbeam/time_mb;/foregroundwidth ;TS: Shouldn't be normed by foreground
emainb=emainb/time_mb;/foregroundwidth ;TS: Shouldn't be normed by foreground
;print, 'Foregroundwidth:', foregroundwidth
;print, 'mainbeam',n_elements(mainbeam)
;print, 'win',n_elements(lwin)
;print, 'MB',n_elements(lmb)
; *** Interpolate the wavelengths of the main beam on to the wavelengths for the data
mb=interpol(mainbeam,lmb,lwin)
emb=interpol(emainb,lmb,lwin)
;mb=mainbeam
;emb=mbdatp.e
; *** Create a workspace with main beam, the same dimensions as the data
if cflag eq 1 then begin
dmb=reform(mb,1,n_elements(mb),/OVERWRITE)
emb=reform(emb,1,n_elements(emb),/OVERWRITE)
endif
winsize=[winsize,0,0]
winsize=winsize[0:2]
tdmb=dmb
temb=emb
for i=1,winsize[0]-1 do begin
dmb=[dmb,tdmb]
emb=[emb,temb]
endfor
mbdiv=dmb
ediv=emb
for i=1,winsize[2]-1 do begin
mbdiv=[[[mbdiv]],[[dmb]]]
ediv=[[[ediv]],[[emb]]]
endfor
; *** Carry out the normalizations
wout=win/mbdiv
eout=sqrt((windatp.e/mbdiv)^2+(ediv*wout/mbdiv^2)^2)
; *** Check for NaN and infinity
wnan=where(finite(wout) ne 1)
if wnan[0] ne -1 then begin
wout[wnan] = 0.0
eout[wnan] = 1.0
endif
; *** Determine whether a lambda range has been set
if n_elements(lrange) eq 2 then begin
lout=where(lwin ge lrange[0] and lwin le lrange[1])
if cflag eq 1 then begin
yout=yout[[lout]]
wout=wout[*,[lout],*]
eout=eout[*,[lout],*]
endif else begin
xout=xout[[lout]]
wout=wout[[lout]]
eout=eout[[lout]]
endelse
endif else d17_errmsg,1,'lrange not set. Output is not truncated'
; *** Modifify the 'datp'
MOD_DATP,woutdatp,'x',xout
MOD_DATP,woutdatp,'y',yout
MOD_DATP,woutdatp,'e',eout
GIVE_DATP,woutdatp
return,wout
end
\ No newline at end of file
......@@ -55,6 +55,9 @@ print,'**********'
;.run d17tof2.pro ;commented by CS (procedure not present)
;.run d17tof3.pro ;commented by CS (procedure not present)
.run d17_2dnorm.pro
.run d17_2dnorm_mask.pro
.run d17_2dnorm_runs.pro
.run d17_2dnorm_win.pro
.run slitwit.pro
print,'**********'
......
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