in5_t2e.pro 9.87 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
function in5_t2e,w_in, w_van, elp=elp, average_elp=average_elp, verbose=verbose
;========================================================================================================
;+
;
;<h2>NAME:</h2>
;
;	t2e 
;
;<h2>PURPOSES:</h2>
; Converts array w_in  to energy in w_out
; Based on the original "tee" procedure instead of on "t2e". 
;
;<h2>COMMAND SYNTAX:</h2>
;	w2 = t2e(w1)
;         convert tof spectrum in w1 to energy and output in w2
;         Energies are in x2.
;   w2 = t2e(w1,/verbose)    ; same with information on the standard output.
;   w2 = t2e(w1,w3,/verbose) ; Uses an elastic peak standard raw data reference for finding the 
;                              elastic peak position (ELP). 
;                              The Vanadium or better the low temperature sample is this 
;                              reference usually. This is required each time the ELP is incorrectly
;                              defined in the sample. E.g. on liquid spectra.
;
;<h2>ARGUMENTS:</h2>
; (all optional)
;
; <b>/average_elp</b>     : take average elp instead of 
; <b>W_van</b>            : Workspace for the vanadium/low temp. sample, that gives the elastic peak reference 
; <b>/verbose</b>         : more info during execution 
; 
;<h2>DIMENSIONS:</h2>
;
; w_in = w_out(nchannels, nspectra)
;
;<h2>VERSION HISTORY:</h2>
;
;	Written by  GJK 1994
; 	Revisions:     JO 2004-11-03 Include other Tof instruments (e.g.: MIBEMOL, LLB)
;		           JO 2006-03-28 Html help, more verbosity
;                  JO 2006-10-23 tee -> t2e with the same properties 
;                                in case the original t2e fails.
;                  JO 2006-11-10 gaussian fit with str_fit instead of fitgauss (due to 
;                                problems when error bars are not computed), with the same restriction
;                                to the elastic peak than before +/- 30 channels from  
;                                mean elastic maximum.
;                  JO 2006-12-14 Reference of the elastic line by the vanadium (keyword: Wvan)
;                                Vanadium data must be in the same format than sample data.
;
;                  JO, 2007-11-05 Seems that the gauss ELP is 1 channel out from the reality.
;                                 Is it always true ???
;                  JO, Sat Oct 11 11:57:30 CEST 2008: If spectra are zeros remove first and last index
;                                                    str_fit gives the first X if no maxi found
;                  JO, Thu Dec  4 10:55:17 CET 2008: correction on the index ind
;                  JO, Fri Nov 19 12:11:00 CET 2010: Change the intensities correction. speed-up
;                                                    S(Q,w) transformation by matricing. 
;                                                    Warning: Intensities corrections with CHW in [micro-sec] 
;                                                    to be concistent with vnorm.          
;                  JO 2011-09-14 Correction error on the ELP. Add user defined ELP when not enough
;                                ELP and no vanadium.  ELP might be average after distance correction.
;
;-
;
;========================================================================================================
;;	@lamp.cbk

; for functions that has to appear later but are unknown from here when saving
; using: save,/routines,'t2e',filename='t2e.sav'
    forward_function str_fit

    COMMON c_lamp_access, inst

    take_datp, datp

	x = 0.0
	y = 0.0
	w = 0.0
	n_in = datp.n
	pv_in= datp.pv     
	p_in = datp.p       
	x_in = datp.x
	y_in = datp.y
	z_in = datp.z
	e_in = datp.e
        if n_elements(e_in) ne n_elements(w_in) then $
            e_in = sqrt(w_in)

	; set defaults            
	x_out = x_in
	e_out = e_in
	
	sizeWin = size(w_in)
	if sizeWin[0] lt 1 then  begin
		mess='T2E: ERROR: Workspace empty!'
		widget_control,bad_id=iii,l_message,set_value=mess
        print, 'T2E: ERROR!: Workspace empty'
		return,w_in
	endif
		
	nchannels = sizeWin[1]
	if sizeWin[0] gt 1 then nspectra = sizeWin[2] else nspectra = 1

	coef  = fltarr(nchannels)
	x_out = fltarr(nchannels,nspectra)
	w_out = w_in*1.0
	if n_elements(e_in)  ne n_elements(w_in) then e_out=0
	
; --------------------------------------------------------------------
;       Convert to Energy
; --------------------------------------------------------------------
    Ei = 81.8049/p_in[21]^2                     ; Ei in meV
    telast = p_in[27]*p_in[21]/3956.035         ; TOF in sec
    cwidth = p_in[18]  *1e-6                    ; chw from mu-sec to sec

	if keyword_set(verbose) then begin
		print, 'T2E: Instrument is                 : ',strupcase(inst)
		print, 'T2E: Sample-detector distance      : ',strtrim(string(p_in[27],format='(F8.2)'),2),' m'
    endif        

; --------------------------------------------------------------------
;       Find the elastic peak position for each group     
;       ELP coming possibily from the vanadium (or low temp) run  
; --------------------------------------------------------------------
    w_elp = w_out
    if  (size(w_van))[0] ne 0 then begin
        if ( (size(w_van))[1] ne (size(w_out))[1] ) or $
            ( (size(w_van))[2] ne (size(w_out))[2] ) then $
                print,'T2E: When calling vanadium, workspaces must have same format!' $
    else begin
        w_elp = w_van    
        if keyword_set(verbose) then $
            print,'T2E: Transformation according to VANADIUM elastic peak' $
        else $
            print,'T2E: Transformation according to SAMPLE elastic peak'   
        endelse    
    endif
    if (size(w_elp))[0] eq 2 then $
        ytot = total(w_elp,2) $
    else $
        ytot = w_elp
    y0    = max(ytot,i0)
    xmin  = x_in[(i0 - 30) > 0]
    xmax  = x_in[(i0 + 30) < (nchannels - 1)]
    ind   = where(x_in gt xmin and x_in lt xmax)  
    epp   = fltarr(nspectra)
    gauss = epp
    gauss = str_fit(w_elp[ind,*],Bkgd=0,Xin=x_in[ind])
; ----------------------------------------------------------------------------
;        epp   = gauss(*,2) 
; JO, 2007-11-05 Seems that the gauss ELP is 1 channel out from the reality.
;                Is it always true ???
; ----------------------------------------------------------------------------------       
    epp   = gauss[*,2] - 1.0
153
  
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175

; ----------------------------------------------------------------------------------       
;       Transform into S(2theta,W) 
; ----------------------------------------------------------------------------------               
; ----------------------------------------------------------------------------
;  JO, Sat Oct 11 11:57:30 CEST 2008
;  If spectra are zeros remove first and last index
;  str_fit gives the first X if no maxi found
; ----------------------------------------------------------------------------
;  before:   ind   = where(x_in+1. gt xmin and x_in lt xmax-1.)  
; ----------------------------------------------------------------------------
;  JO, Thu Dec  4 10:55:17 CET 2008: correction on the index ind
;   before:   ind   = where(x_in+1. gt xmin and x_in lt xmax-1.)  
; ----------------------------------------------------------------------------
    if keyword_set(average_elp) then begin
        ind         = where(epp gt xmin and epp lt xmax)  
        average_elp = mean(epp[ind])
        epp = average_elp # (fltarr(nspectra)+1.)
       if keyword_set(verbose) then  $
            print, 'T2E: Elastic peak position forced to (AVERAGE_ELP): ', strtrim(string( average_elp ),2)   
    endif 
    if keyword_set(elp) then begin
176
177
;      same on the ELP index: should be read indice - 1 (see above)
        average_elp = elp - 1
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
        epp = average_elp # (fltarr(nspectra)+1.)
        if keyword_set(verbose) then $
           print, 'T2E: Elastic peak position forced to user value (ELP) : ', strtrim(string( elp ),2)
    endif
    
   if keyword_set(verbose) then begin
        print, 'T2E: Average elastic peak position : ',strtrim(string( round( mean(epp[ind]) ) ),2)
        print, 'T2E: Incident energy               : ',strtrim( string(Ei) ,2),' meV'
        print, 'T2E: Channel width                 : ',strtrim(string(cwidth*1e6),2),' mu-sec'
   endif

    x_dummy = w_out*0.0
    j   = indgen(nchannels) # (fltarr(nspectra)+1.)
    epp = (fltarr(nchannels)+1.) # epp
    flightTime = telast - (epp - j)*cwidth
    Ef = Ei*(telast/flightTime)^2
    dtdE    = flightTime/(2.*Ef)  * 1e6  ; to be concistent with vnorm absolute calcualtion
    kikf    = SQRT(Ei/Ef)                ; all times above ccalculated in sec.
    
    x_dummy = Ei - Ef
    w_out = w_out * dtdE*kikf
    e_out = e_out * dtdE*kikf
   
; ----------------------------------------------------------------------------------       
;    interpolate to have same x at all angle and data at same energies accordingly
; ----------------------------------------------------------------------------------       
    if (size(x_dummy))[0] gt 1 then begin
       sx    = size(x_dummy)
        x_out = total(x_dummy,2)/sx[2]
        for j = 0,sx[2]-1 do begin
            w_out[*,j] = interpol(w_out[*,j],x_dummy[*,j],x_out) 
            e_out[*,j] = interpol(e_out[*,j],x_dummy[*,j],x_out) 
        endfor 
    endif else x_out = x_dummy

; ----------------------------------------------------------------------------------       
; To avoid data greter than 1000 meV that are not written in the INX format data files
; ----------------------------------------------------------------------------------       
    ind = where(abs(x_out) lt 999.0)
    x_out = x_out[ind]
    w_out = w_out[ind,*]       
    e_out = e_out[ind,*]       

; ----------------------------------------------------------------------------------       
;       Copies work-space stuff from simple arrays back to main LAMP arrays
; ----------------------------------------------------------------------------------       
    mod_datp, datp, "x_tit", ' Energy Transfer (meV)'
   mod_datp, datp, "x", x_out 
   mod_datp, datp, "e", e_out
    give_datp,datp

   return, w_out
   
end