diff --git a/Code/Mantid/Framework/PythonInterface/plugins/algorithms/CylinderPaalmanPingsCorrection.py b/Code/Mantid/Framework/PythonInterface/plugins/algorithms/CylinderPaalmanPingsCorrection.py index 24e4fc0384c0..e56af5905cef 100644 --- a/Code/Mantid/Framework/PythonInterface/plugins/algorithms/CylinderPaalmanPingsCorrection.py +++ b/Code/Mantid/Framework/PythonInterface/plugins/algorithms/CylinderPaalmanPingsCorrection.py @@ -1,13 +1,16 @@ -#pylint: disable=no-init +#pylint: disable=no-init,too-many-locals,too-many-instance-attributes + from mantid.simpleapi import * -from mantid.api import PythonAlgorithm, AlgorithmFactory, PropertyMode, MatrixWorkspaceProperty, \ - WorkspaceGroupProperty, InstrumentValidator, WorkspaceUnitValidator -from mantid.kernel import StringListValidator, StringMandatoryValidator, \ - FloatBoundedValidator, Direction, logger, CompositeValidator +from mantid.api import (PythonAlgorithm, AlgorithmFactory, PropertyMode, MatrixWorkspaceProperty, + WorkspaceGroupProperty, InstrumentValidator, WorkspaceUnitValidator) +from mantid.kernel import (StringListValidator, StringMandatoryValidator, + FloatBoundedValidator, Direction, logger, CompositeValidator) from mantid import config -import math, numpy as np -#pylint: disable=too-many-instance-attributes +import math +import numpy as np + + class CylinderPaalmanPingsCorrection(PythonAlgorithm): _sample_ws_name = None @@ -25,29 +28,34 @@ class CylinderPaalmanPingsCorrection(PythonAlgorithm): _emode = None _efixed = 0.0 _output_ws_name = None - _elastic = None _use_can = None - _angles = None - _beam_width = None _beam_height = None - _waves = None + _beam_width = None _interpolate = None + _angles = None + _waves = None + _elastic = None + +#------------------------------------------------------------------------------ def category(self): return "Workflow\\MIDAS;PythonAlgorithms;CorrectionFunctions\\AbsorptionCorrections" - def summary(self): return "Calculates absorption corrections for a cylindrical or annular sample using Paalman & Pings format." +#------------------------------------------------------------------------------ def PyInit(self): ws_validator = CompositeValidator([WorkspaceUnitValidator('Wavelength'), InstrumentValidator()]) - self.declareProperty(MatrixWorkspaceProperty('SampleWorkspace', '', direction=Direction.Input, validator=ws_validator), + self.declareProperty(MatrixWorkspaceProperty('SampleWorkspace', '', + direction=Direction.Input, + validator=ws_validator), doc='Name for the input sample workspace') - self.declareProperty(name='SampleChemicalFormula', defaultValue='',validator=StringMandatoryValidator(), + self.declareProperty(name='SampleChemicalFormula', defaultValue='', + validator=StringMandatoryValidator(), doc='Sample chemical formula') self.declareProperty(name='SampleNumberDensity', defaultValue=0.1, validator=FloatBoundedValidator(0.0), @@ -58,9 +66,9 @@ def PyInit(self): doc='Sample outer radius') self.declareProperty(MatrixWorkspaceProperty('CanWorkspace', '', - direction=Direction.Input, - optional=PropertyMode.Optional, - validator=ws_validator), + direction=Direction.Input, + optional=PropertyMode.Optional, + validator=ws_validator), doc="Name for the input container workspace") self.declareProperty(name='CanChemicalFormula', defaultValue='', @@ -71,14 +79,13 @@ def PyInit(self): self.declareProperty(name='CanOuterRadius', defaultValue=0.15, doc='Can outer radius') - self.declareProperty(name='BeamHeight', defaultValue=0.1, + self.declareProperty(name='BeamHeight', defaultValue=3.0, doc='Height of the beam at the sample.') - self.declareProperty(name='BeamWidth', defaultValue=0.1, + self.declareProperty(name='BeamWidth', defaultValue=2.0, doc='Width of the beam at the sample.') self.declareProperty(name='StepSize', defaultValue=0.002, doc='Step size for calculation') - self.declareProperty(name='Interpolate', defaultValue=True, doc='Interpolate the correction workspaces to match the sample workspace') @@ -88,10 +95,12 @@ def PyInit(self): self.declareProperty(name='Efixed', defaultValue=1.0, doc='Analyser energy') - self.declareProperty(WorkspaceGroupProperty('OutputWorkspace', '',direction=Direction.Output), + self.declareProperty(WorkspaceGroupProperty('OutputWorkspace', '', + direction=Direction.Output), doc='The output corrections workspace group') - #pylint: disable=too-many-locals +#------------------------------------------------------------------------------ + def PyExec(self): from IndirectImport import is_supported_f2py_platform, import_f2py @@ -156,11 +165,11 @@ def PyExec(self): for angle_idx in range(number_angles): kill, ass, assc, acsc, acc = cylabs.cylabs(self._step_size, beam, ncan, radii, - density, sigs, siga, self._angles[angle_idx], - self._elastic, self._waves, angle_idx, wrk, 0) + density, sigs, siga, self._angles[angle_idx], self._elastic, self._waves, angle_idx, wrk, 0) if kill == 0: logger.information('Angle %d: %f successful' % (angle_idx+1, self._angles[angle_idx])) + data_ass = np.append(data_ass, ass) data_assc = np.append(data_assc, assc) data_acsc = np.append(data_acsc, acsc) @@ -170,7 +179,7 @@ def PyExec(self): raise ValueError('Angle ' + str(angle_idx) + ' : ' + str(self._angles[angle_idx]) + ' *** failed : Error code ' + str(kill)) sample_logs = {'sample_shape': 'cylinder', 'sample_filename': self._sample_ws_name, - 'sample_inner_radius': self._sample_inner_radius, 'sample_outer_radius': self._sample_outer_radius} + 'sample_inner_radius': self._sample_inner_radius, 'sample_outer_radius': self._sample_outer_radius} dataX = self._waves * number_angles # Create the output workspaces @@ -214,6 +223,7 @@ def PyExec(self): GroupWorkspaces(InputWorkspaces=','.join(workspaces), OutputWorkspace=self._output_ws_name) self.setPropertyValue('OutputWorkspace', self._output_ws_name) +#------------------------------------------------------------------------------ def validateInputs(self): self._setup() @@ -232,6 +242,7 @@ def validateInputs(self): return issues +#------------------------------------------------------------------------------ def _setup(self): """ @@ -267,6 +278,7 @@ def _setup(self): self._output_ws_name = self.getPropertyValue('OutputWorkspace') +#------------------------------------------------------------------------------ def _get_angles(self): """ @@ -283,6 +295,7 @@ def _get_angles(self): two_theta = detector.getTwoTheta(sample_pos, beam_pos) * 180.0 / math.pi # calc angle self._angles.append(two_theta) +#------------------------------------------------------------------------------ def _wave_range(self): wave_range = '__WaveRange' @@ -305,6 +318,7 @@ def _wave_range(self): logger.information('Elastic lambda %f' % self._elastic) DeleteWorkspace(wave_range) +#------------------------------------------------------------------------------ def _interpolate_corrections(self, workspaces): """ @@ -314,12 +328,13 @@ def _interpolate_corrections(self, workspaces): @param workspaces List of correction workspaces to interpolate """ - for wks in workspaces: + for wrksp in workspaces: SplineInterpolation(WorkspaceToMatch=self._sample_ws_name, - WorkspaceToInterpolate=wks, - OutputWorkspace=wks, + WorkspaceToInterpolate=wrksp, + OutputWorkspace=wrksp, OutputWorkspaceDeriv='') +#------------------------------------------------------------------------------ def _copy_detector_table(self, workspaces): """ @@ -330,22 +345,23 @@ def _copy_detector_table(self, workspaces): instrument = mtd[self._sample_ws_name].getInstrument().getName() - for wks in workspaces: - LoadInstrument(Workspace=wks, + for wrksp in workspaces: + LoadInstrument(Workspace=wrksp, InstrumentName=instrument) CopyDetectorMapping(WorkspaceToMatch=self._sample_ws_name, - WorkspaceToRemap=wks, + WorkspaceToRemap=wrksp, IndexBySpectrumNumber=True) +#------------------------------------------------------------------------------ - def _add_sample_logs(self, wks, sample_logs): + def _add_sample_logs(self, wrksp, sample_logs): """ Add a dictionary of logs to a workspace. The type of the log is inferred by the type of the value passed to the log. - @param wks - workspace to add logs too. + @param wrksp - workspace to add logs too. @param sample_logs - dictionary of logs to append to the workspace. """ @@ -357,8 +373,9 @@ def _add_sample_logs(self, wks, sample_logs): else: log_type = 'String' - AddSampleLog(Workspace=wks, LogName=key, LogType=log_type, LogText=str(value)) + AddSampleLog(Workspace=wrksp, LogName=key, LogType=log_type, LogText=str(value)) +#------------------------------------------------------------------------------ # Register algorithm with Mantid AlgorithmFactory.subscribe(CylinderPaalmanPingsCorrection) diff --git a/Code/Mantid/Framework/PythonInterface/plugins/algorithms/CylinderPaalmanPingsCorrection2.py b/Code/Mantid/Framework/PythonInterface/plugins/algorithms/CylinderPaalmanPingsCorrection2.py new file mode 100644 index 000000000000..46b7df2f9346 --- /dev/null +++ b/Code/Mantid/Framework/PythonInterface/plugins/algorithms/CylinderPaalmanPingsCorrection2.py @@ -0,0 +1,575 @@ +#pylint: disable=no-init,too-many-locals,too-many-instance-attributes,too-many-arguments,invalid-name + +from mantid.simpleapi import * +from mantid.api import (PythonAlgorithm, AlgorithmFactory, PropertyMode, MatrixWorkspaceProperty, + WorkspaceGroupProperty, InstrumentValidator, WorkspaceUnitValidator) +from mantid.kernel import (StringListValidator, StringMandatoryValidator, IntBoundedValidator, + FloatBoundedValidator, Direction, logger, CompositeValidator) +import math +import numpy as np + + +class CylinderPaalmanPingsCorrection(PythonAlgorithm): + + _sample_ws_name = None + _sample_chemical_formula = None + _sample_number_density = None + _sample_inner_radius = None + _sample_outer_radius = None + _use_can = False + _can_ws_name = None + _can_chemical_formula = None + _can_number_density = None + _can_outer_radius = None + _number_can = 1 + _ms = 1 + _number_wavelengths = 10 + _emode = None + _efixed = 0.0 + _step_size = None + _output_ws_name = None + _beam = list() + _angles = list() + _waves = list() + _elastic = 0.0 + _sig_s = None + _sig_a = None + _density = None + _radii = None + _interpolate = False + +#------------------------------------------------------------------------------ + + def version(self): + return 2 + + def category(self): + return "Workflow\\MIDAS;PythonAlgorithms;CorrectionFunctions\\AbsorptionCorrections" + + def summary(self): + return "Calculates absorption corrections for a cylindrical or annular sample using Paalman & Pings format." + +#------------------------------------------------------------------------------ + + def PyInit(self): + ws_validator = CompositeValidator([WorkspaceUnitValidator('Wavelength'), InstrumentValidator()]) + + self.declareProperty(MatrixWorkspaceProperty('SampleWorkspace', '', + validator=ws_validator, + direction=Direction.Input), + doc="Name for the input Sample workspace.") + + self.declareProperty(name='SampleChemicalFormula', defaultValue='', + validator=StringMandatoryValidator(), + doc='Sample chemical formula') + self.declareProperty(name='SampleNumberDensity', defaultValue=0.1, + validator=FloatBoundedValidator(0.0), + doc='Sample number density') + self.declareProperty(name='SampleInnerRadius', defaultValue=0.05, + validator=FloatBoundedValidator(0.0), + doc='Sample inner radius') + self.declareProperty(name='SampleOuterRadius', defaultValue=0.1, + validator=FloatBoundedValidator(0.0), + doc='Sample outer radius') + + self.declareProperty(MatrixWorkspaceProperty('CanWorkspace', '', + optional=PropertyMode.Optional, + validator=ws_validator, + direction=Direction.Input), + doc="Name for the input Can workspace.") + + self.declareProperty(name='CanChemicalFormula', defaultValue='', + doc='Can chemical formula') + self.declareProperty(name='CanNumberDensity', defaultValue=0.1, + validator=FloatBoundedValidator(0.0), + doc='Can number density') + self.declareProperty(name='CanOuterRadius', defaultValue=0.15, + validator=FloatBoundedValidator(0.0), + doc='Can outer radius') + + self.declareProperty(name='BeamHeight', defaultValue=3.0, + validator=FloatBoundedValidator(0.0), + doc='Beam height') + self.declareProperty(name='BeamWidth', defaultValue=2.0, + validator=FloatBoundedValidator(0.0), + doc='Beam width') + + self.declareProperty(name='StepSize', defaultValue=0.002, + validator=FloatBoundedValidator(0.0), + doc='Step size') + self.declareProperty(name='Interpolate', defaultValue=True, + doc='Interpolate the correction workspaces to match the sample workspace') + self.declareProperty(name='NumberWavelengths', defaultValue=10, + validator=IntBoundedValidator(1), + doc='Number of wavelengths for calculation') + + self.declareProperty(name='Emode', defaultValue='Elastic', + validator=StringListValidator(['Elastic', 'Indirect']), + doc='Emode: Elastic or Indirect') + self.declareProperty(name='Efixed', defaultValue=1.0, + doc='Analyser energy') + + self.declareProperty(WorkspaceGroupProperty('OutputWorkspace', '', + direction=Direction.Output), + doc='The output corrections workspace group') +#------------------------------------------------------------------------------ + + def validateInputs(self): + self._setup() + issues = dict() + + # Ensure that a can chemical formula is given when using a can workspace + if self._use_can: + can_chemical_formula = self.getPropertyValue('CanChemicalFormula') + if can_chemical_formula == '': + issues['CanChemicalFormula'] = 'Must provide a chemical foruma when providing a can workspace' + + # Ensure there are enough steps + number_steps = int((self._sample_outer_radius - self._sample_inner_radius) / self._step_size) + if number_steps < 20: + issues['StepSize'] = 'Number of steps (%d) should be >= 20' % number_steps + + return issues + +#------------------------------------------------------------------------------ + + def PyExec(self): + self._setup() + self._sample() + self._wave_range() + self._get_angles() + self._transmission() + + dataA1 = [] + dataA2 = [] + dataA3 = [] + dataA4 = [] + + for angle in self._angles: + (A1, A2, A3, A4) = self._cyl_abs(angle) + logger.information('Angle : %f * successful' % (angle)) + + dataA1 = np.append(dataA1, A1) + dataA2 = np.append(dataA2, A2) + dataA3 = np.append(dataA3, A3) + dataA4 = np.append(dataA4, A4) + + dataX = self._waves * len(self._angles) + + # Create the output workspaces + ass_ws = self._output_ws_name + '_ass' + CreateWorkspace(OutputWorkspace=ass_ws, DataX=dataX, DataY=dataA1, + NSpec=len(self._angles), UnitX='Wavelength') + workspaces = [ass_ws] + + if self._use_can: + assc_ws = self._output_ws_name + '_assc' + workspaces.append(assc_ws) + CreateWorkspace(OutputWorkspace=assc_ws, DataX=dataX, DataY=dataA2, + NSpec=len(self._angles), UnitX='Wavelength') + + acsc_ws = self._output_ws_name + '_acsc' + workspaces.append(acsc_ws) + CreateWorkspace(OutputWorkspace=acsc_ws, DataX=dataX, DataY=dataA3, + NSpec=len(self._angles), UnitX='Wavelength') + + acc_ws = self._output_ws_name + '_acc' + workspaces.append(acc_ws) + CreateWorkspace(OutputWorkspace=acc_ws, DataX=dataX, DataY=dataA4, + NSpec=len(self._angles), UnitX='Wavelength') + + if self._interpolate: + self._interpolate_corrections(workspaces) + try: + self. _copy_detector_table(workspaces) + except RuntimeError: + logger.warning('Cannot copy spectra mapping. Check input workspace instrument.') + + sample_log_workspaces = workspaces + sample_logs = [('sample_shape', 'cylinder'), + ('sample_filename', self._sample_ws_name), + ('sample_inner', self._sample_inner_radius), + ('sample_outer', self._sample_outer_radius)] + + if self._use_can: + sample_logs.append(('can_filename', self._can_ws_name)) + sample_logs.append(('can_outer', self._can_outer_radius)) + + log_names = [item[0] for item in sample_logs] + log_values = [item[1] for item in sample_logs] + + for ws_name in sample_log_workspaces: + AddSampleLogMultiple(Workspace=ws_name, LogNames=log_names, LogValues=log_values) + + GroupWorkspaces(InputWorkspaces=','.join(workspaces), OutputWorkspace=self._output_ws_name) + self.setPropertyValue('OutputWorkspace', self._output_ws_name) + +#------------------------------------------------------------------------------ + + def _setup(self): + self._sample_ws_name = self.getPropertyValue('SampleWorkspace') + self._sample_chemical_formula = self.getPropertyValue('SampleChemicalFormula') + self._sample_number_density = self.getProperty('SampleNumberDensity').value + self._sample_inner_radius = self.getProperty('SampleInnerRadius').value + self._sample_outer_radius = self.getProperty('SampleOuterRadius').value + self._number_can = 1 + + self._can_ws_name = self.getPropertyValue('CanWorkspace') + self._use_can = self._can_ws_name != '' + self._can_chemical_formula = self.getPropertyValue('CanChemicalFormula') + self._can_number_density = self.getProperty('CanNumberDensity').value + self._can_outer_radius = self.getProperty('CanOuterRadius').value + if self._use_can: + self._number_can = 2 + + self._step_size = self.getProperty('StepSize').value + self._radii = np.zeros(self._number_can +1) + self._radii[0] = self._sample_inner_radius + self._radii[1] = self._sample_outer_radius + if (self._radii[1] - self._radii[0]) < 1e-4: + raise ValueError('Sample outer radius not > inner radius') + else: + logger.information('Sample : inner radius = %f ; outer radius = %f' % ( + self._radii[0], self._radii[1])) + self._ms = int((self._radii[1] - self._radii[0] + 0.0001)/self._step_size) + if self._ms < 20: + raise ValueError('Number of steps ( %i ) should be >= 20' % (self._ms)) + else: + if self._ms < 1: + self._ms=1 + logger.information('Sample : ms = %i ' % (self._ms)) + if self._use_can: + self._radii[2] = self._can_outer_radius + if (self._radii[2] - self._radii[1]) < 1e-4: + raise ValueError('Can outer radius not > sample outer radius') + else: + logger.information('Can : inner radius = %f ; outer radius = %f' % ( + self._radii[1], self._radii[2])) + + beam_width = self.getProperty('BeamWidth').value + beam_height = self.getProperty('BeamHeight').value + 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] + + self._interpolate = self.getProperty('Interpolate').value + self._number_wavelengths = self.getProperty('NumberWavelengths').value + + self._emode = self.getPropertyValue('Emode') + self._efixed = self.getProperty('Efixed').value + + self._output_ws_name = self.getPropertyValue('OutputWorkspace') + +#------------------------------------------------------------------------------ + + def _sample(self): + SetSampleMaterial(self._sample_ws_name , ChemicalFormula=self._sample_chemical_formula, + SampleNumberDensity=self._sample_number_density) + sample = mtd[self._sample_ws_name].sample() + sam_material = sample.getMaterial() + # total scattering x-section + self._sig_s = np.zeros(self._number_can) + self._sig_s[0] = sam_material.totalScatterXSection() + # absorption x-section + self._sig_a = np.zeros(self._number_can) + self._sig_a[0] = sam_material.absorbXSection() + # density + self._density = np.zeros(self._number_can) + self._density[0] = self._sample_number_density + + if self._use_can: + SetSampleMaterial(InputWorkspace=self._can_ws_name, ChemicalFormula=self._can_chemical_formula, + SampleNumberDensity=self._can_number_density) + can_sample = mtd[self._can_ws_name].sample() + can_material = can_sample.getMaterial() + self._sig_s[1] = can_material.totalScatterXSection() + self._sig_a[1] = can_material.absorbXSection() + self._density[1] = self._can_number_density + +#------------------------------------------------------------------------------ + + def _get_angles(self): + num_hist = mtd[self._sample_ws_name].getNumberHistograms() + source_pos = mtd[self._sample_ws_name].getInstrument().getSource().getPos() + sample_pos = mtd[self._sample_ws_name].getInstrument().getSample().getPos() + beam_pos = sample_pos - source_pos + self._angles = list() + for index in range(0, num_hist): + detector = mtd[self._sample_ws_name].getDetector(index) + two_theta = detector.getTwoTheta(sample_pos, beam_pos) * 180.0 / math.pi + self._angles.append(two_theta) + logger.information('Detector angles : %i from %f to %f ' % ( + len(self._angles), self._angles[0], self._angles[-1])) + +#------------------------------------------------------------------------------ + + def _wave_range(self): + wave_range = '__wave_range' + ExtractSingleSpectrum(InputWorkspace=self._sample_ws_name, OutputWorkspace=wave_range, WorkspaceIndex=0) + + Xin = mtd[wave_range].readX(0) + wave_min = mtd[wave_range].readX(0)[0] + wave_max = mtd[wave_range].readX(0)[len(Xin) - 1] + number_waves = self._number_wavelengths + wave_bin = (wave_max - wave_min) / (number_waves-1) + + self._waves = list() + for idx in range(0, number_waves): + self._waves.append(wave_min + idx * wave_bin) + DeleteWorkspace(wave_range) + + if self._emode == 'Elastic': + self._elastic = self._waves[int(len(self._waves) / 2)] + elif self._emode == 'Direct': + self._elastic = math.sqrt(81.787/self._efixed) # elastic wavelength + elif self._emode == 'Indirect': + self._elastic = math.sqrt(81.787/self._efixed) # elastic wavelength + + logger.information('Elastic lambda : %f' % (self._elastic)) + logger.information('Lambda : %i values from %f to %f' % ( + len(self._waves), self._waves[0], self._waves[-1])) + +#------------------------------------------------------------------------------ + + def _transmission(self): + distance = self._radii[1] - self._radii[0] + trans= math.exp(-distance*self._density[0]*(self._sig_s[0] + self._sig_a[0])) + logger.information('Sample transmission : %f' % (trans)) + if self._use_can: + distance = self._radii[2] - self._radii[1] + trans= math.exp(-distance*self._density[1]*(self._sig_s[1] + self._sig_a[1])) + logger.information('Can transmission : %f' % (trans)) + +#------------------------------------------------------------------------------ + + def _interpolate_corrections(self, workspaces): + """ + Performs interpolation on the correction workspaces such that the number of bins + matches that of the input sample workspace. + + @param workspaces List of correction workspaces to interpolate + """ + + for ws in workspaces: + SplineInterpolation(WorkspaceToMatch=self._sample_ws_name, + WorkspaceToInterpolate=ws, + OutputWorkspace=ws, + OutputWorkspaceDeriv='') + +#------------------------------------------------------------------------------ + + def _copy_detector_table(self, workspaces): + """ + Copy the detector table from the sample workspaces to the correction workspaces. + + @param workspaces List of correction workspaces + """ + + instrument = mtd[self._sample_ws_name].getInstrument().getName() + + for ws in workspaces: + LoadInstrument(Workspace=ws, + InstrumentName=instrument) + + CopyDetectorMapping(WorkspaceToMatch=self._sample_ws_name, + WorkspaceToRemap=ws, + IndexBySpectrumNumber=True) + +#------------------------------------------------------------------------------ + + def _cyl_abs(self, angle): +# Parameters : +# self._step_size - step size +# self._beam - beam parameters +# nan - number of annuli +# radii - list of radii (for each annulus) +# density - list of densities (for each annulus) +# sigs - list of scattering cross-sections (for each annulus) +# siga - list of absorption cross-sections (for each annulus) +# angle - list of angles +# wavelas - elastic wavelength +# waves - list of wavelengths +# Output parameters : A1 - Ass ; A2 - Assc ; A3 - Acsc ; A4 - Acc + + amu_scat = np.zeros(self._number_can) + amu_scat = self._density*self._sig_s + sig_abs = np.zeros(self._number_can) + sig_abs = self._density*self._sig_a + amu_tot_i = np.zeros(self._number_can) + amu_tot_s = np.zeros(self._number_can) + + theta = angle*math.pi/180. + A1 = [] + A2 = [] + A3 = [] + A4 = [] + #loop over wavelengths + for wave in self._waves: + #loop over annuli + if self._emode == 'Elastic': + amu_tot_i = amu_scat + sig_abs*self._elastic/1.7979 + amu_tot_s = amu_scat + sig_abs*self._elastic/1.7979 + if self._emode == 'Direct': + amu_tot_i = amu_scat + sig_abs*self._elastic/1.7979 + amu_tot_s = amu_scat + sig_abs*wave/1.7979 + if self._emode == 'Indirect': + amu_tot_i = amu_scat + sig_abs*wave/1.7979 + amu_tot_s = amu_scat + sig_abs*self._elastic/1.7979 + (Ass, Assc, Acsc, Acc) = self._acyl(theta, amu_scat, amu_tot_i, amu_tot_s) + A1.append(Ass) + A2.append(Assc) + A3.append(Acsc) + A4.append(Acc) + return A1, A2, A3, A4 + +#------------------------------------------------------------------------------ + + def _acyl(self, theta, amu_scat, amu_tot_i, amu_tot_s): + A = self._beam[1] + Area_s = 0.0 + Ass = 0.0 + Acc = 0.0 + Acsc = 0.0 + Assc = 0.0 + nan = self._number_can + if self._number_can < 2: +# +# No. STEPS ARE CHOSEN SO THAT STEP WIDTH IS THE SAME FOR ALL ANNULI +# + AAAA, BBBA, Area_A = self._sum_rom(0, 0, A, self._radii[0], self._radii[1], self._ms, + theta, amu_scat, amu_tot_i, amu_tot_s) + AAAB, BBBB, Area_B = self._sum_rom(0, 0, -A, self._radii[0], self._radii[1], self._ms, + theta, amu_scat, amu_tot_i, amu_tot_s) + Area_s += Area_A + Area_B + Ass += AAAA + AAAB + Ass = Ass/Area_s + else: + for i in range(0, self._number_can -1): + radius_1 = self._radii[i] + radius_2 = self._radii[i+1] +# +# No. STEPS ARE CHOSEN SO THAT STEP WIDTH IS THE SAME FOR ALL ANNULI +# + ms = int(self._ms*(radius_2 - radius_1)/(self._radii[1] - self._radii[0])) + if ms < 1: + ms = 1 + AAAA, BBBA, Area_A = self._sum_rom(i, 0, A, radius_1, radius_2, + ms, theta, amu_scat, amu_tot_i, amu_tot_s) + AAAB, BBBB, Area_B = self._sum_rom(i, 0, -A, radius_1, radius_2, + ms, theta, amu_scat, amu_tot_i, amu_tot_s) + Area_s += Area_A + Area_B + Ass += AAAA + AAAB + Assc += BBBA + BBBB + Ass = Ass/Area_s + Assc = Assc/Area_s + radius_1 = self._radii[nan -1] + radius_2 = self._radii[nan] + ms = int(self._ms*(radius_2 - radius_1)/(self._radii[1] - self._radii[0])) + if ms < 1: + ms = 1 + AAAA, BBBA, Area_A = self._sum_rom(nan-1, 1, A, radius_1, radius_2, + ms, theta, amu_scat, amu_tot_i, amu_tot_s) + AAAB, BBBB, Area_B = self._sum_rom(nan-1, 1, -A, radius_1, radius_2, + ms, theta, amu_scat, amu_tot_i, amu_tot_s) + Area_C = Area_A + Area_B + Acsc = (AAAA + AAAB)/Area_C + Acc = (BBBA + BBBB)/Area_C + return Ass, Assc, Acsc, Acc + +#------------------------------------------------------------------------------ + + def _sum_rom(self, n_scat, n_abs, a, r1, r2, ms, theta, amu_scat, amu_tot_i, amu_tot_s): + #n_scat is region for scattering + #n_abs is region for absorption + nan = self._number_can + omega_add = 0. + if a < 0.: + omega_add = math.pi + AAA = 0. + BBB = 0. + Area = 0. + theta_deg = math.pi - theta + num = ms + r_step = (r2 - r1)/ms + r_add = -0.5*r_step + r1 + +# start loop over M + for M in range(1, num+1): + r = M*r_step + r_add + number_omega = int(math.pi*r/r_step) + omega_ster = math.pi/number_omega + omega_deg = -0.5*omega_ster + omega_add + Area_y = r*r_step*omega_ster*amu_scat[n_scat] + sum_1 = 0. + sum_2 = 0. + I = 1 + Area_sum = 0. + for _ in range(1, number_omega +1): + omega = I*omega_ster + omega_deg + distance = r*math.sin(omega) + + skip = True + if abs(distance) <= a: +# +# CALCULATE DISTANCE INCIDENT NEUTRON PASSES THROUGH EACH ANNULUS + LIS = [] + for j in range(0, nan): + LIST = self._distance(r, self._radii[j], omega) + LISN = self._distance(r, self._radii[j+1], omega) + LIS.append(LISN - LIST) +# +# CALCULATE DISTANCE SCATTERED NEUTRON PASSES THROUGH EACH ANNULUS + O = omega + theta_deg + LSS = [] + for j in range(0, nan): + LSST = self._distance(r, self._radii[j], O) + LSSN = self._distance(r, self._radii[j+1], O) + LSS.append(LSSN - LSST) +# +# CALCULATE ABSORBTION FOR PATH THROUGH ALL ANNULI,AND THROUGH INNER ANNULI + path = np.zeros(3) +# split into input (I) and scattered (S) paths + path[0] += amu_tot_i[0]*LIS[0] + amu_tot_s[0]*LSS[0] + if nan == 2: + path[2] += amu_tot_i[1]*LIS[1] + amu_tot_s[1]*LSS[1] + path[1] = path[0] + path[2] + sum_1 += math.exp(-path[n_abs]) + sum_2 += math.exp(-path[n_abs +1]) + Area_sum += 1.0 + skip = False + + if skip: + I = number_omega -I +2 + else: + I += 1 + AAA += sum_1*Area_y + BBB += sum_2*Area_y + Area += Area_sum*Area_y + return AAA, BBB, Area + +#------------------------------------------------------------------------------ + + def _distance(self, r1, radius, omega): + r = r1 + distance = 0. + b = r*math.sin(omega) + if abs(b) < radius: + t = r*math.cos(omega) + c = radius*radius -b*b + d = math.sqrt(c) + if r <= radius: + distance = t + d + else: + distance = d*(1.0 + math.copysign(1.0,t)) + return distance + +#------------------------------------------------------------------------------ + +# Register algorithm with Mantid +AlgorithmFactory.subscribe(CylinderPaalmanPingsCorrection) diff --git a/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt b/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt index f3ea6367c4ea..9e32c9b72529 100644 --- a/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt +++ b/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt @@ -14,6 +14,7 @@ set ( TEST_PY_FILES IndirectCalibrationTest.py CreateWorkspaceTest.py CylinderPaalmanPingsCorrectionTest.py + CylinderPaalmanPingsCorrection2Test.py DakotaChiSquaredTest.py DensityOfStatesTest.py DSFinterpTest.py diff --git a/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/CylinderPaalmanPingsCorrection2Test.py b/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/CylinderPaalmanPingsCorrection2Test.py new file mode 100644 index 000000000000..5e1a38d02496 --- /dev/null +++ b/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/CylinderPaalmanPingsCorrection2Test.py @@ -0,0 +1,177 @@ +import unittest +from mantid.kernel import * +from mantid.api import * +from mantid.simpleapi import (CreateSampleWorkspace, Scale, DeleteWorkspace, + CylinderPaalmanPingsCorrection) + + +class CylinderPaalmanPingsCorrection2Test(unittest.TestCase): + + def setUp(self): + """ + Create sample workspaces. + """ + + # Create some test data + sample = CreateSampleWorkspace(NumBanks=1, + BankPixelWidth=1, + XUnit='Wavelength', + XMin=6.8, + XMax=7.9, + BinWidth=0.1) + self._sample_ws = sample + + can = Scale(InputWorkspace=sample, Factor=1.2) + self._can_ws = can + + self._corrections_ws_name = 'corrections' + + + def tearDown(self): + """ + Remove workspaces from ADS. + """ + + DeleteWorkspace(self._sample_ws) + DeleteWorkspace(self._can_ws) + + if self._corrections_ws_name in mtd: + DeleteWorkspace(self._corrections_ws_name) + + + def _verify_workspace(self, ws_name): + """ + Do validation on a correction workspace. + + @param ws_name Name of workspace to validate + """ + + corrections_ws = mtd[self._corrections_ws_name] + + # Check it is in the corrections workspace group + self.assertTrue(corrections_ws.contains(ws_name)) + + test_ws = mtd[ws_name] + + # Check workspace is in wavelength + self.assertEqual(test_ws.getAxis(0).getUnit().unitID(), + 'Wavelength') + + # Check it has the same number of spectra as the sample + self.assertEqual(test_ws.getNumberHistograms(), + self._sample_ws.getNumberHistograms()) + + # Check it has X binning matching sample workspace + self.assertEqual(test_ws.blocksize(), self._sample_ws.blocksize()) + + + def _verify_workspaces_for_can(self): + """ + Do validation on the additional correction factors for sample and can. + """ + + ass_ws_name = self._corrections_ws_name + '_ass' + assc_ws_name = self._corrections_ws_name + '_assc' + acsc_ws_name = self._corrections_ws_name + '_acsc' + acc_ws_name = self._corrections_ws_name + '_acc' + + workspaces = [ass_ws_name, assc_ws_name, acsc_ws_name, acc_ws_name] + + for workspace in workspaces: + self._verify_workspace(workspace) + + + def test_sampleOnly(self): + """ + Test simple run with sample workspace only. + """ + + CylinderPaalmanPingsCorrection(OutputWorkspace=self._corrections_ws_name, + SampleWorkspace=self._sample_ws, + SampleChemicalFormula='H2-O', + SampleInnerRadius=0.05, + SampleOuterRadius=0.1, + Emode='Indirect', + Efixed=1.845) + + ass_ws_name = self._corrections_ws_name + '_ass' + self. _verify_workspace(ass_ws_name) + + + def test_sampleAndCan(self): + """ + Test simple run with sample and can workspace. + """ + + CylinderPaalmanPingsCorrection(OutputWorkspace=self._corrections_ws_name, + SampleWorkspace=self._sample_ws, + SampleChemicalFormula='H2-O', + SampleInnerRadius=0.05, + SampleOuterRadius=0.1, + CanWorkspace=self._can_ws, + CanChemicalFormula='V', + CanOuterRadius=0.15, + BeamHeight=0.1, + BeamWidth=0.1, + Emode='Indirect', + Efixed=1.845) + + self._verify_workspaces_for_can() + + + def test_sampleAndCanDefaults(self): + """ + Test simple run with sample and can workspace using the default values. + """ + + CylinderPaalmanPingsCorrection(OutputWorkspace=self._corrections_ws_name, + SampleWorkspace=self._sample_ws, + SampleChemicalFormula='H2-O', + CanWorkspace=self._can_ws, + CanChemicalFormula='V') + + self._verify_workspaces_for_can() + + + def test_InterpolateDisabled(self): + """ + Tests that a workspace with a bin count equal to NumberWavelengths is created + when interpolation is disabled. + """ + + CylinderPaalmanPingsCorrection(OutputWorkspace=self._corrections_ws_name, + SampleWorkspace=self._sample_ws, + SampleChemicalFormula='H2-O', + CanWorkspace=self._can_ws, + CanChemicalFormula='V', + Interpolate=False) + + corrections_ws = mtd[self._corrections_ws_name] + + # Check each correction workspace has X binning matching NumberWavelengths + for workspace in corrections_ws: + self.assertEqual(workspace.blocksize(), 10) + + + def test_validationNoCanFormula(self): + """ + Tests validation for no chemical formula for can when a can WS is provided. + """ + + self.assertRaises(RuntimeError, + CylinderPaalmanPingsCorrection, + OutputWorkspace=self._corrections_ws_name, + SampleWorkspace=self._sample_ws, + SampleChemicalFormula='H2-O', + SampleInnerRadius=0.05, + SampleOuterRadius=0.1, + CanWorkspace=self._can_ws, + CanOuterRadius=0.15, + BeamHeight=0.1, + BeamWidth=0.1, + Emode='Indirect', + Efixed=1.845) + + +if __name__=="__main__": + unittest.main() diff --git a/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/CylinderPaalmanPingsCorrectionTest.py b/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/CylinderPaalmanPingsCorrectionTest.py index 27955836344d..a1f760636a72 100644 --- a/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/CylinderPaalmanPingsCorrectionTest.py +++ b/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/CylinderPaalmanPingsCorrectionTest.py @@ -14,11 +14,11 @@ def setUp(self): # Create some test data sample = CreateSampleWorkspace(NumBanks=1, - BankPixelWidth=1, - XUnit='Wavelength', - XMin=6.8, - XMax=7.9, - BinWidth=0.1) + BankPixelWidth=1, + XUnit='Wavelength', + XMin=6.8, + XMax=7.9, + BinWidth=0.1) self._sample_ws = sample can = Scale(InputWorkspace=sample, Factor=1.2) @@ -96,7 +96,8 @@ def test_sampleOnly(self): SampleInnerRadius=0.05, SampleOuterRadius=0.1, Emode='Indirect', - Efixed=1.845) + Efixed=1.845, + Version=1) ass_ws_name = self._corrections_ws_name + '_ass' self. _verify_workspace(ass_ws_name) @@ -122,7 +123,8 @@ def test_sampleAndCan(self): BeamHeight=0.1, BeamWidth=0.1, Emode='Indirect', - Efixed=1.845) + Efixed=1.845, + Version=1) self._verify_workspaces_for_can() @@ -140,7 +142,8 @@ def test_sampleAndCanDefaults(self): SampleWorkspace=self._sample_ws, SampleChemicalFormula='H2-O', CanWorkspace=self._can_ws, - CanChemicalFormula='V') + CanChemicalFormula='V', + Version=1) self._verify_workspaces_for_can() @@ -160,7 +163,8 @@ def test_InterpolateDisabled(self): SampleChemicalFormula='H2-O', CanWorkspace=self._can_ws, CanChemicalFormula='V', - Interpolate=False) + Interpolate=False, + Version=1) corrections_ws = mtd[self._corrections_ws_name] @@ -186,7 +190,8 @@ def test_validationNoCanFormula(self): BeamHeight=0.1, BeamWidth=0.1, Emode='Indirect', - Efixed=1.845) + Efixed=1.845, + Version=1) if __name__=="__main__": diff --git a/Code/Mantid/docs/source/algorithms/CylinderPaalmanPingsCorrection-v1.rst b/Code/Mantid/docs/source/algorithms/CylinderPaalmanPingsCorrection-v1.rst index b1e4a48ccfef..7786c6f41db4 100644 --- a/Code/Mantid/docs/source/algorithms/CylinderPaalmanPingsCorrection-v1.rst +++ b/Code/Mantid/docs/source/algorithms/CylinderPaalmanPingsCorrection-v1.rst @@ -51,7 +51,8 @@ Usage BeamWidth=0.1, StepSize=0.002, Emode='Indirect', - Efixed=1.845) + Efixed=1.845, + Version=1) print 'Correction workspaces: %s' % (', '.join(corr.getNames())) diff --git a/Code/Mantid/docs/source/algorithms/CylinderPaalmanPingsCorrection-v2.rst b/Code/Mantid/docs/source/algorithms/CylinderPaalmanPingsCorrection-v2.rst new file mode 100644 index 000000000000..6c3eac156f0c --- /dev/null +++ b/Code/Mantid/docs/source/algorithms/CylinderPaalmanPingsCorrection-v2.rst @@ -0,0 +1,64 @@ +.. algorithm:: + +.. summary:: + +.. alias:: + +.. properties:: + +Description +----------- + +Calculates absorption corrections for a cylindrical or annular sample giving +output in the Paalman & Pings absorption factors: :math:`A_{s,s}` (correction +factor for scattering and absorption in sample), :math:`A_{s,sc}` (scattering in +sample and absorption in sample and container), :math:`A_{c,sc}` (scattering in +container and absorption in sample and container) and :math:`A_{c,c}` +(scattering and absorption in container). + +Restrictions on the input workspace +################################### + +The input workspace must have a fully defined instrument that has X axis units +of wavelength. + +Usage +----- + +**Example:** + +.. testcode:: ExCylinderPaalmanPingsCorrection + + # Create a sample workspace + sample = CreateSampleWorkspace(NumBanks=1, BankPixelWidth=1, + XUnit='Wavelength', + XMin=6.8, XMax=7.9, + BinWidth=0.1) + + # Copy and scale it to make a can workspace + can = CloneWorkspace(InputWorkspace=sample) + can = Scale(InputWorkspace=can, Factor=1.2) + + # Calculate absorption corrections + corr = CylinderPaalmanPingsCorrection(SampleWorkspace=sample, + SampleChemicalFormula='H2-O', + SampleInnerRadius=0.05, + SampleOuterRadius=0.1, + CanWorkspace=can, + CanChemicalFormula='V', + CanOuterRadius=0.15, + BeamHeight=0.1, + BeamWidth=0.1, + StepSize=0.002, + Emode='Indirect', + Efixed=1.845) + + print 'Correction workspaces: %s' % (', '.join(corr.getNames())) + +Output: + +.. testoutput:: ExCylinderPaalmanPingsCorrection + + Correction workspaces: corr_ass, corr_assc, corr_acsc, corr_acc + +.. categories::