1 | # IDA F2PY Absorption Corrections Wrapper |
---|
2 | ## Handle selection of .pyd files for absorption corrections |
---|
3 | import platform, sys |
---|
4 | from IndirectImport import * |
---|
5 | if is_supported_f2py_platform(): |
---|
6 | cylabs = import_f2py("cylabs") |
---|
7 | else: |
---|
8 | unsupported_message() |
---|
9 | |
---|
10 | from IndirectCommon import * |
---|
11 | from mantid.simpleapi import * |
---|
12 | from mantid import config, logger, mtd |
---|
13 | import math, os.path, numpy as np |
---|
14 | mp = import_mantidplot() |
---|
15 | |
---|
16 | def WaveRange(inWS, efixed): |
---|
17 | # create a list of 10 equi-spaced wavelengths spanning the input data |
---|
18 | oWS = '__WaveRange' |
---|
19 | ExtractSingleSpectrum(InputWorkspace=inWS, OutputWorkspace=oWS, WorkspaceIndex=0) |
---|
20 | ConvertUnits(InputWorkspace=oWS, OutputWorkspace=oWS, Target='Wavelength', |
---|
21 | EMode='Indirect', EFixed=efixed) |
---|
22 | Xin = mtd[oWS].readX(0) |
---|
23 | xmin = mtd[oWS].readX(0)[0] |
---|
24 | xmax = mtd[oWS].readX(0)[len(Xin)-1] |
---|
25 | ebin = 0.5 |
---|
26 | nw1 = int(xmin/ebin) |
---|
27 | nw2 = int(xmax/ebin)+1 |
---|
28 | w1 = nw1*ebin |
---|
29 | w2 = nw2*ebin |
---|
30 | wave = [] |
---|
31 | nw = 10 |
---|
32 | ebin = (w2-w1)/(nw-1) |
---|
33 | for l in range(0,nw): |
---|
34 | wave.append(w1+l*ebin) |
---|
35 | DeleteWorkspace(oWS) |
---|
36 | return wave |
---|
37 | |
---|
38 | def CheckSize(size,geom,ncan,Verbose): |
---|
39 | if geom == 'cyl': |
---|
40 | if (size[1] - size[0]) < 1e-4: |
---|
41 | error = 'Sample outer radius not > inner radius' |
---|
42 | logger.notice('ERROR *** '+error) |
---|
43 | sys.exit(error) |
---|
44 | else: |
---|
45 | if Verbose: |
---|
46 | message = 'Sam : inner radius = '+str(size[0])+' ; outer radius = '+str(size[1]) |
---|
47 | logger.notice(message) |
---|
48 | if geom == 'flt': |
---|
49 | if size[0] < 1e-4: |
---|
50 | error = 'Sample thickness is zero' |
---|
51 | logger.notice('ERROR *** '+error) |
---|
52 | sys.exit(error) |
---|
53 | else: |
---|
54 | if Verbose: |
---|
55 | logger.notice('Sam : thickness = '+str(size[0])) |
---|
56 | if ncan == 2: |
---|
57 | if geom == 'cyl': |
---|
58 | if (size[2] - size[1]) < 1e-4: |
---|
59 | error = 'Can inner radius not > sample outer radius' |
---|
60 | logger.notice('ERROR *** '+error) |
---|
61 | sys.exit(error) |
---|
62 | else: |
---|
63 | if Verbose: |
---|
64 | message = 'Can : inner radius = '+str(size[1])+' ; outer radius = '+str(size[2]) |
---|
65 | logger.notice(message) |
---|
66 | if geom == 'flt': |
---|
67 | if size[1] < 1e-4: |
---|
68 | error = 'Can thickness is zero' |
---|
69 | logger.notice('ERROR *** '+error) |
---|
70 | sys.exit(error) |
---|
71 | else: |
---|
72 | if Verbose: |
---|
73 | logger.notice('Can : thickness = '+str(size[1])) |
---|
74 | |
---|
75 | def CheckDensity(density,ncan): |
---|
76 | if density[0] < 1e-5: |
---|
77 | error = 'Sample density is zero' |
---|
78 | logger.notice('ERROR *** '+error) |
---|
79 | sys.exit(error) |
---|
80 | if ncan == 2: |
---|
81 | if density[1] < 1e-5: |
---|
82 | error = 'Can density is zero' |
---|
83 | logger.notice('ERROR *** '+error) |
---|
84 | sys.exit(error) |
---|
85 | |
---|
86 | def AbsRun(inputWS, geom, beam, ncan, size, density, sigs, siga, avar, Verbose, Save): |
---|
87 | workdir = config['defaultsave.directory'] |
---|
88 | if Verbose: |
---|
89 | logger.notice('Sample run : '+inputWS) |
---|
90 | Xin = mtd[inputWS].readX(0) |
---|
91 | if len(Xin) == 0: # check that there is data |
---|
92 | error = 'Sample file has no data' |
---|
93 | logger.notice('ERROR *** '+error) |
---|
94 | sys.exit(error) |
---|
95 | det = GetWSangles(inputWS) |
---|
96 | ndet = len(det) |
---|
97 | efixed = getEfixed(inputWS) |
---|
98 | wavelas = math.sqrt(81.787/efixed) # elastic wavelength |
---|
99 | waves = WaveRange(inputWS, efixed) # get wavelengths |
---|
100 | nw = len(waves) |
---|
101 | CheckSize(size,geom,ncan,Verbose) |
---|
102 | CheckDensity(density,ncan) |
---|
103 | run_name = getWSprefix(inputWS) |
---|
104 | if Verbose: |
---|
105 | message = 'Sam : sigt = '+str(sigs[0])+' ; siga = '+str(siga[0])+' ; rho = '+str(density[0]) |
---|
106 | logger.notice(message) |
---|
107 | if ncan == 2: |
---|
108 | message = 'Can : sigt = '+str(sigs[1])+' ; siga = '+str(siga[1])+' ; rho = '+str(density[1]) |
---|
109 | logger.notice(message) |
---|
110 | logger.notice('Elastic lambda : '+str(wavelas)) |
---|
111 | message = 'Lambda : '+str(nw)+' values from '+str(waves[0])+' to '+str(waves[nw-1]) |
---|
112 | logger.notice(message) |
---|
113 | message = 'Detector angles : '+str(ndet)+' from '+str(det[0])+' to '+str(det[ndet-1]) |
---|
114 | logger.notice(message) |
---|
115 | eZ = np.zeros(nw) # set errors to zero |
---|
116 | name = run_name + geom |
---|
117 | assWS = name + '_ass' |
---|
118 | asscWS = name + '_assc' |
---|
119 | acscWS = name + '_acsc' |
---|
120 | accWS = name + '_acc' |
---|
121 | fname = name +'_Abs' |
---|
122 | for n in range(0,ndet): |
---|
123 | if geom == 'flt': |
---|
124 | angles = [avar, det[n]] |
---|
125 | (A1,A2,A3,A4) = FlatAbs(ncan, size, density, sigs, siga, angles, waves) |
---|
126 | kill = 0 |
---|
127 | if geom == 'cyl': |
---|
128 | wrk = workdir + run_name |
---|
129 | wrk.ljust(120,' ') |
---|
130 | astep = avar |
---|
131 | if (astep) < 1e-5: |
---|
132 | error = 'Step size is zero' |
---|
133 | logger.notice('ERROR *** '+error) |
---|
134 | sys.exit(error) |
---|
135 | nstep = int((size[1] - size[0])/astep) |
---|
136 | if nstep < 20: |
---|
137 | error = 'Number of steps ( '+str(nstep)+' ) should be >= 20' |
---|
138 | logger.notice('ERROR *** '+error) |
---|
139 | sys.exit(error) |
---|
140 | angle = det[n] |
---|
141 | kill, A1, A2, A3, A4 = cylabs.cylabs(astep, beam, ncan, size, |
---|
142 | density, sigs, siga, angle, wavelas, waves, n, wrk, 0) |
---|
143 | if kill == 0: |
---|
144 | if Verbose: |
---|
145 | logger.notice('Detector '+str(n)+' at angle : '+str(det[n])+' * successful') |
---|
146 | if n == 0: |
---|
147 | dataA1 = A1 |
---|
148 | dataA2 = A2 |
---|
149 | dataA3 = A3 |
---|
150 | dataA4 = A4 |
---|
151 | eZero =eZ |
---|
152 | else: |
---|
153 | dataA1 = np.append(dataA1,A1) |
---|
154 | dataA2 = np.append(dataA2,A2) |
---|
155 | dataA3 = np.append(dataA3,A3) |
---|
156 | dataA4 = np.append(dataA4,A4) |
---|
157 | eZero = np.append(eZero,eZ) |
---|
158 | else: |
---|
159 | error = 'Detector '+str(n)+' at angle : '+str(det[n])+' *** failed : Error code '+str(kill) |
---|
160 | logger.notice('ERROR *** '+error) |
---|
161 | sys.exit(error) |
---|
162 | ## Create the workspaces |
---|
163 | dataX = waves * ndet |
---|
164 | qAxis = createQaxis(inputWS) |
---|
165 | CreateWorkspace(OutputWorkspace=assWS, DataX=dataX, DataY=dataA1, DataE=eZero, |
---|
166 | NSpec=ndet, UnitX='Wavelength', |
---|
167 | VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=qAxis) |
---|
168 | CreateWorkspace(OutputWorkspace=asscWS, DataX=dataX, DataY=dataA2, DataE=eZero, |
---|
169 | NSpec=ndet, UnitX='Wavelength', |
---|
170 | VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=qAxis) |
---|
171 | CreateWorkspace(OutputWorkspace=acscWS, DataX=dataX, DataY=dataA3, DataE=eZero, |
---|
172 | NSpec=ndet, UnitX='Wavelength', |
---|
173 | VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=qAxis) |
---|
174 | CreateWorkspace(OutputWorkspace=accWS, DataX=dataX, DataY=dataA4, DataE=eZero, |
---|
175 | NSpec=ndet, UnitX='Wavelength', |
---|
176 | VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=qAxis) |
---|
177 | ## Save output |
---|
178 | group = assWS +','+ asscWS +','+ acscWS +','+ accWS |
---|
179 | GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=fname) |
---|
180 | if Save: |
---|
181 | opath = os.path.join(workdir,fname+'.nxs') |
---|
182 | SaveNexusProcessed(InputWorkspace=fname, Filename=opath) |
---|
183 | if Verbose: |
---|
184 | logger.notice('Output file created : '+opath) |
---|
185 | if ncan > 1: |
---|
186 | return [assWS, asscWS, acscWS, accWS] |
---|
187 | else: |
---|
188 | return [assWS] |
---|
189 | |
---|
190 | def AbsRunFeeder(inputWS, geom, beam, ncan, size, density, sigs, siga, avar, |
---|
191 | plotOpt='None', Verbose=False, Save=False): |
---|
192 | Verbose = True |
---|
193 | Save = True |
---|
194 | StartTime('CalculateCorrections') |
---|
195 | '''Handles the feeding of input and plotting of output for the F2PY |
---|
196 | absorption correction routine.''' |
---|
197 | workspaces = AbsRun(inputWS, geom, beam, ncan, size, density, |
---|
198 | sigs, siga, avar, Verbose, Save) |
---|
199 | EndTime('CalculateCorrections') |
---|
200 | if ( plotOpt == 'None' ): |
---|
201 | return |
---|
202 | if ( plotOpt == 'Wavelength' or plotOpt == 'Both' ): |
---|
203 | graph = mp.plotSpectrum(workspaces, 0) |
---|
204 | if ( plotOpt == 'Angle' or plotOpt == 'Both' ): |
---|
205 | graph = mp.plotTimeBin(workspaces, 0) |
---|
206 | graph.activeLayer().setAxisTitle(mp.Layer.Bottom, 'Angle') |
---|
207 | |
---|
208 | # FlatAbs - CALCULATE FLAT PLATE ABSORPTION FACTORS |
---|
209 | # |
---|
210 | # Input parameters : |
---|
211 | # sigs - list of scattering cross-sections |
---|
212 | # siga - list of absorption cross-sections |
---|
213 | # density - list of density |
---|
214 | # ncan - =0 no can, >1 with can |
---|
215 | # thick - list of thicknesses ts,t1,t2 |
---|
216 | # angles - list of angles |
---|
217 | # waves - list of wavelengths |
---|
218 | # Output parameters : |
---|
219 | # A1 - Ass ; A2 - Assc ; A3 - Acsc ; A4 - Acc |
---|
220 | |
---|
221 | def Fact(AMU,T,SEC1,SEC2): |
---|
222 | S = AMU*T*(SEC1-SEC2) |
---|
223 | F = 1.0 |
---|
224 | if (S == 0.): |
---|
225 | F = T |
---|
226 | else: |
---|
227 | S = (1-math.exp(-S))/S |
---|
228 | F = T*S |
---|
229 | return F |
---|
230 | |
---|
231 | def FlatAbs(ncan, thick, density, sigs, siga, angles, waves): |
---|
232 | PICONV = math.pi/180. |
---|
233 | ssigs = sigs[0] #sam scatt x-sect |
---|
234 | ssiga = siga[0] #sam abs x-sect |
---|
235 | rhos = density[0] #sam density |
---|
236 | TS = thick[0] #sam thicknes |
---|
237 | T1 = thick[1] #can thickness 1 |
---|
238 | T2 = thick[2] #can thickness 2 |
---|
239 | csigs = sigs[1] #can scatt x-sect |
---|
240 | csiga = siga[1] #can abs x-sect |
---|
241 | rhoc = density[1] #can density |
---|
242 | TCAN1 = angles[0] #angle can to beam |
---|
243 | TCAN = TCAN1*PICONV |
---|
244 | THETA1 = angles[1] #THETAB value - detector angle |
---|
245 | THETA = PICONV*THETA1 |
---|
246 | |
---|
247 | AmuS1 = [] # sample & can cross sections |
---|
248 | AmuC1 = [] |
---|
249 | nlam = len(waves) |
---|
250 | for n in range(0,nlam): |
---|
251 | AS1 = ssigs + ssiga*waves[n]/1.8 |
---|
252 | AmuS1.append(AS1*rhos) |
---|
253 | if (ncan > 1): |
---|
254 | for n in range(0,nlam): |
---|
255 | AC1 = csigs + csiga*waves[n]/1.8 |
---|
256 | AmuC1.append(AC1*rhoc) |
---|
257 | else: |
---|
258 | rhoc=0. |
---|
259 | |
---|
260 | SEC1 = 1./math.cos(TCAN) |
---|
261 | TSEC=THETA1-TCAN1 # TSEC IS THE ANGLE THE SCATTERED BEAM MAKES WITH THE NORMAL TO THE SAMPLE SURFACE. |
---|
262 | A1 = [] |
---|
263 | A2 = [] |
---|
264 | A3 = [] |
---|
265 | A4 = [] |
---|
266 | if (abs(abs(TSEC)-90.0) < 1.0): # case where TSEC is close to 90. CALCULATION IS UNRELIABLE |
---|
267 | ASS = 1.0 |
---|
268 | for n in range(0,nlam): #start loop over wavelengths |
---|
269 | A1.append(ASS) |
---|
270 | A2.append(ASS) |
---|
271 | A3.append(ASS) |
---|
272 | A4.append(ASS) |
---|
273 | else: |
---|
274 | TSEC = TSEC*PICONV |
---|
275 | SEC2 = 1./math.cos(TSEC) |
---|
276 | for n in range(0,nlam): #start loop over wavelengths |
---|
277 | AMUS = AmuS1[n] |
---|
278 | FS = Fact(AMUS,TS,SEC1,SEC2) |
---|
279 | ES1=AMUS*TS*SEC1 |
---|
280 | ES2=AMUS*TS*SEC2 |
---|
281 | if (ncan > 1): |
---|
282 | AMUC = AmuC1[n] |
---|
283 | F1 = Fact(AMUC,T1,SEC1,SEC2) |
---|
284 | F2 = Fact(AMUC,T2,SEC1,SEC2) |
---|
285 | E11 = AMUC*T1*SEC1 |
---|
286 | E12 = AMUC*T1*SEC2 |
---|
287 | E21 = AMUC*T2*SEC1 |
---|
288 | E22 = AMUC*T2*SEC2 |
---|
289 | if (SEC2 < 0.): |
---|
290 | ASS=FS/TS |
---|
291 | if(ncan > 1): |
---|
292 | ASSC = ASS*math.exp(-(E11-E12)) |
---|
293 | ACC1 = F1 |
---|
294 | ACC2 = F2*math.exp(-(E11-E12)) |
---|
295 | ACSC1 = ACC1 |
---|
296 | ACSC2 = ACC2*math.exp(-(ES1-ES2)) |
---|
297 | else: |
---|
298 | ASSC = 1.0 |
---|
299 | ACSC = 1.0 |
---|
300 | ACC = 1.0 |
---|
301 | else: |
---|
302 | ASS=math.exp(-ES2)*FS/TS |
---|
303 | if(ncan > 1): |
---|
304 | ASSC = math.exp(-(E11+E22))*ASS |
---|
305 | ACC1 = math.exp(-(E12+E22))*F1 |
---|
306 | ACC2 = math.exp(-(E11+E22))*F2 |
---|
307 | ACSC1 = ACC1*math.exp(-ES2) |
---|
308 | ACSC2 = ACC2*math.exp(-ES1) |
---|
309 | else: |
---|
310 | ASSC = 1.0 |
---|
311 | ACSC = 1.0 |
---|
312 | ACC = 1.0 |
---|
313 | tsum = T1+T2 |
---|
314 | if(tsum > 0.): |
---|
315 | ACC = (ACC1+ACC2)/tsum |
---|
316 | ACSC = (ACSC1+ACSC2)/tsum |
---|
317 | else: |
---|
318 | ACC = 1.0 |
---|
319 | ACSC = 1.0 |
---|
320 | A1.append(ASS) |
---|
321 | A2.append(ASSC) |
---|
322 | A3.append(ACSC) |
---|
323 | A4.append(ACC) |
---|
324 | return A1, A2, A3, A4 |
---|