Ticket #11413: IndirectCylinderAbsorption.py

File IndirectCylinderAbsorption.py, 9.4 KB (added by Dan Nixon, 6 years ago)
Line 
1from mantid.simpleapi import *
2from mantid.api import DataProcessorAlgorithm, AlgorithmFactory, MatrixWorkspaceProperty, PropertyMode
3from mantid.kernel import StringMandatoryValidator, Direction, logger
4
5
6class IndirectCylinderAbsorption(DataProcessorAlgorithm):
7
8    def category(self):
9        return "Workflow\\Inelastic;PythonAlgorithms;CorrectionFunctions\\AbsorptionCorrections;Workflow\\MIDAS"
10
11
12    def summary(self):
13        return "Calculates indirect absorption corrections for a cylinder sample shape."
14
15
16    def PyInit(self):
17        self.declareProperty(MatrixWorkspaceProperty('SampleWorkspace', '', direction=Direction.Input),
18                             doc='Sample workspace.')
19        self.declareProperty(name='SampleChemicalFormula', defaultValue='', validator=StringMandatoryValidator(),
20                             doc='Sample chemical formula')
21        self.declareProperty(name='SampleRadius', defaultValue=0.0, doc='Sample radius')
22        self.declareProperty(name='SampleNumberDensity', defaultValue=0.1, doc='Sample number density')
23
24        self.declareProperty(MatrixWorkspaceProperty('CanWorkspace', '', optional=PropertyMode.Optional,
25                                                     direction=Direction.Input),
26                             doc='Container workspace.')
27        self.declareProperty(name='UseCanCorrections', defaultValue=False, doc='Use can corrections in subtraction')
28        self.declareProperty(name='CanChemicalFormula', defaultValue='', validator=StringMandatoryValidator(),
29                             doc='Can chemical formula')
30        self.declareProperty(name='CanRadius', defaultValue=0.0, doc='Can radius')
31        self.declareProperty(name='CanNumberDensity', defaultValue=0.1, doc='Can number density')
32        self.declareProperty(name='CanScaleFactor', defaultValue=1.0, doc='Scale factor to multiply can data')
33        self.declareProperty(name='Events', defaultValue=5000, doc = 'Number of neutron events')
34
35        self.declareProperty(name='Plot', defaultValue=False, doc='Plot options')
36
37        self.declareProperty(MatrixWorkspaceProperty('OutputWorkspace', '', direction=Direction.Output),
38                             doc='The output corrected workspace.')
39
40        self.declareProperty(MatrixWorkspaceProperty('CorrectionsWorkspace', '', direction=Direction.Output,
41                                                     optional=PropertyMode.Optional),
42                             doc='The corrections workspace for scattering and absorptions in sample.')
43
44
45    def PyExec(self):
46        from IndirectCommon import getEfixed, addSampleLogs
47
48        self._setup()
49        efixed = getEfixed(self._sample_ws_name)
50
51        sample_wave_ws = '__sam_wave'
52        ConvertUnits(InputWorkspace=self._sample_ws_name, OutputWorkspace=sample_wave_ws,
53                     Target='Wavelength', EMode='Indirect', EFixed=efixed)
54
55        SetSampleMaterial(sample_wave_ws, ChemicalFormula=self._sample_chemical_formula, SampleNumberDensity=self._sample_number_density)
56
57        CylinderAbsorption(InputWorkspace=sample_wave_ws,
58                           OutputWorkspace=self._ass_ws,
59                           SampleNumberDensity=self._sample_number_density,
60                           NumberOfWavelengthPoints=10,
61                           CylinderSampleHeight=3.0,
62                           CylinderSampleRadius=self._sample_radius,
63                           NumberOfSlices=1,
64                           NumberOfAnnuli=10)
65
66        plot_data = [self._output_ws, self._sample_ws_name]
67        plot_corr = [self._ass_ws]
68        group = self._ass_ws
69
70        if self._can_ws_name is not None:
71            can_wave_ws = '__can_wave'
72            ConvertUnits(InputWorkspace=self._can_ws_name, OutputWorkspace=can_wave_ws,
73                         Target='Wavelength', EMode='Indirect', EFixed=efixed)
74            if self._can_scale != 1.0:
75                logger.information('Scaling can by: ' + str(self._can_scale))
76                Scale(InputWorkspace=can_wave_ws, OutputWorkspace=can_wave_ws, Factor=self._can_scale, Operation='Multiply')
77
78            can_thickness = self._can_radius - self._sample_radius
79            logger.information('Can thickness: ' + str(can_thickness))
80
81            if self._use_can_corrections:
82                Divide(LHSWorkspace=sample_wave_ws, RHSWorkspace=self._ass_ws, OutputWorkspace=sample_wave_ws)
83
84                SetSampleMaterial(can_wave_ws, ChemicalFormula=self._can_chemical_formula, SampleNumberDensity=self._can_number_density)
85                AnnularRingAbsorption(InputWorkspace=can_wave_ws,
86                              OutputWorkspace=self._acc_ws,
87                              SampleHeight=3.0,
88                              SampleThickness=can_thickness,
89                              CanInnerRadius=0.9*self._sample_radius,
90                              CanOuterRadius=1.1*self._can_radius,
91                              SampleChemicalFormula=self._can_chemical_formula,
92                              SampleNumberDensity=self._can_number_density,
93                              NumberOfWavelengthPoints=10,
94                              EventsPerPoint=self._events)
95
96                Divide(LHSWorkspace=can_wave_ws, RHSWorkspace=self._acc_ws, OutputWorkspace=can_wave_ws)
97                Minus(LHSWorkspace=sample_wave_ws, RHSWorkspace=can_wave_ws, OutputWorkspace=sample_wave_ws)
98                plot_corr.append(self._acc_ws)
99                group += ',' + self._acc_ws
100
101            else:
102                Minus(LHSWorkspace=sample_wave_ws, RHSWorkspace=can_wave_ws, OutputWorkspace=sample_wave_ws)
103                Divide(LHSWorkspace=sample_wave_ws, RHSWorkspace=self._ass_ws, OutputWorkspace=sample_wave_ws)
104
105                DeleteWorkspace(can_wave_ws)
106                plot_data.append(self._can_ws_name)
107
108        else:
109            Divide(LHSWorkspace=sample_wave_ws, RHSWorkspace=self._ass_ws, OutputWorkspace=sample_wave_ws)
110
111        ConvertUnits(InputWorkspace=sample_wave_ws, OutputWorkspace=self._output_ws,
112                     Target='DeltaE', EMode='Indirect', EFixed=efixed)
113        DeleteWorkspace(sample_wave_ws)
114
115        sample_logs = {'sample_shape': 'cylinder',
116                       'sample_filename': self._sample_ws_name,
117                       'sample_radius': self._sample_radius}
118        addSampleLogs(self._ass_ws, sample_logs)
119        addSampleLogs(self._output_ws, sample_logs)
120
121        if self._can_ws_name is not None:
122            AddSampleLog(Workspace=self._output_ws, LogName='can_filename', LogType='String', LogText=str(self._can_ws_name))
123            AddSampleLog(Workspace=self._output_ws, LogName='can_scale', LogType='String', LogText=str(self._can_scale))
124            if self._use_can_corrections:
125                addSampleLogs(self._acc_ws, sample_logs)
126                AddSampleLog(Workspace=self._acc_ws, LogName='can_filename', LogType='String', LogText=str(self._can_ws_name))
127                AddSampleLog(Workspace=self._acc_ws, LogName='can_scale', LogType='String', LogText=str(self._can_scale))
128                AddSampleLog(Workspace=self._output_ws, LogName='can_thickness', LogType='String', LogText=str(can_thickness))
129
130        self.setProperty('OutputWorkspace', self._output_ws)
131
132        # Output the Abs group workspace if it is wanted, delete if not
133        if self._abs_ws == '':
134            DeleteWorkspace(self._ass_ws)
135            if self._can_ws_name is not None:
136                if self._use_can_corrections:
137                    DeleteWorkspace(self._acc_ws)
138        else:
139            group_name = self._abs_ws + '_abs'
140            GroupWorkspaces(InputWorkspaces=group, OutputWorkspace=group_name)
141            self.setProperty('CorrectionsWorkspace', group_name)
142
143        if self._plot:
144            from IndirectImport import import_mantidplot
145            mantid_plot = import_mantidplot()
146            mantid_plot.plotSpectrum(plot_data, 0)
147            mantid_plot.plotSpectrum(plot_corr, 0)
148
149    def _setup(self):
150        """
151        Get algorithm properties.
152        """
153
154        self._sample_ws_name = self.getPropertyValue('SampleWorkspace')
155        self._sample_chemical_formula = self.getPropertyValue('SampleChemicalFormula')
156        self._sample_number_density = self.getProperty('SampleNumberDensity').value
157        self._sample_radius = self.getProperty('SampleRadius').value
158
159        self._can_ws_name = self.getPropertyValue('CanWorkspace')
160        if self._can_ws_name == '':
161            self._can_ws_name = None
162        self._use_can_corrections = self.getProperty('UseCanCorrections').value
163        self._can_chemical_formula = self.getPropertyValue('CanChemicalFormula')
164        self._can_number_density = self.getProperty('CanNumberDensity').value
165        self._can_radius = self.getProperty('CanRadius').value
166        self._can_scale = self.getProperty('CanScaleFactor').value
167        self._events = self.getPropertyValue('Events')
168
169        self._plot = self.getProperty('Plot').value
170        self._output_ws = self.getPropertyValue('OutputWorkspace')
171
172        self._abs_ws = self.getPropertyValue('CorrectionsWorkspace')
173        if self._abs_ws == '':
174            self._ass_ws = '__ass'
175            self._acc_ws = '__acc'
176        else:
177            self._ass_ws = self._abs_ws + '_ass'
178            self._acc_ws = self._abs_ws + '_acc'
179
180
181# Register algorithm with Mantid
182AlgorithmFactory.subscribe(IndirectCylinderAbsorption)