1 | # Algorithm to start Decon |
---|
2 | from mantid.simpleapi import * |
---|
3 | from mantid.api import PythonAlgorithm, AlgorithmFactory, MatrixWorkspaceProperty, FileProperty, FileAction, PropertyMode |
---|
4 | from mantid.kernel import StringListValidator, StringMandatoryValidator, Direction, logger |
---|
5 | from mantid import config |
---|
6 | import math, os.path, numpy as np |
---|
7 | |
---|
8 | class 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 |
---|
328 | AlgorithmFactory.subscribe(IndirectCylAbs) |
---|
329 | # |
---|