Ticket #10896: IndirectCylAbs.py

File IndirectCylAbs.py, 14.2 KB (added by Dan Nixon, 6 years ago)
Line 
1# Algorithm to start Decon
2from mantid.simpleapi import *
3from mantid.api import PythonAlgorithm, AlgorithmFactory, MatrixWorkspaceProperty, FileProperty, FileAction, PropertyMode
4from mantid.kernel import StringListValidator, StringMandatoryValidator, Direction, logger
5from mantid import config
6import math, os.path, numpy as np
7
8class IndirectCylAbs(PythonAlgorithm):
9 
10    def category(self):
11        return "Workflow\\MIDAS;PythonAlgorithms"
12
13    def PyInit(self):
14        self.declareProperty(name='Sample Input', defaultValue='Workspace', validator=StringListValidator(['Workspace','File']),
15            doc='Sample input type')
16        self.declareProperty(MatrixWorkspaceProperty('Sample Workspace', '', optional=PropertyMode.Optional, 
17            direction=Direction.Input), doc="Name for the input Sample workspace.")
18        self.declareProperty(FileProperty('Sample File', '', action=FileAction.OptionalLoad, extensions=["_red.nxs"]),
19                             doc='File path for Sample file')
20        self.declareProperty(name='Sample chemical formula', defaultValue='', doc = 'Sample chemical formula')
21        self.declareProperty(name='Sample number density', defaultValue='', doc = 'Sample number density')
22        self.declareProperty(name='Sample radius1', defaultValue='', doc = 'Sample radius1')
23        self.declareProperty(name='Sample radius2', defaultValue='', doc = 'Sample radius2')
24
25        self.declareProperty(name='Use Can', defaultValue=False, doc = 'Use Can')
26        self.declareProperty(name='Can Input', defaultValue='Workspace', validator=StringListValidator(['Workspace','File']),
27            doc='Can input type')
28        self.declareProperty(MatrixWorkspaceProperty('Can Workspace', '', optional=PropertyMode.Optional, 
29            direction=Direction.Input), doc="Name for the input Can workspace.")
30        self.declareProperty(FileProperty('Can File', '', action=FileAction.OptionalLoad, extensions=["_red.nxs"]),
31                             doc='File path for Can file')
32        self.declareProperty(name='Can chemical formula', defaultValue='', doc = 'Can chemical formula')
33        self.declareProperty(name='Can number density', defaultValue='', doc = 'Can number density')
34        self.declareProperty(name='Can radius3', defaultValue='', doc = 'Can radius3')
35        self.declareProperty(name='Can scale factor', defaultValue='1.0', doc = 'Scale factor to multiply can data')
36
37        self.declareProperty(name='Step size', defaultValue='', doc = 'Step size')
38        self.declareProperty(name='Beam width', defaultValue=2.0, doc = 'Beam width')
39        self.declareProperty(name='Beam height', defaultValue=3.0, doc = 'Beam height')
40
41        self.declareProperty(name='Verbose', defaultValue=False, doc = 'Switch Verbose Off/On')
42        self.declareProperty(name='Plot', defaultValue=False, doc = 'Plot options')
43        self.declareProperty(name='Save', defaultValue=False, doc = 'Switch Save result to nxs file Off/On')
44 
45    def PyExec(self):
46
47        from IndirectCommon import StartTime, EndTime, getEfixed, addSampleLogs, getDefaultWorkingDirectory
48        from IndirectAbsCor import FlatAbs
49        from IndirectImport import import_mantidplot, is_supported_f2py_platform, import_f2py, unsupported_message
50        if is_supported_f2py_platform():
51            cylabs = import_f2py("cylabs")
52        else:
53            unsupported_message()
54        mp = import_mantidplot()
55
56        StartTime('Cylinder Absorption')
57        workdir = config['defaultsave.directory']
58        self._setup()
59        self._waveRange()
60        swaveWS = '__sam_wave'
61        if self._diffraction:
62            ConvertUnits(InputWorkspace=self._sam, OutputWorkspace=swaveWS, Target='Wavelength')
63        else:
64            ConvertUnits(InputWorkspace=self._sam, OutputWorkspace=swaveWS, Target='Wavelength',
65            EMode='Indirect', EFixed=self._efixed)
66
67        name = self._sam[:-4] + '_cyl'
68#        wrk = getDefaultWorkingDirectory() + name
69        wrk = 'CylAbs'
70        wrk = wrk.ljust(120,' ')
71        len_wrk = 0
72        assWS = name + '_ass'
73        SetSampleMaterial(swaveWS, ChemicalFormula=self._sam_chem, SampleNumberDensity=self._sam_density)
74        sample = mtd[swaveWS].sample()
75        sam_mat = sample.getMaterial()
76        # total scattering x-section
77        sigs = [sam_mat.totalScatterXSection()]
78        # absorption x-section
79        siga = [sam_mat.absorbXSection()]
80        radii = [self._sam_rad1, self._sam_rad2]
81        density = [self._sam_density]
82        ndet = len(self._det)
83
84        if self._usecan:
85            cwaveWS = '__can_wave'
86            if self._diffraction:
87                ConvertUnits(InputWorkspace=self._can, OutputWorkspace=cwaveWS, Target='Wavelength')
88            else:
89                ConvertUnits(InputWorkspace=self._can, OutputWorkspace=cwaveWS, Target='Wavelength',
90                    EMode='Indirect', EFixed=self._efixed)
91            SetSampleMaterial(InputWorkspace=cwaveWS, ChemicalFormula=self._can_chem, SampleNumberDensity=self._can_density)
92            can_sample = mtd[cwaveWS].sample()
93            can_mat = can_sample.getMaterial()
94
95        # total scattering x-section for can
96            sigs.append(can_mat.totalScatterXSection())
97            sigs.append(can_mat.totalScatterXSection())
98        # absorption x-section for can
99            siga.append(can_mat.absorbXSection())
100            siga.append(can_mat.absorbXSection())
101            radii.append(self._can_rad3)
102            radii.append(self._can_rad3)
103            density.append(self._can_density)
104            density.append(self._can_density)
105            ncan = 2
106        else:
107            ncan = 0
108            sigs.append(0.0)
109            sigs.append(0.0)
110            siga.append(0.0)
111            siga.append(0.0)
112            radii.append(0.0)
113            radii.append(0.0)
114            density.append(0.0)
115            density.append(0.0)
116
117        dataA1 = []
118        dataA2 = []
119        dataA3 = []
120        dataA4 = []
121
122    #initially set errors to zero
123        eZero = np.zeros(len(self._waves))
124
125        for n in range(ndet):
126            kill, A1, A2, A3, A4 = cylabs.cylabs(self._step_size, self._beam, ncan, radii,
127                density, sigs, siga, self._det[n], self._wavelas, self._waves, n, wrk, len_wrk)
128#      real step, beam(9), radii(4)
129#      real density(3),sigs(3),siga(3)
130#      real wavelas, angle, waves(10)
131#      integer ncan, n, len_wrk
132#      character*120 wrk
133
134            if self._verbose:
135                logger.notice('Detector ' + str(n) + ' at angle : ' + str(self._det[n]) + ' * successful')
136
137            dataA1 = np.append(dataA1, A1)
138            dataA2 = np.append(dataA2, A2)
139            dataA3 = np.append(dataA3, A3)
140            dataA4 = np.append(dataA4, A4)
141
142        sample_logs = {'sample_shape': 'cylinder', 'sample_filename': self._sam,
143                        'sample_rad1': self._sam_rad1, 'sample_rad2': self._sam_rad2}
144        dataX = self._waves * ndet
145
146        if self._diffraction:
147            v_axis_unit = 'dSpacing'
148            v_axis_values = [1.0]
149        else:
150            v_axis_unit = 'MomentumTransfer'
151            v_axis_values = self._q
152
153    # Create the output workspaces
154        assWS = name + '_ass'
155        asscWS = name + '_assc'
156        acscWS = name + '_acsc'
157        accWS = name + '_acc'
158        fname = name + '_abs'
159
160        CreateWorkspace(OutputWorkspace=assWS, DataX=dataX, DataY=dataA1,
161                    NSpec=ndet, UnitX='Wavelength',
162                    VerticalAxisUnit=v_axis_unit, VerticalAxisValues=v_axis_values)
163        addSampleLogs(assWS, sample_logs)
164
165        CreateWorkspace(OutputWorkspace=asscWS, DataX=dataX, DataY=dataA2,
166                    NSpec=ndet, UnitX='Wavelength',
167                    VerticalAxisUnit=v_axis_unit, VerticalAxisValues=v_axis_values)
168        addSampleLogs(asscWS, sample_logs)
169
170        CreateWorkspace(OutputWorkspace=acscWS, DataX=dataX, DataY=dataA3,
171                    NSpec=ndet, UnitX='Wavelength',
172                    VerticalAxisUnit=v_axis_unit, VerticalAxisValues=v_axis_values)
173        addSampleLogs(acscWS, sample_logs)
174
175        CreateWorkspace(OutputWorkspace=accWS, DataX=dataX, DataY=dataA4,
176                    NSpec=ndet, UnitX='Wavelength',
177                    VerticalAxisUnit=v_axis_unit, VerticalAxisValues=v_axis_values)
178        addSampleLogs(accWS, sample_logs)
179
180        if self._usecan:
181            workspaces = [assWS, asscWS, acscWS, accWS]
182            AddSampleLog(Workspace=assWS, LogName='can_filename', LogType='String', LogText=str(self._can))
183            AddSampleLog(Workspace=asscWS, LogName='can_filename', LogType='String', LogText=str(self._can))
184            AddSampleLog(Workspace=acscWS, LogName='can_filename', LogType='String', LogText=str(self._can))
185            AddSampleLog(Workspace=accWS, LogName='can_filename', LogType='String', LogText=str(self._can))
186        else:
187            workspaces = [assWS]
188
189        group = assWS + ',' + asscWS + ',' + acscWS + ',' + accWS
190        GroupWorkspaces(InputWorkspaces=group, OutputWorkspace=fname)
191
192        if self._plot:
193            graph1 = mp.plotSpectrum(workspaces, 0)
194            graph2 = mp.plotTimeBin(workspaces, 0)
195            graph2.activeLayer().setAxisTitle(mp.Layer.Bottom, 'Angle')
196
197        if self._save:
198            path = os.path.join(workdir,assWS + '.nxs')
199            SaveNexusProcessed(InputWorkspace=assWS, Filename=path)
200            if self._verbose:
201                logger.notice('Output file created : '+path)
202
203        EndTime('Cylinder Absorption')
204
205    def _setup(self):
206        self._verbose = self.getProperty('Verbose').value
207        sInput = self.getPropertyValue('Sample Input')
208        if sInput == 'Workspace':
209            s_ws = self.getPropertyValue('Sample Workspace')
210        else:
211            s_ws = ''
212        if sInput == 'File':
213            s_file = self.getPropertyValue('Sample File')
214        else:
215            s_file = ''
216        self._input = sInput
217        self._path = s_file
218        self._ws = s_ws
219        self._getData()
220        self._sam = self._name
221        self._sam_chem = self.getPropertyValue('Sample chemical formula')
222        self._sam_density = float(self.getPropertyValue('Sample number density'))
223        self._sam_rad1 = float(self.getPropertyValue('Sample radius1'))
224        self._sam_rad2 = float(self.getPropertyValue('Sample radius2'))
225
226        self._usecan = self.getProperty('Use can').value
227        if self._usecan:
228            cInput = self.getPropertyValue('Can Input')
229            if cInput == 'Workspace':
230                c_ws = self.getPropertyValue('Can Workspace')
231            else:
232                c_ws = ''
233            if cInput == 'File':
234                c_file = self.getPropertyValue('Can File')
235            else:
236                c_file = ''
237            self._input = cInput
238            self._path = c_file
239            self._ws = c_ws
240            self._getData()
241            self._can = self._name
242            self._can_chem = self.getPropertyValue('Can chemical formula')
243            self._can_density = float(self.getPropertyValue('Can number density'))
244            self._can_rad3 = float(self.getPropertyValue('Can radius3'))
245            self._can_scale = self.getPropertyValue('Can scale factor')
246
247        self._step_size = float(self.getPropertyValue('Step size'))
248        beam_width = float(self.getPropertyValue('Beam width'))
249        beam_height = float(self.getPropertyValue('Beam height'))
250        self._beam = [beam_height, 0.5 * beam_width, -0.5 * beam_width, (beam_width / 2), -(beam_width / 2), 0.0, beam_height, 0.0, beam_height]
251
252        self._plot = self.getProperty('Plot').value
253        self._save = self.getProperty('Save').value
254               
255    def _getData(self):   #get data
256        if self._input == 'Workspace':
257            inWS = self._ws
258            self._name = inWS
259            if self._verbose:
260                logger.notice('Input from Workspace : '+inWS)
261        elif self._input == 'File':
262            self._getFileName()
263            inWS = self._name
264            LoadNexus(Filename=self._path, OutputWorkspace=inWS)
265        else:
266            raise ValueError('Input type not defined')
267
268    def _getFileName(self):
269        import os.path
270        path = self._path
271        if(os.path.isfile(path)): 
272            base = os.path.basename(path)
273            self._name = os.path.splitext(base)[0]
274            ext = os.path.splitext(base)[1]
275            if self._verbose:
276                logger.notice('Input file : '+path)
277        else:
278            raise ValueError('Could not find file: ' + path)
279
280    def _waveRange(self):
281        from IndirectCommon import checkUnitIs, GetWSangles, getEfixed, GetThetaQ
282# create a list of 10 equi-spaced wavelengths spanning the input data
283        oWS = '__WaveRange'
284        ExtractSingleSpectrum(InputWorkspace=self._sam, OutputWorkspace=oWS, WorkspaceIndex=0)
285        self._diffraction = checkUnitIs(self._sam, 'dSpacing')
286
287        if self._diffraction:
288            self._det = GetWSangles(self._sam)
289            self._efixed = 0.0
290            ConvertUnits(InputWorkspace=oWS, OutputWorkspace=oWS, Target='Wavelength',
291                     EMode='Elastic')
292        else:
293            self._det, self._q = GetThetaQ(self._sam)
294            self._efixed = getEfixed(self._sam)
295            ConvertUnits(InputWorkspace=oWS, OutputWorkspace=oWS, Target='Wavelength',
296                     EMode='Indirect', EFixed=self._efixed)
297        Xin = mtd[oWS].readX(0)
298        xmin = mtd[oWS].readX(0)[0]
299        xmax = mtd[oWS].readX(0)[len(Xin) - 1]
300        ebin = 0.5
301        nw1 = int(xmin/ebin)
302        nw2 = int(xmax/ebin)+1
303        w1 = nw1*ebin
304        w2 = nw2*ebin
305        waves = []
306        nw = 10
307        ebin = (w2-w1)/(nw-1)
308        for l in range(0,nw):
309            waves.append(w1+l*ebin)
310        DeleteWorkspace(oWS)
311        self._waves = waves
312
313        if self._diffraction:
314            self._wavelas = waves[int(nw / 2)]
315        else:
316            self._wavelas = math.sqrt(81.787/self._efixed) # elastic wavelength
317
318        if self._verbose:
319            logger.notice('Elastic lambda : ' + str(self._wavelas))
320            nw = len(self._waves)
321            message = 'Lambda : ' + str(nw) + ' values from ' + str(self._waves[0]) + ' to ' + str(self._waves[nw - 1])
322            logger.notice(message)
323            ndet = len(self._det)
324            message = 'Detector angles : ' + str(ndet) + ' from ' + str(self._det[0]) + ' to ' + str(self._det[ndet - 1])
325            logger.notice(message)
326
327# Register algorithm with Mantid
328AlgorithmFactory.subscribe(IndirectCylAbs)
329#