Ticket #7976: IndirectAbsCor.py

File IndirectAbsCor.py, 12.2 KB (added by Samuel Jackson, 7 years ago)

Modified version of the python script

Line 
1# IDA F2PY Absorption Corrections Wrapper
2## Handle selection of .pyd files for absorption corrections
3import platform, sys
4from IndirectImport import *
5if is_supported_f2py_platform():
6    cylabs = import_f2py("cylabs")
7else:
8    unsupported_message()
9
10from IndirectCommon import *
11from mantid.simpleapi import *
12from mantid import config, logger, mtd
13import math, os.path, numpy as np
14mp = import_mantidplot()
15
16def 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
38def 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
75def 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
86def 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
190def 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
221def 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
231def 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