Ticket #10896: IndirectFlatAbs.py

File IndirectFlatAbs.py, 13.1 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 IndirectFlatAbs(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 thickness', defaultValue='', doc = 'Sample thickness')
23        self.declareProperty(name='Sample angle', defaultValue=0.1, doc = 'Sample angle')
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 thickness1', defaultValue='', doc = 'Can thickness1 front')
35        self.declareProperty(name='Can thickness2', defaultValue='', doc = 'Can thickness2 back')
36        self.declareProperty(name='Can scale factor', defaultValue='1.0', doc = 'Scale factor to multiply can data')
37
38        self.declareProperty(name='Verbose', defaultValue=False, doc = 'Switch Verbose Off/On')
39        self.declareProperty(name='Plot', defaultValue=False, doc = 'Plot options')
40        self.declareProperty(name='Save', defaultValue=False, doc = 'Switch Save result to nxs file Off/On')
41 
42    def PyExec(self):
43
44        from IndirectCommon import StartTime, EndTime, getEfixed, addSampleLogs
45        from IndirectAbsCor import FlatAbs
46        from IndirectImport import import_mantidplot
47        mp = import_mantidplot()
48
49        StartTime('FlatPlate Absorption')
50        workdir = config['defaultsave.directory']
51        self._setup()
52        self._waveRange()
53        swaveWS = '__sam_wave'
54        if self._diffraction:
55            ConvertUnits(InputWorkspace=self._sam, OutputWorkspace=swaveWS, Target='Wavelength')
56        else:
57            ConvertUnits(InputWorkspace=self._sam, OutputWorkspace=swaveWS, Target='Wavelength',
58            EMode='Indirect', EFixed=self._efixed)
59
60        name = self._sam[:-4] + '_flt'
61        assWS = name + '_ass'
62        SetSampleMaterial(swaveWS, ChemicalFormula=self._sam_chem, SampleNumberDensity=self._sam_density)
63        sample = mtd[swaveWS].sample()
64        sam_mat = sample.getMaterial()
65        # total scattering x-section
66        sigs = [sam_mat.totalScatterXSection()]
67        # absorption x-section
68        siga = [sam_mat.absorbXSection()]
69        size = [self._sam_thickness]
70        density = [self._sam_density]
71        ncan = 0
72        ndet = len(self._det)
73
74        if self._usecan:
75            cwaveWS = '__can_wave'
76            if self._diffraction:
77                ConvertUnits(InputWorkspace=self._can, OutputWorkspace=cwaveWS, Target='Wavelength')
78            else:
79                ConvertUnits(InputWorkspace=self._can, OutputWorkspace=cwaveWS, Target='Wavelength',
80                    EMode='Indirect', EFixed=self._efixed)
81            SetSampleMaterial(InputWorkspace=cwaveWS, ChemicalFormula=self._can_chem, SampleNumberDensity=self._can_density)
82            can_sample = mtd[cwaveWS].sample()
83            can_mat = can_sample.getMaterial()
84
85        # total scattering x-section for can
86            sigs.append(can_mat.totalScatterXSection())
87            sigs.append(can_mat.totalScatterXSection())
88        # absorption x-section for can
89            siga.append(can_mat.absorbXSection())
90            siga.append(can_mat.absorbXSection())
91            size.append(self._can_thickness1)
92            size.append(self._can_thickness2)
93            density.append(self._can_density)
94            density.append(self._can_density)
95            ncan = 2
96
97        dataA1 = []
98        dataA2 = []
99        dataA3 = []
100        dataA4 = []
101
102    #initially set errors to zero
103        eZero = np.zeros(len(self._waves))
104
105        for n in range(ndet):
106            angles = [self._sam_angle, self._det[n]]
107            (A1, A2, A3, A4) = FlatAbs(ncan, size, density, sigs, siga, angles, self._waves)
108
109            if self._verbose:
110                logger.notice('Detector ' + str(n) + ' at angle : ' + str(self._det[n]) + ' * successful')
111
112            dataA1 = np.append(dataA1, A1)
113            dataA2 = np.append(dataA2, A2)
114            dataA3 = np.append(dataA3, A3)
115            dataA4 = np.append(dataA4, A4)
116
117        sample_logs = {'sample_shape': 'flatplate', 'sample_filename': self._sam,
118                        'sample_thickness': self._sam_thickness, 'sample_angle': self._sam_angle}
119        dataX = self._waves * ndet
120
121        if self._diffraction:
122            v_axis_unit = 'dSpacing'
123            v_axis_values = [1.0]
124        else:
125            v_axis_unit = 'MomentumTransfer'
126            v_axis_values = self._q
127
128    # Create the output workspaces
129        assWS = name + '_ass'
130        asscWS = name + '_assc'
131        acscWS = name + '_acsc'
132        accWS = name + '_acc'
133        fname = name + '_abs'
134
135        CreateWorkspace(OutputWorkspace=assWS, DataX=dataX, DataY=dataA1,
136                    NSpec=ndet, UnitX='Wavelength',
137                    VerticalAxisUnit=v_axis_unit, VerticalAxisValues=v_axis_values)
138        addSampleLogs(assWS, sample_logs)
139
140        CreateWorkspace(OutputWorkspace=asscWS, DataX=dataX, DataY=dataA2,
141                    NSpec=ndet, UnitX='Wavelength',
142                    VerticalAxisUnit=v_axis_unit, VerticalAxisValues=v_axis_values)
143        addSampleLogs(asscWS, sample_logs)
144
145        CreateWorkspace(OutputWorkspace=acscWS, DataX=dataX, DataY=dataA3,
146                    NSpec=ndet, UnitX='Wavelength',
147                    VerticalAxisUnit=v_axis_unit, VerticalAxisValues=v_axis_values)
148        addSampleLogs(acscWS, sample_logs)
149
150        CreateWorkspace(OutputWorkspace=accWS, DataX=dataX, DataY=dataA4,
151                    NSpec=ndet, UnitX='Wavelength',
152                    VerticalAxisUnit=v_axis_unit, VerticalAxisValues=v_axis_values)
153        addSampleLogs(accWS, sample_logs)
154
155        if self._usecan:
156            workspaces = [assWS, asscWS, acscWS, accWS]
157            AddSampleLog(Workspace=assWS, LogName='can_filename', LogType='String', LogText=str(self._can))
158            AddSampleLog(Workspace=asscWS, LogName='can_filename', LogType='String', LogText=str(self._can))
159            AddSampleLog(Workspace=acscWS, LogName='can_filename', LogType='String', LogText=str(self._can))
160            AddSampleLog(Workspace=accWS, LogName='can_filename', LogType='String', LogText=str(self._can))
161        else:
162            workspaces = [assWS]
163
164        group = assWS + ',' + asscWS + ',' + acscWS + ',' + accWS
165        GroupWorkspaces(InputWorkspaces=group, OutputWorkspace=fname)
166
167        if self._plot:
168            graph1 = mp.plotSpectrum(workspaces, 0)
169            graph2 = mp.plotTimeBin(workspaces, 0)
170            graph2.activeLayer().setAxisTitle(mp.Layer.Bottom, 'Angle')
171
172        if self._save:
173            path = os.path.join(workdir,assWS + '.nxs')
174            SaveNexusProcessed(InputWorkspace=assWS, Filename=path)
175            if self._verbose:
176                logger.notice('Output file created : '+path)
177
178        EndTime('FlatPlate Absorption')
179
180    def _setup(self):
181        self._verbose = self.getProperty('Verbose').value
182        sInput = self.getPropertyValue('Sample Input')
183        if sInput == 'Workspace':
184            s_ws = self.getPropertyValue('Sample Workspace')
185        else:
186            s_ws = ''
187        if sInput == 'File':
188            s_file = self.getPropertyValue('Sample File')
189        else:
190            s_file = ''
191        self._input = sInput
192        self._path = s_file
193        self._ws = s_ws
194        self._getData()
195        self._sam = self._name
196        self._sam_chem = self.getPropertyValue('Sample chemical formula')
197        self._sam_density = float(self.getPropertyValue('Sample number density'))
198        self._sam_thickness = float(self.getPropertyValue('Sample thickness'))
199        self._sam_angle = float(self.getPropertyValue('Sample angle'))
200
201        self._usecan = self.getProperty('Use can').value
202        if self._usecan:
203            cInput = self.getPropertyValue('Can Input')
204            if cInput == 'Workspace':
205                c_ws = self.getPropertyValue('Can Workspace')
206            else:
207                c_ws = ''
208            if cInput == 'File':
209                c_file = self.getPropertyValue('Can File')
210            else:
211                c_file = ''
212            self._input = cInput
213            self._path = c_file
214            self._ws = c_ws
215            self._getData()
216            self._can = self._name
217            self._can_chem = self.getPropertyValue('Can chemical formula')
218            self._can_density = float(self.getPropertyValue('Can number density'))
219            self._can_thickness1 = float(self.getPropertyValue('Can thickness1'))
220            self._can_thickness2 = float(self.getPropertyValue('Can thickness2'))
221            self._can_scale = self.getPropertyValue('Can scale factor')
222
223        self._plot = self.getProperty('Plot').value
224        self._save = self.getProperty('Save').value
225               
226    def _getData(self):   #get data
227        if self._input == 'Workspace':
228            inWS = self._ws
229            self._name = inWS
230            if self._verbose:
231                logger.notice('Input from Workspace : '+inWS)
232        elif self._input == 'File':
233            self._getFileName()
234            inWS = self._name
235            LoadNexus(Filename=self._path, OutputWorkspace=inWS)
236        else:
237            raise ValueError('Input type not defined')
238
239    def _getFileName(self):
240        import os.path
241        path = self._path
242        if(os.path.isfile(path)): 
243            base = os.path.basename(path)
244            self._name = os.path.splitext(base)[0]
245            ext = os.path.splitext(base)[1]
246            if self._verbose:
247                logger.notice('Input file : '+path)
248        else:
249            raise ValueError('Could not find file: ' + path)
250
251    def _waveRange(self):
252        from IndirectCommon import checkUnitIs, GetWSangles, getEfixed, GetThetaQ
253# create a list of 10 equi-spaced wavelengths spanning the input data
254        oWS = '__WaveRange'
255        ExtractSingleSpectrum(InputWorkspace=self._sam, OutputWorkspace=oWS, WorkspaceIndex=0)
256        self._diffraction = checkUnitIs(self._sam, 'dSpacing')
257
258        if self._diffraction:
259            self._det = GetWSangles(self._sam)
260            self._efixed = 0.0
261            ConvertUnits(InputWorkspace=oWS, OutputWorkspace=oWS, Target='Wavelength',
262                     EMode='Elastic')
263        else:
264            self._det, self._q = GetThetaQ(self._sam)
265            self._efixed = getEfixed(self._sam)
266            ConvertUnits(InputWorkspace=oWS, OutputWorkspace=oWS, Target='Wavelength',
267                     EMode='Indirect', EFixed=self._efixed)
268        Xin = mtd[oWS].readX(0)
269        xmin = mtd[oWS].readX(0)[0]
270        xmax = mtd[oWS].readX(0)[len(Xin) - 1]
271        ebin = 0.5
272        nw1 = int(xmin/ebin)
273        nw2 = int(xmax/ebin)+1
274        w1 = nw1*ebin
275        w2 = nw2*ebin
276        waves = []
277        nw = 10
278        ebin = (w2-w1)/(nw-1)
279        for l in range(0,nw):
280            waves.append(w1+l*ebin)
281        DeleteWorkspace(oWS)
282        self._waves = waves
283
284        if self._diffraction:
285            self._wavelas = waves[int(nw / 2)]
286        else:
287            self._wavelas = math.sqrt(81.787/self._efixed) # elastic wavelength
288
289        if self._verbose:
290            logger.notice('Elastic lambda : ' + str(self._wavelas))
291            nw = len(self._waves)
292            message = 'Lambda : ' + str(nw) + ' values from ' + str(self._waves[0]) + ' to ' + str(self._waves[nw - 1])
293            logger.notice(message)
294            ndet = len(self._det)
295            message = 'Detector angles : ' + str(ndet) + ' from ' + str(self._det[0]) + ' to ' + str(self._det[ndet - 1])
296            logger.notice(message)
297
298# Register algorithm with Mantid
299AlgorithmFactory.subscribe(IndirectFlatAbs)
300#